Linear elasticity: vector PDEs

Finite Elements

Poisson and the heat equation each have one unknown per point. Elasticity has two — the displacement vector $u = (u_{x}, u_{y})$ . Every node has two degrees of freedom, every element matrix doubles in each dimension (2×2 → 6×6 for triangles), and a new constitutive law connects strain to stress. The assembly recipe is identical; the inside of each element does more work.

The PDE

For static elasticity with no body forces, equilibrium says the divergence of the stress tensor is zero:

- \nabla \cdot σ = f, σ = D ε, ε = \frac{1}{2} (\nabla u + \nabla u^{⊤})

Three equations linked together: equilibrium relates stress and body force; Hooke's law relates stress and strain; the strain-displacement relation defines strain as the symmetric part of the displacement gradient. Bake them into a single weak form and the FEM machinery from chapter 5 carries over.

Voigt notation

In 2D, strain has three independent components: $ε_{x} = \partial u_{x} / \partial x$ , $ε_{y} = \partial u_{y} / \partial y$ , and the shear $γ_{x y} = \partial u_{x} / \partial y + \partial u_{y} / \partial x$ . Stack them as a 3-vector. Stress stacks the same way: $(σ_{x}, σ_{y}, τ_{x y})$ . This is Voigt notation — the trick that turns the symmetric 2×2 stress and strain tensors into 3-vectors so we can use ordinary matrix-vector products instead of double-contractions.

Hooke's law in plane stress (the assumption: stress in the out-of-plane direction is zero — appropriate for thin plates):

D = \frac{E}{1 - ν ^{2}} 1 ν 0 ν 10 00 (1 - ν) /2

$E$ is Young's modulus (stiffness), $ν$ is Poisson's ratio (lateral contraction under axial load). Plane strain — the assumption for thick bodies — has a slightly different $D$ ; the rest of the machinery is unchanged.

import numpy as np

def D_plane_stress(E: float, nu: float) -> np.ndarray:
    """Plane-stress elasticity matrix (3x3) — relates strain (eps_x, eps_y, gamma_xy)
    to stress (sigma_x, sigma_y, tau_xy) via sigma = D eps."""
    c = E / (1 - nu * nu)
    return c * np.array([
        [1.0, nu,  0.0],
        [nu,  1.0, 0.0],
        [0.0, 0.0, (1 - nu) / 2],
    ])

The B matrix

Strain is a linear function of nodal displacements. Encode that relationship as a 3×6 matrix $B$ per triangle: each row is one strain component, each column is one nodal DOF.

ε = B u_{e}, B = b_{1} 0 c_{1} 0 c_{1} b_{1} b_{2} 0 c_{2} 0 c_{2} b_{2} b_{3} 0 c_{3} 0 c_{3} b_{3}

where $u_{e} = [u_{1 x}, u_{1 y}, u_{2 x}, u_{2 y}, u_{3 x}, u_{3 y}]$ and $(b_{i}, c_{i})$ is the gradient of basis function $ϕ_{i}$ from chapter 5. The first row reads " $ε_{x}$ equals each $b_{i}$ times the corresponding $u_{i x}$ " — exactly the chain rule on $u_{x} = \sum c_{i} ϕ_{i}$ . The third row's $(c_{i}, b_{i})$ pairs encode the cross-derivative shear.

def strain_displacement(verts: np.ndarray) -> tuple[float, np.ndarray]:
    """B matrix (3x6) for one triangle: relates the 6 nodal displacements
    [u1x, u1y, u2x, u2y, u3x, u3y] to the 3 strains [eps_x, eps_y, gamma_xy].
    Returns (area, B). For P1 elements B is constant within the triangle."""
    (x1, y1), (x2, y2), (x3, y3) = verts
    two_A = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)
    A = abs(two_A) * 0.5
    inv = 1.0 / abs(two_A)

    # Each phi_i has gradient (b_i, c_i):
    b1, c1 = (y2 - y3) * inv, (x3 - x2) * inv
    b2, c2 = (y3 - y1) * inv, (x1 - x3) * inv
    b3, c3 = (y1 - y2) * inv, (x2 - x1) * inv

    B = np.array([
        [b1, 0,  b2, 0,  b3, 0 ],
        [0,  c1, 0,  c2, 0,  c3],
        [c1, b1, c2, b2, c3, b3],
    ])
    return A, B

The element stiffness matrix

With $B$ and $D$ in hand, the element stiffness is a one-line matrix expression:

K^{e} = A B^{⊤} D B

A 6×6 symmetric matrix. Compare to chapter 5's scalar Poisson case: $K_{ij}^{e} = A (\nabla ϕ_{i} \cdot \nabla ϕ_{j})$ . Same shape, same $A \cdot (stuff)$ structure; the "stuff" is just bigger because we now have two unknowns per node and a tensor constitutive law.

def element_stiffness_elastic(verts: np.ndarray, D: np.ndarray) -> np.ndarray:
    """K_e = A * B^T D B  (6x6 element stiffness for plane-stress P1 triangle)."""
    A, B = strain_displacement(verts)
    return A * B.T @ D @ B

Assembly

Identical to Poisson, only twice as wide. The DOF ordering convention is "interleaved": node $i$ owns DOFs $2 i$ (its $x$ displacement) and $2 i + 1$ (its $y$ displacement). Element DOFs: $[2 a, 2 a + 1, 2 b, 2 b + 1, 2 c, 2 c + 1]$ .

def assemble_elastic(nodes, elements, D, body_force=None):
    """Assemble K (2N x 2N) and f (2N) for a P1 plane-stress mesh."""
    n = len(nodes)
    K = np.zeros((2 * n, 2 * n))
    f = np.zeros(2 * n)

    for tri in elements:
        verts = nodes[tri]
        K_e = element_stiffness_elastic(verts, D)

        # DOF ordering: [u_a_x, u_a_y, u_b_x, u_b_y, u_c_x, u_c_y].
        dofs = [2 * tri[0], 2 * tri[0] + 1,
                2 * tri[1], 2 * tri[1] + 1,
                2 * tri[2], 2 * tri[2] + 1]
        for r in range(6):
            for s in range(6):
                K[dofs[r], dofs[s]] += K_e[r, s]
        # body force / load handled separately

    return K, f

Boundary conditions split the same way as before. Dirichlet (specified displacement) is essential — bake into the function space by zeroing out the corresponding DOF rows/columns. Neumann (specified traction on a boundary edge) is natural — distribute the traction integral across the edge nodes as a force vector contribution.

See it

A cantilever beam — fixed along its left edge, loaded with a distributed downward traction on its right edge. Real displacements are tiny (the load is small relative to the stiffness), so the visualizer multiplies them by a magnification factor so the bending is visible.

mesh

magnify

fillundeformed outlinemesh overlay

175 nodes · 288 triangles · 350 DOFsactual tip droop: 1.95e-1magnified droop: 9.750

Things to read off the plot:

The bend itself. At magnification 50× or 100× the beam droops noticeably under the tip load. The deflection is largest at the free end and zero at the fixed end — exactly the classical Euler-Bernoulli cantilever profile, recovered without ever assuming beam theory.
Stress concentration at the fixed end. Switch the fill to "von Mises stress" and the colormap shows hottest values along the top and bottom fibers near the wall — that's the bending stress, tension on top, compression on bottom. Refine the mesh and the band gets sharper.
Mesh dependence. The coarse mesh (16×4) gives a noticeably stiffer beam than the fine mesh (32×8) — P1 triangles famously overestimate stiffness in bending ("shear locking"). Refining the mesh fixes it; in production you'd use higher-order elements (P2) or a different formulation to avoid the locking.
The undeformed outline (dashed) is drawn at the original mesh positions. The deformed mesh is the same mesh translated by $mag \cdot u$ per node. The two sit on top of each other at magnification 0.

What was conceptually new

Two new things this chapter introduced, neither hard once you've seen them:

Vector unknowns and the $B$ matrix. Two displacement components per node, packed into a 3×6 strain-displacement matrix. The B matrix is to elasticity what the gradient was to Poisson — the linear operator that turns "values at nodes" into "field quantities of physical interest" (strain, then stress).
The constitutive law as a matrix. Hooke's law in Voigt form is just $σ = D ε$ — same shape as Ohm's law ( $V = R I$ ), just with vectors and a 3×3 matrix. Anisotropic materials, plane-strain assumptions, and 3D extensions all just swap $D$ out for the appropriate version.

Everything else is chapter-5 machinery: assemble per-element contributions into a global system, apply BCs, solve, post-process. The same code structure handles Poisson (1 DOF per node), heat (1 DOF per node, plus a mass matrix), and elasticity (2 or 3 DOFs per node, plus a constitutive law). That genericity is the whole point of building FEM as a recipe rather than a single solver.

Where this is going

Natural next stops, none of which we'll write yet but worth knowing exist: dynamic elasticity (add a mass matrix, get the wave equation, time-step with Newmark); plasticity (the constitutive law becomes nonlinear — incremental assembly, Newton's method per step); 3D (4-node tetrahedra or 8-node hexahedra, same recipe with bigger element matrices); and contact (the boundary condition itself becomes an unknown — variational inequalities, augmented Lagrangians, the rabbit hole of the field).