Finite Elements

The finite element method solves partial differential equations by chopping the domain into small pieces ("elements"), pretending the solution is a simple function on each piece, and demanding that this piecewise solution satisfy the equation in an averaged sense. The averaging trick — multiplying by a "test function" and integrating — turns a differential equation into a linear system you can solve on a computer.

This series builds the method from scratch. We write our own mesher, our own assembly, and our own visualizers. Solutions get plotted alongside their analytic counterparts so you can see convergence with your own eyes as the mesh refines. By the end we have a working solver for Poisson's equation on a 2D square, in well under 200 lines of code, with the assembly routine structured exactly like a textbook FEM library — just smaller.

Chapters

  1. The weak form — Why we don't solve the PDE directly. Multiplying by a test function and integrating turns "second derivatives" into "first derivatives of pairs," which is exactly the trade we want to make.
  2. A 1D mesh from scratch — Nodes, elements, and the piecewise-linear basis functions ("hat functions") that live on them. With a visualizer that lets you refine the mesh and see how the basis adapts.
  3. Solving 1D Poisson — Element stiffness matrices, global assembly, boundary conditions, solve, plot. Compare to the analytic solution and watch the error shrink as the mesh refines.
  4. A 2D triangular mesh — Structured triangulation of the unit square. Nodes, elements, boundary nodes, and a visualizer that shows the whole thing.
  5. Solving 2D Poisson — Same recipe in 2D: per-triangle stiffness, assembly into the global system, solve, render the solution as a filled-triangle heatmap with a color bar.
  6. Convergence — How fast does FEM actually get there? The O(h²) theorem stated, then measured on the visualizers we built. Log-log plots showing the empirical slope lands at 2.
  7. The heat equation — Adds time. The mass matrix, three time-stepping schemes (forward Euler, backward Euler, Crank-Nicolson), and a live animation that lets you push past the CFL bound and watch forward Euler explode.
  8. Linear elasticity — Vector PDE. Two displacement components per node, the B matrix, plane-stress Hooke's law, and a deformed cantilever beam with von Mises stress as a heatmap.
  9. Delaunay triangulation — Real meshes for non-rectangular domains. The empty-circle property, the flip operation, Bowyer-Watson incremental insertion, and a click-to-add visualizer where the empty-circle property is verifiable by eye when circumcircles are toggled on.
  10. Poisson on a Delaunay mesh — End-to-end. Click points to define a domain; the program Delaunay-triangulates them, identifies the boundary as the convex hull, runs the chapter-5 assembly loop on the resulting unstructured mesh, and renders the solution as a heatmap. Same code, different geometry.

What we're solving

Poisson's equation in one and two dimensions, with zero boundary values:

−u′′(x)=f(x),x∈(0,1),u(0)=u(1)=0-u''(x) = f(x), \quad x \in (0, 1), \qquad u(0) = u(1) = 0
−∇2u(x,y)=f(x,y),(x,y)∈(0,1)2,u=0 on the boundary-\nabla^2 u(x, y) = f(x, y), \quad (x, y) \in (0, 1)^2, \qquad u = 0 \text{ on the boundary}

Both pick up the entire FEM machinery (weak form, basis functions, element stiffness, assembly, boundary handling) but keep the analytic solutions cheap to write down — so we can overlay "what FEM thinks" on "what's actually true" and watch them converge. Nothing here is exotic; the same code structure extends to elasticity, heat transfer, and convection–diffusion by swapping the integrand.

Prerequisites

Calculus through partial derivatives, comfort with linear algebra (matrices, solving Ax = b), and Python with numpy for the implementation. Helpful but not required: any prior exposure to PDEs or to a finite-difference method, since we'll spend time contrasting FEM with FD where it matters.