About
This script simulates the uniaxial tensile deformation of a rectangular 3D bar using the finite element method (FEM) with hexahedral elements in scikit-fem.
It creates a structured mesh of the bar, applies a fixed boundary condition on one end and a uniform traction load on the opposite end, then solves for the displacement field under linear elasticity.
Finally, it compares the numerical displacement from FEM with the analytical solution derived from Hooke’s law.
Tensile Test
\delta = \varepsilon , L = \frac{\sigma}{E}\\L = \frac{F L}{E A}In short, this function computes the basic linear elastic tensile displacement from the Young’s modulus (E), cross-sectional area (A), length (L), and applied force (F).
Implementation
import numpy as np
import scipy.sparse.linalg
import scipy
import skfem
from skfem import *
from skfem.helpers import ddot, sym_grad, eye, trace
from skfem.models.elasticity import lame_parameters
from skfem.models.elasticity import linear_elasticity
from skfem.helpers import dot, sym_grad
# Parameters
E = 210e3 # Young [Pa]
nu = 0.3 # Poisson
x_len, y_len, z_len = 10.0, 10.0, 50.0 # The size of Box
mesh_size = 1.0
displacement_scale = 1
# mesh = skfem.MeshHex().refined(3).with_defaults()
mesh = create_box_hex(x_len, y_len, z_len, mesh_size)
e = skfem.ElementVector(skfem.ElementHex1())
basis = skfem.Basis(mesh, e, intorder=1)
lam, mu = lame_parameters(E, nu)
fixed_dofs = basis.get_dofs(
mesh.facets_satisfying(lambda x: np.isclose(x[2], 0.0))
).nodal['u^3']
F_facet_ids = mesh.facets_satisfying(
lambda x: np.isclose(x[2], z_len)
)
fbasis = skfem.FacetBasis(mesh, e, facets=F_facet_ids)
total_force = 100.0 # [N]
@skfem.Functional
def surface_measure(w):
return 1.0
A_proj_z = surface_measure.assemble(fbasis)
pressure = total_force / A_proj_z
# ---- traction: t = p * n （l(v) = ∫_Γ (p n) · v ds）----
# @skfem.LinearForm
# def l_comp(v, w):
# return pressure * v[2]
@skfem.LinearForm
def l_comp(v, w):
# dim, local_node, quad_point
# print(v.shape) # (3, 4, 16)
return pressure * dot(w.n, v)
K = linear_elasticity(lam, mu).assemble(basis)
F = l_comp.assemble(fbasis)
K_c, F_c, U_c, I = skfem.condense(K, F, D=fixed_dofs)
u = skfem.solve(K_c, F_c, U_c, I)
reaction_force = K[fixed_dofs, :] @ u - F[fixed_dofs]
print("reaction_force:", np.sum(reaction_force))
# print(np.sum(U_e - U_c))
# ---------------------------
# print(U, np.max(U), np.min(U))
# m = mesh.translated(displacement_scale * U_e[basis.nodal_dofs])
# m.save('box_deformed.vtk')
# mesh.save('box_org.vtk')
def tensile_displacement(E, nu, x_len, y_len, z_len, F):
A = x_len * y_len
stress = F / A
strain = stress / E
displacement = strain * z_len
return displacement, stress, strain
displacement, stress, strain = tensile_displacement(
E, nu, x_len, y_len, z_len, total_force
)
print(f"maximum displacement-fem: {np.max(u)}")
print(f"displacement: {displacement}")
Tools
from typing import Optional
import numpy as np
import skfem
import meshio
def create_box_hex(
x_len: int, y_len: int, z_len: int, mesh_size: int
):
"""
Create a hexahedral mesh box with given dimensions and element size.
Parameters
----------
x_len, y_len, z_len : float
Dimensions of the box in x, y, z directions.
mesh_size : float
Desired approximate size of each hex element.
Returns
-------
mesh : MeshHex
A scaled MeshHex object with the specified dimensions.
"""
nx = int(np.ceil(x_len / mesh_size))
ny = int(np.ceil(y_len / mesh_size))
nz = int(np.ceil(z_len / mesh_size))
x = np.linspace(0, x_len, nx + 1)
y = np.linspace(0, y_len, ny + 1)
z = np.linspace(0, z_len, nz + 1)
mesh = skfem.MeshHex.init_tensor(x, y, z)
mesh_fixed = skfem.MeshHex(mesh.p, mesh.t)
return mesh_fixed
def export_mesh_with_info(
mesh: skfem.Mesh,
point_data_values: Optional[list[np.ndarray]] = None,
point_data_names: Optional[list[str]] = None,
cell_data_values: Optional[list[np.ndarray]] = None,
cell_data_names: Optional[list[str]] = None,
filepath: str = "output.vtu"
):
"""
Export a skfem.Mesh object and its data to a VTU file via meshio.
Parameters
----------
mesh : skfem.Mesh
The finite element mesh object (MeshTet, MeshTri, MeshHex, etc.).
point_data_values : list of np.ndarray, optional
List of arrays of point-wise data (length = n_nodes).
point_data_names : list of str, optional
Names for each point-wise data array.
cell_data_values : list of np.ndarray, optional
List of arrays of cell-wise data (length = n_elements).
cell_data_names : list of str, optional
Names for each cell-wise data array.
filepath : str
Output filename (e.g., "result.vtu").
"""
# Determine element type
if isinstance(mesh, skfem.MeshTet):
cell_type = "tetra"
elif isinstance(mesh, skfem.MeshHex):
cell_type = "hexahedron"
else:
raise ValueError(f"Unsupported mesh type: {type(mesh)}")
# Convert skfem's (dim, n) shape to meshio format
points = mesh.p.T
# cells = [(cell_type, mesh.t.T)]
# Point data
point_data = {}
if point_data_values and point_data_names:
for name, val in zip(point_data_names, point_data_values):
point_data[name] = val
# Cell data (wrap in a list for each cell block)
cell_data = {}
if cell_data_values and cell_data_names:
for name, val in zip(cell_data_names, cell_data_values):
cell_data[name] = [val]
# Build and write meshio.Mesh
meshio_mesh = meshio.Mesh(
points=points,
cells=[meshio.CellBlock(cell_type, mesh.t.T)],
point_data=point_data,
cell_data=cell_data
)
meshio_mesh.write(filepath)
