Delaunay triangulation: real meshes

Finite Elements

Our 2D mesher only handles rectangles. Real FEM problems live on L-shapes, circles with holes, turbine blades, and crash-test dummies. Time to upgrade.

Why structured isn't enough

Three things the structured mesher from chapter 4 can't do:

Curved or non-rectangular domains. A circle, an airfoil, an interior cavity — none of those tile naturally with a Cartesian grid. Forcing one to fit gives staircase boundaries that contaminate the solution.
Local refinement. If the solution has a steep gradient near a corner, you want small elements there and big elements where the field is smooth. Structured grids are uniform by construction.
Element-quality control. A skinny "sliver" triangle wrecks conditioning of the stiffness matrix and slows convergence. We want a way to guarantee well-shaped elements.

The algorithm that fixes all three at once is Delaunay triangulation. It produces the "fattest possible" triangulation of any point set — and it's the starting point of essentially every production mesher.

The empty-circle property

A triangulation is Delaunay if, for every triangle, the circumcircle (the unique circle passing through all three vertices) contains no other point from the set. That single property is what we test for, and what the triangulation algorithm enforces.

The remarkable fact: this property is equivalent to "maximize the minimum angle." Of all possible triangulations of a given point set, the Delaunay one has the largest minimum angle — i.e., the fewest skinny slivers. So choosing Delaunay isn't an arbitrary preference; it's what FEM quality wants.

import numpy as np

def in_circle(a, b, c, d) -> bool:
    """Returns True iff point d is strictly inside the circumcircle of
    counter-clockwise triangle abc. Standard 4x4 determinant predicate."""
    ax, ay = a[0] - d[0], a[1] - d[1]
    bx, by = b[0] - d[0], b[1] - d[1]
    cx, cy = c[0] - d[0], c[1] - d[1]
    det = ((ax**2 + ay**2) * (bx*cy - cx*by)
         - (bx**2 + by**2) * (ax*cy - cx*ay)
         + (cx**2 + cy**2) * (ax*by - bx*ay))
    return det > 1e-12

The in_circle test is a 4×4 determinant; you can derive it from "the circumcircle is the locus of points equidistant from $a$ , $b$ , $c$ " plus a sign trick. The interesting fact about it is that exactly this predicate, as a sign-of-determinant test, has been studied to death by computational-geometry people who care about robustness in the face of floating-point error. Production code uses Shewchuk's adaptive predicates; the textbook version above is fine for pedagogy.

The flip operation

Two adjacent triangles sharing an edge can always be retriangulated by flipping that edge — replacing the shared edge with the other diagonal of the four-vertex polygon they form. Of the two possible diagonals, exactly one yields a Delaunay-correct local triangulation (the empty-circle property is satisfied for both new triangles).

A Lawson flip algorithm starts from any triangulation and repeatedly flips edges that violate Delaunay until none remain. It always terminates. The resulting triangulation is uniquely Delaunay (assuming general position — no four points cocircular). This is one of the two main construction strategies; the other is incremental insertion.

Bowyer-Watson incremental insertion

The algorithm we'll actually implement. Start with a "super- triangle" big enough to contain every point you'll ever insert. Then insert points one at a time:

Find every existing triangle whose circumcircle contains the new point. Call them "bad."
The bad triangles form a connected region (a "cavity") whose boundary is the polygon we need to retriangulate.
Delete the bad triangles and connect each boundary edge to the new point — fan-triangulate the cavity.

After all points are inserted, drop any triangle that still touches the super-triangle vertices. The result is the Delaunay triangulation of your point set. The whole algorithm is short enough to fit on a screen:

def bowyer_watson(points: list[tuple[float, float]]) -> list[tuple[int, int, int]]:
    """Incremental Delaunay triangulation. Returns triangle vertex-index triples
    over the input points (super-triangle vertices are stripped before return)."""
    if len(points) < 3:
        return []

    # Super-triangle large enough to contain anything in [0, 1]^2.
    super_pts = [(-5.0, -3.0), (7.0, -3.0), (1.0, 9.0)]
    all_pts = list(points) + super_pts
    super_idx0 = len(points)

    triangles = [(super_idx0, super_idx0 + 1, super_idx0 + 2)]

    for pi, p in enumerate(points):
        # 1. Find every triangle whose circumcircle contains p.
        bad = [i for i, (a, b, c) in enumerate(triangles)
               if in_circle(all_pts[a], all_pts[b], all_pts[c], p)]
        if not bad:
            continue

        # 2. Find boundary edges of the cavity — edges of bad triangles that
        #    appear in only one of them. (Internal edges appear twice.)
        from collections import Counter
        def edge(u, v): return (u, v) if u < v else (v, u)
        counts = Counter(
            edge(t[k], t[(k + 1) % 3])
            for i in bad for t, k in [(triangles[i], 0), (triangles[i], 1), (triangles[i], 2)]
        )
        boundary = []
        for i in bad:
            t = triangles[i]
            for k in range(3):
                e = edge(t[k], t[(k + 1) % 3])
                if counts[e] == 1:
                    boundary.append((t[k], t[(k + 1) % 3]))

        # 3. Remove bad triangles, fan-triangulate the cavity from p.
        bad_set = set(bad)
        triangles = [t for i, t in enumerate(triangles) if i not in bad_set]
        for u, v in boundary:
            triangles.append((u, v, pi))

    # Drop any triangle touching the super-triangle.
    return [t for t in triangles if all(i < super_idx0 for i in t)]

Naïve runtime is $O (n^{2})$ — for each of $n$ points, we scan all current triangles. With a spatial index (location of the containing triangle in $O (lo g n)$ ) you can get to $O (n lo g n)$ , which is what production libraries do. For meshes up to a few thousand points, the naive version is fast enough.

See it

Click anywhere in the box to add a point. The triangulation rebuilds incrementally. Toggle "show circumcircles" and notice: no point ever lies inside any circumcircle. That's the empty- circle property, made visible.

add

presets

show circumcircles6 points · 5 triangles

Things to try:

Hit the 16 random button a few times. Each run produces a different point set, but every triangulation is well-shaped — no slivers, no points stranded inside a circumcircle. That's the maximize-minimum-angle property in action.
Hit ring of 12. The result is a fan from one of the points to all the others — Delaunay's choice for points on a circle, since any interior triangulation would put the missing center inside a circumcircle.
With circumcircles on, click in the middle of an existing large triangle. Watch which triangles light up red (their circumcircles contain your new point) and disappear; the new mesh fills in cleanly.
Add 100+ points by clicking + 20 random a few times. The triangulation stays valid; the algorithm really does scale.

From points to a mesh

Delaunay alone doesn't give you a finite-element mesh — it gives you a triangulation of an arbitrary point cloud, with no notion of "the boundary of my domain" or "I want elements no larger than 0.05 in this region." Two more pieces, both built on Delaunay, do the rest:

Constrained Delaunay respects specified edges (typically the domain boundary). The triangulation contains every constrained edge, even if that means slightly violating the empty-circle property locally. Implementation: insert all domain-boundary segments first, then fill the interior with Delaunay-respecting triangles.
Refinement (Ruppert's algorithm). While any triangle violates a quality threshold (minimum angle too small OR longest edge too long), insert a new point at its circumcenter and re-triangulate locally. Provably terminates with all triangles of bounded aspect ratio. This is what Triangle, Gmsh, and CGAL all do under the hood.

Combine the three — constrained Delaunay for the geometry, Ruppert's refinement for the quality, Bowyer-Watson at the core — and you have a real-world 2D mesher. 3D follows the same playbook with tetrahedra instead of triangles, but tetrahedral Delaunay is famously more delicate (slivers in 3D can have arbitrarily small radius-edge ratio while still being Delaunay-valid).

Adjacent algorithms worth knowing

Frontal advancing. Start from the boundary and march inward, generating one element at a time at the "front." Gmsh's default. Produces very high-quality elements near boundaries; harder to reason about than Delaunay.
Quadtree-based. Recursively subdivide the domain into cells until the local feature size is resolved, then triangulate within each cell. Triangle's "rough conforming" path. Fast, mediocre quality.
Hex meshing. 3D hexahedral elements give better convergence than tetrahedra for the same DOF count, but generating them automatically is famously hard. Most hex meshes are still produced by hand decomposition or sweep/extrusion of a 2D quad mesh.
Mesh adaptation. Given a solution and an error estimator, refine elements where the error is high (h-adaptivity), increase the polynomial order where the solution is smooth (p-adaptivity), or move existing nodes toward the action (r-adaptivity). All three live on top of whatever base mesher you start from.

Where this leaves us

The series has gone from "what is a weak form" to "here's a Delaunay mesher you can stand up against any production library, modulo robustness and speed." The three things that separate this code from FEniCS or deal.II are sparse storage, higher-order elements, and parallelism. None of those change the conceptual recipe. Every production FEM solver is, at its core, the same loop you've now written several times: assemble per-element contributions, apply boundary conditions, solve.