RasmahRasmah

The finite element method

Overview

So far every shape was described exactly — as a signed distance field, a boundary representation, or a mesh. Now we ask the question CAD models exist to answer: what happens to the part when it is loaded?

The finite element method (FEM) is the standard way to answer that question numerically. It replaces a continuous partial differential equation with a finite linear system $K u = f$, which a computer solves. The idea is simple enough to state in a sentence: divide the solid into small pieces, write the physics on each piece, and assemble the pieces into one big system.

This chapter builds that system from the ground up for the two simplest physics — linear elasticity and steady heat conduction — using the low-level assembler functions. The next chapter, Solving with physics, shows the one-line solve(...) wrapper that hides all of this.

From a differential equation to a linear system

A solid at equilibrium obeys the momentum balance

\[\nabla \cdot \boldsymbol{\sigma} + f = 0 ,\]

where $\boldsymbol{\sigma}$ is the stress tensor and $f$ is the body force. For a linear-elastic, isotropic material the stress is proportional to the strain (Hooke's law), and the strain is the symmetric gradient of the displacement field $u$:

\[\boldsymbol{\sigma} = C\,\varepsilon(u), \qquad \varepsilon(u) = \tfrac{1}{2}\big(\nabla u + \nabla u^{T}\big) .\]

Substituting gives a second-order PDE for $u$. The finite element method turns this PDE into algebra in four steps:

  1. Discretize — replace the solid with a tetrahedral mesh (nodes and elements), and approximate $u$ by piecewise-linear functions anchored at the nodes.
  2. Weaken — multiply by a test function and integrate by parts, turning the second derivatives into first derivatives.
  3. Assemble — sum each element's contribution into a global sparse matrix $K$ and a load vector $f$.
  4. Solve — impose the supports and solve $K u = f$ for the nodal displacements.

For a linear tetrahedron the strain is constant inside each element, so the element contributions have a closed form — no quadrature error.

The constitutive matrix

Hooke's law for an isotropic material needs two numbers: Young's modulus $E$ (stiffness) and Poisson's ratio $\nu$ (how much it bulges sideways when squeezed). elasticity_matrix packs them into the 6×6 matrix that relates stress to strain:

Rasmah.elasticity_matrix — Function
elasticity_matrix(E, ν) -> FixedSizeMatrix

The 6×6 isotropic linear-elastic constitutive matrix in Voigt form, from Young's modulus E and Poisson's ratio ν.

Theory

Hooke's law for an isotropic solid relates stress $\boldsymbol{\sigma}$ to strain $\boldsymbol{\varepsilon}$ through two independent constants. Rasmah uses the Lamé pair

\[\mu = \frac{E}{2(1+\nu)}, \qquad \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)} ,\]

where $\mu$ is the shear modulus and $\lambda$ is the first Lamé constant. In Voigt notation the stress and strain are written as 6-vectors (normal components first, then the three engineering shears), so the constitutive law $\sigma = C \varepsilon$ becomes a plain 6×6 matrix product:

\[C = \begin{pmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{pmatrix}.\]

The diagonal shear block is $\mu$ (not $2\mu$) because the Voigt strain stores engineering shear strains.

Arguments

  • E::Real — Young's modulus (a pressure).
  • ν::Real — Poisson's ratio (dimensionless, $-1 < \nu < 0.5$).

Returns

A FixedSizeMatrix (6×6). Passing ν = 0.5 (the incompressible limit) throws, because the Lamé formulation is singular there; use a mixed/hybrid formulation instead.

Example

The shear diagonal entry is the Lamé modulus $\mu$:

julia> C = elasticity_matrix(200e9, 0.3);

julia> C[4, 4] == 200e9 / (2 * (1 + 0.3))
true

julia> C[1, 1] - 2 * C[4, 4] == C[1, 2]
true
source
elasticity_matrix(200e9, 0.3)
6×6 FixedSizeArrays.FixedSizeArray{Float64, 2, Memory{Float64}}:
 2.69231e11  1.15385e11  1.15385e11  0.0         0.0         0.0
 1.15385e11  2.69231e11  1.15385e11  0.0         0.0         0.0
 1.15385e11  1.15385e11  2.69231e11  0.0         0.0         0.0
 0.0         0.0         0.0         7.69231e10  0.0         0.0
 0.0         0.0         0.0         0.0         7.69231e10  0.0
 0.0         0.0         0.0         0.0         0.0         7.69231e10

The diagonal shear block is the Lamé modulus $\mu = E/(2(1+\nu)) \approx 76.9$ GPa.

Element and global stiffness

Each element contributes a 12×12 stiffness $k_e = V\,B^{T} C B$. The matrix $B$ (the strain–displacement matrix) holds the gradients of the shape functions; the volume $V$ is the element's signed volume. The weak-form elasticity_form builds the global stiffness $K$ by scattering every $k_e$ into the shared degrees of freedom, and assemble_matrix turns it into a sparse matrix:

Rasmah.elasticity_form — Function
elasticity_form(U, C) -> BilinearForm

The linear-elasticity BilinearForm on the (vector) finite-element space U with Voigt stiffness matrix C.

source
m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2)
U = FESpace(m, ReferenceFE(3, 1); valuetype = VectorValue{3})
K = assemble_matrix(elasticity_form(U, elasticity_matrix(200e9, 0.3)))
size(K)
(81, 81)

The result is a sparse $81 \times 81$ matrix for a $3 \times 3 \times 3$ node grid ($27$ nodes × $3$ displacement components).

Solving for the displacement

With the stiffness assembled, the displacement is the solution of $K u = f$. Because a body is free to move rigidly until it is supported, the system is only invertible on the free degrees of freedom — the ones not fixed by a support. Name the physics and pass a material, and solve does the assembly and the constrained solve; the Elasticity marker behind that call is documented in the next chapter, Solving with physics.

Rasmah.compliance — Function
compliance(m::TetMesh, material, fixed_dofs, forces[, φ]; p = 3) -> number

The elastic compliance $f^{T} u$ of a structure under its loads, computed from a material. The material supplies its stiffness tensor $C$; the load and support are given by fixed_dofs (the constrained dof indices) and forces (the full-length load vector). An optional per-element density φ SIMP-scales the stiffness ($C \rightarrow φ^{p} C$) for topology optimization.

Example

julia> m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2);

julia> steel = Material(youngs_modulus = 200e9, poisson_ratio = 0.3);

julia> compliance(m, steel, collect(1:3), zeros(3 * size(m.nodes, 2))) == 0.0
true
source

Stress

A displacement field is only useful once it is turned into stress — the quantity that tells you whether the part will fail.

Rasmah.von_mises_stress — Function
von_mises_stress(m, u, material) -> vector
von_mises_stress(m, u, C::AbstractMatrix) -> vector

Per-element von Mises (equivalent) stress from a displacement field.

Theory

The von Mises stress is a single scalar that condenses the 6 stress components into an "equivalent" uniaxial stress — the value a uniaxial test would need to reach the same distortional strain energy. In Voigt components it is

\[\sigma_v = \sqrt{\sigma_1^2 + \sigma_2^2 + \sigma_3^2 - \sigma_1\sigma_2 - \sigma_2\sigma_3 - \sigma_3\sigma_1 + 3(\sigma_4^2 + \sigma_5^2 + \sigma_6^2)} .\]

It is the standard yield criterion (a material yields where $\sigma_v$ reaches its yield strength) and the usual scalar field plotted on stress results.

Arguments

  • m — a TetMesh or TriangleMesh (plane stress).
  • u — the full nodal displacement vector (from solve(Elasticity(), …)).
  • material — the material supplying the stiffness tensor, or C — a Voigt constitutive matrix.

Returns

A per-element vector of von Mises stresses (length size(m.elements, 2)).

Example

A mesh loaded uniformly in one direction has the same von Mises stress in every element:

julia> m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2);

julia> steel = Material(youngs_modulus = 200e9, poisson_ratio = 0.3);

julia> sv = von_mises_stress(m, zeros(3 * size(m.nodes, 2)), steel);

julia> length(sv) == size(m.elements, 2) && all(iszero, sv)
true
source

Heat conduction

Heat conduction is the scalar cousin of elasticity: the temperature $T$ replaces the displacement, and the conductivity $k$ replaces the elastic moduli. The same mesh, the same weak-form idea, but one degree of freedom per node instead of three. The diffusion stiffness assembles exactly like the elasticity stiffness:

Rasmah.heat_stiffness — Function
heat_stiffness(m, k[, φ]; p = 3) -> sparse matrix
heat_stiffness(m, K::AbstractMatrix[, φ]; p = 3) -> sparse matrix

Assemble the heat-conduction (diffusion) stiffness matrix of a mesh.

Theory

Steady heat conduction is Poisson's equation $\nabla\cdot(k\nabla T) = 0$. Its weak form over a test function $v$ is $\int k\, \nabla T \cdot \nabla v\, dV$, which discretizes to the sparse matrix $K_T$ with entries $K_{ij} = \int k\, \nabla N_i \cdot \nabla N_j\, dV$ for the finite-element shape functions $N_i$. This is the exact scalar analogue of the elasticity stiffness (one dof per node instead of three).

Arguments

  • m — a TetMesh or TriangleMesh.
  • k — a scalar conductivity, or K — a 3×3 anisotropic conductivity tensor.
  • φ — optional per-element density (SIMP-scaled, for thermal topology optimization).

Keyword options

  • p = 3 — SIMP penalty exponent: each element scales by $\varphi_e^p$.

Returns

A sparse matrix of size n × n (one dof per node).

Example

julia> m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2);

julia> size(heat_stiffness(m, 1.0)) == (27, 27)
true
source

The one-line solve wrapper for this equation is the HeatConduction marker, documented in Solving with physics.

Post-processing

Raw fields from an optimizer or a coarse solve can be noisy. Laplacian smoothing averages each value with its neighbours to produce a clean field for display.

Rasmah.smooth_field — Function
smooth_field(m, values; iterations = 10) -> vector

Laplacian smoothing of a nodal field over the mesh adjacency.

Theory

The discrete Laplacian of a nodal field replaces each value with the average of its neighbours; iterating this diffuses the field across the mesh. It is the standard post-processing pass for checkerboard-free density/stress plots in topology optimization, where raw element fields are noisy from the SIMP penalization.

Arguments

  • m — a TetMesh or TriangleMesh (only its adjacency is used).
  • values — the nodal (per-node) field to smooth.

Keyword options

  • iterations = 10 — number of diffusion passes.

Returns

A copy of values after iterations smoothing passes (isolated nodes keep their value).

Example

Smoothing a constant field leaves it unchanged:

julia> m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2);

julia> smooth_field(m, ones(size(m.nodes, 2))) == ones(size(m.nodes, 2))
true
source

A worked example: a bar in tension

Put the pieces together: mesh a unit cube, clamp its $x = 0$ face, pull the $x = 1$ face with a total force of $10^6$ N, and solve for the displacement.

m = tetrahedralize_box(1.0, 1.0, 1.0, 6, 6, 6)
left  = nodes_in_box(m, [0.0, 0.0, 0.0], [0.0, 1.0, 1.0])
right = nodes_in_box(m, [1.0, 0.0, 0.0], [1.0, 1.0, 1.0])

steel = Material(youngs_modulus = 200e9, poisson_ratio = 0.3)
fixed_dofs = reduce(vcat, [[3i - 2, 3i - 1, 3i] for i in left.ids])
forces = zeros(3, size(m.nodes, 2))
for i in right.ids
    forces[1, i] = 1e6 / length(right.ids)
end

u = solve(Elasticity(), m, steel; fixed = fixed_dofs, loads = vec(forces))
1029-element QuantityVector{Float64}:
 0.0 m
 0.0 m
 0.0 m
 9.759903467333744e-7 m
 4.830757127215885e-7 m
 4.830757127215873e-7 m
 1.8315593804842913e-6 m
 6.429205866873113e-7 m
 6.429205866873084e-7 m
 2.7286592873081523e-6 m
 ⋮
 3.450216801373272e-6 m
 -6.515228466500711e-7 m
 -6.515228466500599e-7 m
 4.5866630847911804e-6 m
 -7.123302576520053e-7 m
 -7.123302576519914e-7 m
 5.96986166244414e-6 m
 -8.916886025407006e-7 m
 -8.916886025406818e-7 m

The maximum displacement is:

round(maximum(abs.(ustrip(u))); digits = 10)
7.5298e-6

The bar stretches $\approx 7.5 \times 10^{-6}$ m. The elementary rod formula predicts $u = F L/(E A) = 5 \times 10^{-6}$ m; the finite-element answer is a little larger because clamping the whole $x = 0$ face constrains lateral Poisson contraction, which stiffens a thin boundary layer. Refining the mesh narrows the gap (Saint-Venant's principle).

The compliance $f^{T} u$ measures the work done by the load:

round(compliance(m, steel, fixed_dofs, vec(forces)); digits = 3)
5.066

The analytic value $F^{2}L/(EA) = 5.0$ J is within $1.5\%$ of the FEM answer, as expected on a coarse mesh.

Finally, colour the bar by its von Mises stress — uniform across the interior, as a rod under pure tension should be:

sv = von_mises_stress(m, ustrip(u), steel)
vtk(m; field = sv, fieldname = "von Mises stress", edges = false)

Next steps

The assemblers here are the engine. To run a simulation in one line — naming the physics and handing over a material and boundary conditions — continue to Solving with physics.