Finite elements

Deep dive

How it works

Parts:

• mesh — nodes and element connectivity (triangles, tetrahedra)
• basis — hat functions φ_a, linear on each element
• element routine — 3×3 stiffness K_e and load F_e per triangle
• assembler — scatter-adds K_e into the global sparse matrix (CSR)
• boundary step — eliminates Dirichlet nodes, adds Neumann terms
• solver — CG with multigrid, or sparse Cholesky
• post-processor — gradients, fluxes, error indicators for refinement

Operations: weak form, mesh and basis, element matrices, assembly, boundary conditions, solve (steps below).

1. Weak form. Multiply −Δu = f by a test function v vanishing on the Dirichlet boundary, integrate by parts: ∫∇u·∇v dx = ∫f v dx (+ boundary flux terms for Neumann conditions).
2. Mesh and basis. Triangulate the domain. φ_a = hat function equal to 1 at node a, 0 at all other nodes, linear on each triangle. u_h = Σ u_a φ_a.
3. Element matrices. For each triangle, ∇φ_a is constant: ∇φ_a = (y_b − y_c, x_c − x_b)/det for (a, b, c) cyclic, det = twice the signed area. K_e[a][b] = area·∇φ_a·∇φ_b; load F_e[a] ≈ f(centroid)·area/3.
4. Assembly. Scatter-add K_e into the global matrix at (node_a, node_b) for all element pairs; collect triplets, sort, merge duplicates into CSR.
5. Boundary conditions. Dirichlet: eliminate those nodes (only interior nodes are unknowns; move K_ab·g_b to the right-hand side for nonzero g). Neumann: add ∫_Γ g·φ_a.
6. Solve K·u = F; then post-process (gradients, fluxes, error indicators for adaptive refinement).

  node a has hat φ_a       one triangle (a, b, c):
      /\                     K_e = area · [∇φ_i · ∇φ_j]  (3×3)
     /  \                    scatter into rows/cols a, b, c
  __/    \__

Code

pub struct Mesh { pub p: Vec<[f64; 2]>, pub tri: Vec<[usize; 3]>, pub boundary: Vec<bool> }

/// Unit square split into 2 m^2 right triangles, (m+1)^2 nodes.
pub fn unit_square(m: usize) -> Mesh {
    let (mut p, mut boundary, mut tri) = (Vec::new(), Vec::new(), Vec::new());
    for j in 0..=m {
        for i in 0..=m {
            p.push([i as f64 / m as f64, j as f64 / m as f64]);
            boundary.push(i == 0 || j == 0 || i == m || j == m);
        }
    }
    let id = |i: usize, j: usize| j * (m + 1) + i;
    for j in 0..m {
        for i in 0..m {
            tri.push([id(i, j), id(i + 1, j), id(i + 1, j + 1)]); // counter-clockwise
            tri.push([id(i, j), id(i + 1, j + 1), id(i, j + 1)]);
        }
    }
    Mesh { p, tri, boundary }
}

/// P1 element: K_ab = area * grad(phi_a) . grad(phi_b). Returns (K_e, area).
pub fn p1_stiffness(x: [[f64; 2]; 3]) -> ([[f64; 3]; 3], f64) {
    let det = (x[1][0] - x[0][0]) * (x[2][1] - x[0][1]) - (x[2][0] - x[0][0]) * (x[1][1] - x[0][1]);
    let area = 0.5 * det.abs();
    let g: [[f64; 2]; 3] = std::array::from_fn(|a| {
        let (b, c) = ((a + 1) % 3, (a + 2) % 3);
        [(x[b][1] - x[c][1]) / det, (x[c][0] - x[b][0]) / det]
    });
    (std::array::from_fn(|a| std::array::from_fn(|b| area * (g[a][0] * g[b][0] + g[a][1] * g[b][1]))), area)
}

pub struct Csr { pub n: usize, pub row_ptr: Vec<usize>, pub col: Vec<u32>, pub val: Vec<f64> }

impl Csr {
    fn from_triplets(n: usize, mut t: Vec<(u32, u32, f64)>) -> Csr {
        t.sort_unstable_by_key(|&(r, c, _)| (r, c));
        let (mut row_ptr, mut col, mut val) = (vec![0usize; n + 1], Vec::new(), Vec::<f64>::new());
        for (k, &(r, c, v)) in t.iter().enumerate() {
            if k > 0 && t[k - 1].0 == r && t[k - 1].1 == c { *val.last_mut().unwrap() += v; continue; } // merge duplicates
            col.push(c);
            val.push(v);
            row_ptr[r as usize + 1] += 1;
        }
        for i in 0..n { row_ptr[i + 1] += row_ptr[i]; }
        Csr { n, row_ptr, col, val }
    }
    fn mul(&self, x: &[f64], y: &mut [f64]) {
        for i in 0..self.n {
            y[i] = (self.row_ptr[i]..self.row_ptr[i + 1]).map(|k| self.val[k] * x[self.col[k] as usize]).sum();
        }
    }
}

/// Assemble -Laplace(u) = f with u = 0 on the boundary; unknowns are interior nodes only.
pub fn assemble(m: &Mesh, f: impl Fn(f64, f64) -> f64) -> (Csr, Vec<f64>, Vec<Option<u32>>) {
    let mut dof = vec![None; m.p.len()];
    let mut n = 0u32;
    for (a, &b) in m.boundary.iter().enumerate() { if !b { dof[a] = Some(n); n += 1; } }
    let (mut trip, mut rhs) = (Vec::with_capacity(9 * m.tri.len()), vec![0.0; n as usize]);
    for t in &m.tri {
        let x = [m.p[t[0]], m.p[t[1]], m.p[t[2]]];
        let (ke, area) = p1_stiffness(x);
        let fc = f((x[0][0] + x[1][0] + x[2][0]) / 3.0, (x[0][1] + x[1][1] + x[2][1]) / 3.0);
        for a in 0..3 {
            let Some(ra) = dof[t[a]] else { continue };
            rhs[ra as usize] += fc * area / 3.0; // one-point (centroid) quadrature
            for b in 0..3 {
                if let Some(cb) = dof[t[b]] { trip.push((ra, cb, ke[a][b])); }
            }
        }
    }
    (Csr::from_triplets(n as usize, trip), rhs, dof)
}

/// Plain conjugate gradient for SPD systems.
pub fn cg(a: &Csr, b: &[f64], tol: f64) -> (Vec<f64>, usize) {
    let n = b.len();
    let (mut x, mut r, mut p, mut ap) = (vec![0.0; n], b.to_vec(), b.to_vec(), vec![0.0; n]);
    let dot = |u: &[f64], v: &[f64]| u.iter().zip(v).map(|(a, b)| a * b).sum::<f64>();
    let (mut rr, bb) = (dot(&r, &r), dot(b, b));
    for it in 0..10 * n {
        if rr <= tol * tol * bb { return (x, it); }
        a.mul(&p, &mut ap);
        let alpha = rr / dot(&p, &ap);
        for i in 0..n { x[i] += alpha * p[i]; r[i] -= alpha * ap[i]; }
        let rr_new = dot(&r, &r);
        for i in 0..n { p[i] = r[i] + rr_new / rr * p[i]; }
        rr = rr_new;
    }
    (x, 10 * n)
}

Complexity and constants

• Assembly: O(elements), 9 triplets per P1 triangle (16 per tetrahedron); sorting triplets is O(nnz log nnz), fine to order of 10⁷ entries. Production codes precompute the sparsity pattern once and scatter directly.
• Matrix: on this structured mesh each interior row stores 7 entries, and the two diagonal-neighbour couplings cancel to explicit zeros, leaving exactly the 5-point stencil [−1, −1, 4, −1, −1]; unstructured 3D P1 meshes average order of 15 nonzeros per row.
• Measured with the code above on −Δu = 2π² sin(πx) sin(πy): max nodal error 5.4·10⁻³ at m = 16, 1.3·10⁻³ at m = 32, 3.3·10⁻⁴ at m = 64 (ratio 4, second order); CG to 10⁻¹² relative residual takes 26, 57, 114 iterations (doubling per halving of h).
• Condition number of K grows as O(h⁻²): CG iterations O(h⁻¹); multigrid-preconditioned CG keeps them near constant.

Efficiency

• Cost: assembly is O(elements), 9 triplets per P1 triangle (16 per tetrahedron); sorting triplets is O(nnz log nnz). Structured 2D P1 rows hold 5 non-zeros, unstructured 3D P1 meshes order of 15; a CSR matrix streams order of 12 bytes per non-zero per SpMV. Measured with the code above: CG to 10⁻¹² takes 26, 57, 114 iterations at m = 16, 32, 64. Bound: memory bandwidth in SpMV; assembly scatter for small problems.
• Levers:
  • precomputed sparsity pattern, direct scatter — skips the triplet sort (page above) — repeated assembly in time stepping and Newton loops
  • element coloring or per-thread assembly — scatter without atomics — multicore assembly
  • node reordering (reverse Cuthill-McKee, nested dissection) — less fill-in, better SpMV locality — direct solvers and large meshes
  • multigrid-preconditioned CG — iterations near constant instead of O(h⁻¹) (page above) — fine meshes
  • matrix-free sum factorization — less memory traffic than a stored matrix (page above) — polynomial order ≥ 3
  • lumped mass matrix — explicit steps need no solve — wave and dynamics problems
• Naive → standard → state of the art:
  • dense K and dense LU · O(N²) memory, O(N³) time
  • sparse assembled K, CG or sparse Cholesky · 114 iterations at m = 64 (measured above)
  • matrix-free high-order operators, multigrid-preconditioned CG, MPI or GPU · iterations near constant in mesh size

Implementation notes

• Element loop is embarrassingly parallel; the scatter is not. Colour elements so no two in a colour share a node, or assemble per-thread and reduce, or use atomics.
• Quadrature: P1 Poisson needs only the constant gradients; the load and variable coefficients need Gauss rules on the reference triangle (1-point centroid is exact for linear integrands).
• Reference element mapping: compute the Jacobian of the affine (or isoparametric) map per element; store det J and J⁻¹ instead of recomputing for every quadrature point.
• Numbering: reorder nodes (reverse Cuthill-McKee or nested dissection) to reduce fill-in for direct solvers and improve cache locality of SpMV.
• Matrix-free high-order: evaluate K·v on the fly with sum factorization instead of storing K; at order ≥ 3 this beats assembled matrices in memory traffic.
• Mass matrix for time-dependent problems: P1 triangle M_e = area/12·[[2,1,1],[1,2,1],[1,1,2]]; lumping it (row sums on the diagonal) makes explicit time stepping cheap.
• Mesh generation (Gmsh, TetGen, Triangle) and mesh quality checks are half the work in practice.

Pitfalls

• Clockwise elements produce a negative det: using det instead of |det| for the area flips the sign of K_e.
• Imposing Dirichlet conditions by overwriting rows without fixing columns: the matrix loses symmetry and CG fails; eliminate, or zero rows and columns and move terms to the right-hand side.
• Linear triangles/tetrahedra for nearly incompressible elasticity or thin bending: locking (far too stiff); use quadratic, mixed or reduced-integration elements.
• Sliver elements (near-zero volume): huge condition numbers and local garbage.
• Measuring error only at nodes: some meshes are superconvergent there; check L² and energy norms.

Variants and alternatives

In production

Engineering: FEniCS/DOLFINx, Firedrake, deal.II, MFEM, libMesh, FreeFEM, NGSolve; commercial Abaqus, ANSYS Mechanical, COMSOL, LS-DYNA (explicit dynamics), Nastran. Geophysics: ASPECT (mantle convection) is built on deal.II. Graphics and animation: the Incremental Potential Contact method simulates deformable solids with contact on FEM meshes.

Implementations

Sources

• Gilbert Strang, George J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, 1973.
• Susanne C. Brenner, L. Ridgway Scott, The Mathematical Theory of Finite Element Methods, 3rd ed., Springer, 2008. book
• Thomas J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover, 2000.
• Mats G. Larson, Fredrik Bengzon, The Finite Element Method: Theory, Implementation, and Applications, Springer, 2013. book
• Anders Logg, Kent-Andre Mardal, Garth N. Wells (eds.), Automated Solution of Differential Equations by the Finite Element Method, Springer, 2012. book
• Minchen Li et al., Incremental Potential Contact: Intersection- and Inversion-free, Large-Deformation Dynamics, ACM TOG 39(4), 2020. paper · code
• Martin Kronbichler, Timo Heister, Wolfgang Bangerth, High accuracy mantle convection simulation through modern numerical methods, Geophysical Journal International 191(1), 2012. paper · code