Poisson on a Delaunay mesh: end-to-end

Nine chapters of FEM machinery. One chapter for the mesher. Now we run them together.

The point

The chapter-5 assembly loop never knew it was running on a structured grid. Every line of it took an array of nodes and an array of elements and did its work without asking how those came to exist. Swap in a Delaunay triangulation and the same code solves Poisson on whatever shape your points happen to enclose.

That genericity is the whole reason FEM exists. Finite-difference methods are wedded to grids; FEM works on any triangulation. We wrote a structured mesher first because it was easy. We wrote a Delaunay mesher because it's what real problems need. They both plug into the same downstream pipeline.

Identifying the boundary

The chapter-4 mesher gave us boundary nodes for free — they were the nodes on the edges of the unit square, identifiable by index. A Delaunay triangulation has no such structure. We have to extract the boundary from the triangulation itself.

The trick: a triangulation edge sits on the boundary if and only if it appears in exactly one triangle. Internal edges are shared by two triangles (the two on either side); boundary edges have no triangle on the outside. Walk all the edges with a counter:

def find_boundary_nodes(triangles: list[tuple[int, int, int]]) -> set[int]:
    """The boundary of a triangulation is the set of edges that appear in
    exactly one triangle. The nodes on those edges are boundary nodes.

    For a Delaunay triangulation of a point set with no constraints, the
    boundary equals the convex hull of the points."""
    from collections import Counter
    edge_count = Counter()
    for a, b, c in triangles:
        for u, v in [(a, b), (b, c), (c, a)]:
            edge_count[(min(u, v), max(u, v))] += 1

    boundary = set()
    for a, b, c in triangles:
        for u, v in [(a, b), (b, c), (c, a)]:
            if edge_count[(min(u, v), max(u, v))] == 1:
                boundary.add(u)
                boundary.add(v)
    return boundary

For an unconstrained Delaunay triangulation of a point set, the boundary nodes form exactly the convex hull. So this routine is also a roundabout way to compute the convex hull — a side benefit. If the point cloud has any concavity in it (an L-shape, say), Delaunay alone produces the convex hull anyway, filling in what should have been the concave part with extra triangles. Real production code uses constrained Delaunay to honor a specified domain boundary; that's its own algorithm and the natural follow-up to this chapter.

The pipeline

Wire it together: triangulate, find the boundary, assemble, apply boundary conditions, solve. The middle three steps are untouched chapter-5 code:

def poisson_on_delaunay(points: list[tuple[float, float]],
                       f: callable) -> tuple[np.ndarray, list, set]:
    """End-to-end pipeline: triangulate the point cloud, identify the
    boundary, assemble the Poisson system, apply Dirichlet u=0, solve."""
    triangles = bowyer_watson(points)              # chapter 9
    boundary = find_boundary_nodes(triangles)

    # The chapter-5 assembly loop, unchanged. The only thing that's different
    # is where 'nodes' and 'elements' came from. Same per-element 3x3
    # stiffness, same 3-point quadrature, same row-replacement BCs.
    nodes = np.asarray(points)
    elements = np.asarray(triangles)
    K, b = assemble_2d(nodes, elements, f)
    apply_dirichlet_2d(K, b, list(boundary))

    u = np.linalg.solve(K, b)
    return u, triangles, boundary

Less than ten lines of new logic on top of everything we'd already written. The architectural payoff: the FEM library and the mesher have nothing to say to each other besides "here are nodes and elements." Either side can be swapped without touching the other.

See it

Click anywhere to add a point. Each click triggers the whole pipeline — re-triangulate, find the new boundary, re-assemble, re-solve, re-color. The forcing-function buttons switch the right-hand side of the PDE. Boundary nodes are drawn hollow, interior nodes are drawn filled.

Things to try:

• Hit 4 corners, then watch the solution as you add points one at a time. Two triangles is barely enough to represent any non-trivial solution; the heatmap is dominated by whatever single value lands at the few interior nodes. Keep adding — the heatmap fills in as the resolution grows.
• Use + 20 (disk) a few times to fill the interior. The heatmap should look smooth, with the solution peaking near the center for the exp(-25 r²) forcing — exactly where you'd expect.
• Switch between the three forcings on the same point set. Each one solves a different PDE on the same mesh; the assembly loop runs from scratch each time. f = 1 gives a smooth dome; exp gives a localized peak; sin·sin gives a tidy bell shape — but only if your point cloud actually fills the unit square so the boundary conditions are imposed where the analytic solution wants them to be.
• Cluster a bunch of points in one corner. The triangulation adapts — denser triangles where points are concentrated, sparser where they're not. That's a primitive form of mesh grading. Production code lets you specify a sizing field directly; here you do it by hand-clicking the desired density.

What's still missing

Compare what's running on this page to a real FEM library (FEniCS, deal.II, libMesh):

• Sparse storage. We allocated a dense $N\times N$ matrix. That breaks for $N>10^4$ or so. Production code uses CSR or similar; the assembly loop is unchanged, only the stiffness-matrix data structure differs.
• Iterative solvers. Our Gauss elimination is $O(N^3)$; conjugate gradient on a Poisson SPD matrix is closer to $O(N^{1.5})$ with a good preconditioner. That's the difference between solving 100 nodes and 100,000.
• Higher-order elements. P2 (quadratic) basis functions converge faster per DOF. Same recipe, bigger element matrices, more nodes per element (6 per triangle for P2 vs. 3 for P1).
• Constrained Delaunay. Our mesher gives the convex hull. Real domains have specified boundaries — potentially curved, potentially with holes. The fix is one algorithmic step (Lawson flips that respect constrained edges) on top of what we have.
• Quality refinement. Ruppert's algorithm iteratively splits poor-quality triangles until every triangle meets a target minimum angle. Couples to Bowyer-Watson at essentially zero cost.
• 3D. Tetrahedra instead of triangles, 4-node element matrices instead of 3-node. The empty-circle property generalizes to "empty-sphere," but slivers are much harder to avoid.
• Parallelism. Modern FEM problems are solved on hundreds of cores. Domain decomposition + parallel sparse solvers are well-developed but a chapter unto themselves.

None of those changes the conceptual recipe. Every production FEM solver, at its core, is the loop you've now seen end-to-end: define a mesh, identify the boundary, assemble per-element contributions, apply boundary conditions, solve. Everything else is engineering for scale, robustness, and accuracy — important work, but downstream of the ideas this series is about.

Where this leaves the series

The series has gone from "here's the weak form" to "here's a working mesher driving a working solver, end-to-end, on arbitrary point clouds." That's a complete arc — every concept has been introduced, made code, made visible, and then exercised. Whatever you build next on top of this — a particular physics problem, a particular geometry, a faster solver — will be a refinement of what's here, not a different game.