Finite Element Method

The setup

Take an elliptic PDE on a 2D domain Ω with boundary ∂Ω — for example Poisson's equation:

−Δu(x, y) = f(x, y)        in Ω
       u  = 0               on ∂Ω

The strong form requires u to be twice differentiable everywhere, which is too restrictive for general engineering geometries with corners, holes, and material interfaces. The first move is to derive a weak form: multiply both sides by any sufficiently smooth test function v vanishing on the boundary, integrate over Ω, and integrate the left side by parts:

∫_Ω ∇u · ∇v dx = ∫_Ω f · v dx        for all admissible test v

This needs only first derivatives of u and v (and only in L²) instead of second derivatives. The space of admissible functions is H¹₀(Ω) — Sobolev space of square-integrable functions with square-integrable gradients vanishing on ∂Ω.

Meshing the domain

Partition Ω into a finite mesh T_h of simple elements:

• 1D: intervals [x_i, x_{i+1}]
• 2D: triangles or quadrilaterals
• 3D: tetrahedra or hexahedra

The mesh size h is the diameter of the largest element. The mesh must conform — elements meet face-to-face, not in T-junctions — and avoid degenerate (sliver) elements with very small angles, which break the convergence rate. Mesh generators like Gmsh, Triangle, TetGen, and CGAL produce conforming meshes for any reasonable CAD geometry; this preprocessing is often the most labour-intensive step of an FEM workflow.

Hat basis functions

At each mesh node i, define a basis function φ_i that is:

• piecewise-polynomial on each element (the simplest case is linear),
• equal to 1 at node i,
• equal to 0 at every other mesh node,
• continuous across element boundaries.

In 1D, φ_i is a hat: linear on [x_{i−1}, x_i], peaking at x_i, linear back down on [x_i, x_{i+1}], and zero outside [x_{i−1}, x_{i+1}]. In 2D on a triangular mesh, φ_i is a pyramidal hat supported on the cluster of triangles meeting at node i. The key property is compact support: each φ_i is nonzero on only a handful of elements, so when two nodes i and j are not connected in the mesh graph, their supports do not overlap and ∫∇φ_i·∇φ_j dx = 0. The resulting stiffness matrix is therefore sparse — at most ~10 nonzeros per row in 2D triangular meshes, ~30 in 3D tetrahedral meshes.

Galerkin's substitution

Approximate u in the finite-dimensional subspace V_h = span{φ_1, …, φ_N}:

u_h(x) = Σ_{i=1}^N c_i · φ_i(x)

The unknowns are the N coefficients c_i — the values of u at the mesh nodes. Substitute into the weak form, requiring it to hold for every basis function v = φ_j:

Σ_i c_i · ∫_Ω ∇φ_i · ∇φ_j dx = ∫_Ω f · φ_j dx        for j = 1, …, N

That is a linear system K c = f where the stiffness matrix K and load vector f have entries

K_ij = ∫_Ω ∇φ_i · ∇φ_j dx      f_j = ∫_Ω f · φ_j dx

K is N × N, sparse (~10–30 nonzeros per row), symmetric (since ∇φ_i·∇φ_j = ∇φ_j·∇φ_i), and positive-definite (for our elliptic problem). The system has a unique solution which can be computed by sparse direct factorisation (CHOLMOD, MUMPS) or iterative methods (preconditioned conjugate gradient, multigrid). The solve cost is O(N^{3/2}) for direct in 2D, O(N²) in 3D, but O(N) for optimal multigrid — which is why multigrid dominates large-scale FEM in production.

Element-by-element assembly

The integrals defining K_ij and f_j decompose into sums over individual elements:

K_ij = Σ_{e ∈ T_h} ∫_{K_e} ∇φ_i · ∇φ_j dx

For each element K_e, only a few basis functions are nonzero on K_e — the ones associated with the element's own nodes. Compute the small local stiffness matrix (3 × 3 for a linear triangle, 10 × 10 for a quadratic tetrahedron) by Gauss quadrature over K_e, then scatter it into K via the connectivity map that translates between local element-node indices and global mesh-node indices. Loop over all elements, accumulate, and you have K. The loop is embarrassingly parallel; production codes (deal.II, FEniCS) parallelise over elements with OpenMP or MPI.

Worked example: 1D Poisson on [0, 1]

Solve −u'' = 1 on [0, 1] with u(0) = u(1) = 0. The exact solution is u(x) = x(1 − x)/2 with peak u(1/2) = 1/8 = 0.125.

Mesh: 4 equal intervals, h = 0.25, with nodes x_0 = 0, x_1 = 0.25, x_2 = 0.5, x_3 = 0.75, x_4 = 1. After removing the Dirichlet nodes (x_0 and x_4 = 0 are fixed), we have 3 free unknowns c_1, c_2, c_3 — the values of u_h at the interior nodes.

The local stiffness on each interval of length h is (1/h)·[[1, −1], [−1, 1]]. Assembling globally:

K = (1/h) · |  2  -1   0 |       f = h · | 1 |
            | -1   2  -1 |               | 1 |
            |  0  -1   2 |               | 1 |

For h = 0.25 this is K = 4·tridiag(−1, 2, −1) and f = 0.25·[1, 1, 1]ᵀ. Solving:

c_1 = 0.09375   (true value 0.09375)
c_2 = 0.12500   (true value 0.12500)
c_3 = 0.09375   (true value 0.09375)

For this problem, linear FEM with h = 0.25 is exact at the nodes — a special feature of 1D Poisson because the true solution is a quadratic and the basis hat functions interpolate it perfectly at the nodes. In general, the error at the nodes scales as O(h^{p+1}) where p is the polynomial degree of the basis — for linear hats in 2D and 3D, O(h²); for quadratic basis functions, O(h³); for cubic, O(h⁴).

Convergence rates by element type

Element	Degree p	Local DoFs	L² error	H¹ error	Use
Linear (P1) triangle	1	3	O(h²)	O(h)	Default, easy meshing
Quadratic (P2) triangle	2	6	O(h³)	O(h²)	Smooth solutions, default in ABAQUS
Cubic (P3) triangle	3	10	O(h⁴)	O(h³)	Fluid dynamics, high accuracy
Bilinear quad (Q1)	1	4	O(h²)	O(h)	Structured meshes, simpler than P2
Linear tetrahedron	1	4	O(h²)	O(h)	3D structural, easy to mesh
Quadratic tetrahedron	2	10	O(h³)	O(h²)	3D engineering analysis, ANSYS default
hp-FEM	variable	variable	exp(-c√N)	exp(-c√N)	Smooth solutions, spectral-like
DG (Discontinuous Galerkin)	p	(p+1)(p+2)/2 per elem	O(h^{p+1})	O(h^p)	Conservation laws, advection

The pattern: each additional polynomial degree adds one order of convergence in L² and one in H¹ (energy norm). Higher-degree elements have more local degrees of freedom and denser local matrices — a quadratic triangle has six DoFs vs the linear triangle's three — but reach the same accuracy on coarser meshes. For smooth solutions, the trade-off favours higher order; for solutions with corner singularities or shocks, lower order with adaptive refinement is more robust.

Where FEM shows up

• Structural engineering. Every modern building's structural model is an FEM mesh of tens of millions of elements. SAP2000, ETABS, RAM, RFEM, and ABAQUS run the analysis. The Burj Khalifa, the Millau Viaduct, and the One World Trade Center were all dimensioned via FEM solves on full 3D geometries.
• Automotive crash simulation. ABAQUS Explicit, LS-DYNA, and PAM-CRASH simulate frontal-impact crashes on 10–50 million element vehicle meshes in 24–72 hours of cluster time. Every modern car undergoes virtual crash testing before any physical prototype.
• Aerospace stress and flutter analysis. NASTRAN was developed by NASA in the 1960s for the Apollo programme; today it sizes every airliner wing (Boeing 787, Airbus A350), every turbine blade (GE9X, Rolls-Royce Trent), and predicts flutter velocities on every flight envelope.
• Computational fluid dynamics. COMSOL, ANSYS Fluent, OpenFOAM, and SU2 all offer FEM (or finite-volume on FEM-like meshes) for incompressible and compressible flow. Stabilised Galerkin methods like SUPG handle the advection-diffusion-reaction equations in chemical reactors and combustion chambers.
• Electromagnetics. Maxwell's equations on intricate 3D geometries — antenna design, transformer modelling, MRI gradient coil simulation — use edge-element FEM (Nédélec elements) to preserve the curl structure. Codes: ANSYS HFSS, COMSOL Multiphysics, FreeFEM++.
• Biomedical engineering. Pre-surgical planning for orthopaedic implants, dental restorations, cardiac valve replacements, and craniomaxillofacial reconstructions all use FEM stress analysis on patient-specific CT/MRI-derived meshes. ANSYS Discovery, FEBio, and SimVascular are the leading codes.
• Geomechanics and reservoir simulation. Oil and gas reservoir modelling (Schlumberger Eclipse, CMG GEM, Roxar Tempest), tunnel and dam stability analysis, earthquake-induced soil liquefaction — all run on FEM solvers for coupled solid-fluid PDEs.

Variants and extensions

• hp-FEM. Combines h-refinement (smaller elements) and p-refinement (higher polynomial degree) adaptively to chase singularities efficiently. Achieves exponential convergence on smooth problems; standard in commercial codes (ANSYS, COMSOL).
• Discontinuous Galerkin (DG). Basis functions are not enforced to be continuous between elements; instead, fluxes across element faces are penalised. Combines FEM's geometric flexibility with the conservation properties of finite-volume methods. Dominant in shock-capturing aerodynamics codes.
• Spectral element method. Uses very high polynomial degree (8 or 16) on a coarse mesh with Gauss-Lobatto-Legendre node positions. Achieves spectral accuracy on smooth solutions. Used in turbulence DNS (Nek5000, Nektar++) and meteorology (HOMME).
• Mixed FEM. Treats different solution variables (e.g. velocity and pressure in Stokes flow) with different basis spaces satisfying the inf-sup (LBB) condition. Standard in incompressible CFD (Taylor-Hood elements).
• Extended FEM (XFEM, GFEM). Enriches the basis with non-polynomial functions to represent cracks, discontinuities, and singularities without remeshing. Used in fracture mechanics and composite materials.
• Isogeometric analysis (IGA). Replaces piecewise-polynomial basis with the NURBS basis functions used in CAD — eliminating the geometry-to-mesh approximation step. Hughes, Cottrell, Bazilevs (2005); growing adoption in industrial codes.

Common pitfalls

• Locking. Linear elements applied to nearly incompressible elasticity or thin shells become artificially stiff (volumetric locking, shear locking). Use mixed formulations or higher-order elements; or reduced/selective integration with hourglass stabilisation.
• Distorted elements. Mesh elements with extreme aspect ratios or small interior angles destroy the constants in the error estimate, sometimes pre-asymptotically. Mesh quality metrics (skewness, Jacobian determinant) flag bad elements before solving.
• Boundary condition treatment. Mixed (Robin) and Neumann conditions show up in the load vector via boundary integrals; ignoring them is a common bug. Dirichlet conditions are enforced by removing the corresponding rows and columns of K or by Lagrange multipliers.
• Singularities at re-entrant corners. Solutions of Poisson's equation on L-shaped domains have unbounded gradient at the inner corner; uniform refinement converges very slowly. Adaptive refinement (a-posteriori error estimators driving local refinement) is essential.
• Inverting the wrong system. The mass matrix M_ij = ∫φ_iφ_j dx and the stiffness matrix K_ij = ∫∇φ_i·∇φ_j dx are routinely confused. M appears in time-dependent and eigenvalue problems; K appears in static equilibrium. Mixing them up produces nonsensical results.

Why FEM dominates engineering

FEM was systematised by Argyris (1954), Clough (who coined the name in 1960), Zienkiewicz (whose 1967 textbook is still in print), and Gilbert Strang in the 1970s. By 1985, NASTRAN, ABAQUS, ANSYS, and ADINA were the dominant commercial FEM codes; by 2000, they were running on every engineer's workstation. Today every major piece of physical engineering on the planet — buildings, bridges, cars, aeroplanes, ships, pipelines, dams, turbines, implants — passes through an FEM solver before being built. The reason is geometric flexibility: triangular and tetrahedral meshes fit any shape, the Galerkin framework handles any second-order PDE, and the resulting sparse linear systems are solvable to engineering accuracy on commodity hardware. Finite differences are easier for slab geometries; finite volumes are sharper for conservation laws; spectral methods win on smooth periodic problems — but FEM is the universal default whenever geometry matters more than spectral accuracy, which is most of engineering simulation.