Linear Elasticity: Tensile Bar#

This tutorial walks through the minimal tensile-bar example in tutorials/linearelastic_tensile_bar.py and explains the weak form, assembly, and boundary conditions with the key equations.

Run the example#

python tutorials/linearelastic_tensile_bar.py

Problem statement#

We solve a small-strain linear elasticity problem on a 3D bar:

  • Unknown displacement field: u

  • Test function: v

  • Material: isotropic, given by D(E, nu)

  • Boundary conditions: - Dirichlet (clamped) on x = xmin - Uniform traction on x = xmax along +x

Weak form#

Find u such that for all v:

\[\int_{\Omega} \varepsilon(v) : D : \varepsilon(u)\, d\Omega = \int_{\Gamma_t} v \cdot t \, ds\]

where \(\varepsilon(u) = \tfrac{1}{2}(\nabla u + \nabla u^T)\).

In this example, t = (traction, 0, 0) is applied on x = xmax.

Implementation flow (FluxFEM)#

1) Build mesh and space#

The script builds a structured hexahedral mesh and creates a vector-valued FEM space:

mesh = ff.StructuredHexBox(nx=nx, ny=ny, nz=nz, lx=lx, ly=ly, lz=lz).build()
space = ff.make_hex_space(mesh, dim=3, intorder=intorder)

2) Assemble the bilinear form (weak form)#

\[a(u, v) = \int_{\Omega} \varepsilon(v) : D : \varepsilon(u)\, d\Omega\]

Single-space assembly stays shortest for this standard Galerkin problem:

import fluxfem.helpers_wf as h_wf

bilinear_form = ff.BilinearForm.volume(
    lambda u, v, D_mat: h_wf.ddot(v.sym_grad, h_wf.matmul_std(D_mat, u.sym_grad))
    * h_wf.dOmega()
)
K = space.assemble(bilinear_form, params=D)

The same operator can also be assembled through explicit test/trial roles:

U = ff.NamedSpace("U", space)
V = ff.NamedSpace("V", space)

bilinear = ff.BilinearForm.volume(
    lambda u, v, D_mat: h_wf.ddot(v.sym_grad, h_wf.matmul_std(D_mat, u.sym_grad))
    * h_wf.dOmega()
)
K = ff.assemble_bilinear_form(
    ff.BilinearSpaces(test=V, trial=U),
    bilinear,
    D,
)

3) Assemble the surface traction#

\[\ell(v) = \int_{\Gamma_t} v \cdot t \, ds\]
surface_form = ff.LinearForm.surface(lambda v, p: (v | p) * h_wf.ds())
traction_vec = np.array([traction, 0.0, 0.0], dtype=float)
F = surface.assemble_linear_form_on_space(space, surface_form, params=traction_vec)

For volume load vectors, the matching role-explicit entry point is ff.assemble_linear_form(ff.LinearSpaces(test=V), ...).

4) Apply Dirichlet clamp#

The clamp fixes all components on the x = xmin face:

dir_dofs = ff.DirichletBC.from_boundary_dofs(
    mesh,
    lambda pts: np.isclose(pts[:, 0], xmin, atol=1e-8),
    components="xyz",
).dofs

5) Solve the linear system#

solver = ff.LinearSolver(method="spsolve")
u, _ = solver.solve(
    K, F, dirichlet=ff.DirichletBC(dir_dofs, None), dirichlet_mode="condense"
)

The script also prints the maximum axial displacement at x = xmax and compares it with the 1D bar theory u_x(L) = traction * L / E.