Complete Examples (All Libraries)

This page gives practical examples for every SILEXlight library:

SILEXlight.silex_lib_gmsh
SILEXlight.silex_lib_tri3_python
SILEXlight.silex_lib_tet4_python
SILEXlight.silex_lib_tri3_fortran
SILEXlight.silex_lib_tet4_fortran

Conventions used in all examples

Node numbering in element connectivity starts at 1.
nodes shapes: - 2D: (n_nodes, 2) - 3D: (n_nodes, 3)
elements shapes: - TRI3: (n_elem, 3) - TET4: (n_elem, 4)
Material arrays: - TRI3 plane stress: [Young, nu, thickness] - TET4 3D isotropic: [Young, nu]

Gmsh IO Library

Write a minimal TRI3 mesh and nodal displacement field to .msh:

import numpy as np
from SILEXlight import silex_lib_gmsh as gmsh

nodes = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 1.0],
])
elements = np.array([[1, 2, 3]], dtype=int)

# Example nodal displacement Ux, Uy for each node
U = np.array([
    [0.0, 0.0],
    [1e-5, 0.0],
    [0.0, 1e-5],
])

fields = [
    [U, 'nodal', 2, 'Displacement'],
]

gmsh.WriteResults('tri3_demo', nodes, elements, 2, fields)

Read back nodes and elements from a gmsh file:

import numpy as np
from SILEXlight import silex_lib_gmsh as gmsh

nodes2d = gmsh.ReadGmshNodes('tri3_demo.msh', nbcoord=2)
tri3, used_node_ids = gmsh.ReadGmshElements('tri3_demo.msh', elmttype=2, prop=1)

print(nodes2d.shape)       # (n_nodes, 2)
print(tri3.shape)          # (n_triangles, 3)
print(used_node_ids[:10])

TRI3 Python Library

Compute element assembly triplets (I, J, V) for stiffness matrix:

import numpy as np
from SILEXlight import silex_lib_tri3_python as tri3

nodes = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 1.0],
], dtype=float)
elements = np.array([[1, 2, 3]], dtype=int)
material = np.array([210e9, 0.30, 0.01], dtype=float)

Ik, Jk, Vk = tri3.stiffnessmatrix(nodes, elements, material)
print(Ik.shape, Jk.shape, Vk.shape)  # (36,) each for one TRI3 element

Compute stress, strain, and error indicator from a displacement vector:

import numpy as np
from SILEXlight import silex_lib_tri3_python as tri3

nodes = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 1.0],
], dtype=float)
elements = np.array([[1, 2, 3]], dtype=int)
material = np.array([210e9, 0.30, 0.01], dtype=float)

# 2 dof per node -> 6 dof total
QQ = np.array([0.0, 0.0, 1e-5, 0.0, 0.0, 1e-5], dtype=float)

Sigma, sig_smooth, Epsilon, EpsilonNodes, ErrElem, ErrGlob = tri3.compute_stress_strain_error(
    nodes, elements, material, QQ
)

print(Sigma.shape)         # (n_elem, 4)  : [sx, sy, txy, von_mises]
print(sig_smooth.shape)    # (n_nodes, 4)
print(Epsilon.shape)       # (n_elem, 3)
print(EpsilonNodes.shape)  # (n_nodes, 3)
print(float(ErrGlob))

Apply distributed force on line elements:

import numpy as np
from SILEXlight import silex_lib_tri3_python as tri3

nodes = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 1.0],
], dtype=float)

# Load applied on edge (node 1 -> node 2)
edge_elements = np.array([[1, 2]], dtype=int)

# [fx(pt1), fy(pt1), fx(pt2), fy(pt2)] in N/m
fs = np.array([0.0, -1000.0, 0.0, -1000.0], dtype=float)

# Coordinates of pt1 and pt2
pts = np.array([0.0, 0.0, 1.0, 0.0], dtype=float)

F = tri3.forceonline(nodes, edge_elements, fs, pts)
print(F.shape)  # (2*n_nodes,)

TET4 Python Library

Compute 3D stiffness assembly triplets:

import numpy as np
from SILEXlight import silex_lib_tet4_python as tet4

nodes = np.array([
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
], dtype=float)
elements = np.array([[1, 2, 3, 4]], dtype=int)
material = np.array([210e9, 0.30], dtype=float)

Ik, Jk, Vk = tet4.stiffnessmatrix(nodes, elements, material)
print(Ik.shape, Jk.shape, Vk.shape)  # (144,) each for one TET4 element

Compute consistent mass assembly triplets:

import numpy as np
from SILEXlight import silex_lib_tet4_python as tet4

nodes = np.array([
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
], dtype=float)
elements = np.array([[1, 2, 3, 4]], dtype=int)
rho = 7800.0

Ik, Jk, Vk = tet4.massmatrix(nodes, elements, rho)
print(Ik.shape, Jk.shape, Vk.shape)  # (144,) each

Apply pressure (or directional traction) on triangular surface elements:

import numpy as np
from SILEXlight import silex_lib_tet4_python as tet4

nodes = np.array([
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
], dtype=float)

# Surface triangle on nodes (1,2,3)
surface_tri = np.array([[1, 2, 3]], dtype=int)

press = 2.0e6

# direction=[0,0,0] -> pressure normal to face
Fp = tet4.forceonsurface(nodes, surface_tri, press, direction=[0.0, 0.0, 0.0])
print(Fp.shape)  # (3*n_nodes,)

Compute 3D stress/strain and global error estimator:

import numpy as np
from SILEXlight import silex_lib_tet4_python as tet4

nodes = np.array([
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
], dtype=float)
elements = np.array([[1, 2, 3, 4]], dtype=int)
material = np.array([210e9, 0.30], dtype=float)

# 3 dof per node -> 12 dof total
QQ = np.array([
    0.0, 0.0, 0.0,
    1e-5, 0.0, 0.0,
    0.0, 1e-5, 0.0,
    0.0, 0.0, 1e-5,
], dtype=float)

sigma, sig_smooth, eps_elem, eps_nodes, err_elem, err_glob = tet4.compute_stress_strain_error(
    nodes, elements, material, QQ
)

print(sigma.shape)       # (n_elem, 7)
print(sig_smooth.shape)  # (n_nodes, 7)
print(eps_elem.shape)    # (n_elem, 6)
print(eps_nodes.shape)   # (n_nodes, 6)
print(float(err_glob))

Fortran Libraries (Optional, Compiled)

The Fortran modules provide the same core finite-element routines with a compiled backend:

SILEXlight.silex_lib_tri3_fortran
SILEXlight.silex_lib_tet4_fortran

If the package is installed with the Fortran extension built, usage is analogous.

TRI3 Fortran example:

import numpy as np
from SILEXlight import silex_lib_tri3_fortran as tri3f

nodes = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 1.0],
], dtype=float)
elements = np.array([[1, 2, 3]], dtype=int)
material = np.array([210e9, 0.30, 0.01], dtype=float)

Ik, Jk, Vk = tri3f.stiffnessmatrix(nodes, elements, material)
print(Ik.shape, Jk.shape, Vk.shape)

TET4 Fortran example:

import numpy as np
from SILEXlight import silex_lib_tet4_fortran as tet4f

nodes = np.array([
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
], dtype=float)
elements = np.array([[1, 2, 3, 4]], dtype=int)
material = np.array([210e9, 0.30], dtype=float)

Ik, Jk, Vk = tet4f.stiffnessmatrix(nodes, elements, material)
print(Ik.shape, Jk.shape, Vk.shape)

Notes on Solver Integration

All stiffnessmatrix and massmatrix functions return assembly triplets (Ik, Jk, Vk). A common sparse assembly pattern with SciPy is:

import scipy.sparse as sp

ndof = 2 * n_nodes   # TRI3
# ndof = 3 * n_nodes # TET4

K = sp.coo_matrix((Vk, (Ik, Jk)), shape=(ndof, ndof)).tocsr()

This assembled matrix can then be used with your preferred linear solver.