mesh-generator

How it is built

Architecture overview

A mesh generator turns a geometry (a CAD solid, or a scanned surface) and a size specification into a mesh: a partition of the domain into simple cells (triangles and tetrahedra, sometimes quadrilaterals and hexahedra) that fit the boundary exactly, are small where the physics needs resolution and large elsewhere, and are well shaped, because solver accuracy and conditioning degrade on thin or flat cells. The mesh then feeds a finite-element or finite-volume solver, which is why the output includes numbering, element groups, boundary tags and often a partition for parallel runs.

The one idea that makes it tractable is meshing by dimension: place vertices on the model's corners, mesh each curve with edges, mesh each surface with triangles that conform to the curve meshes, then mesh each volume with tetrahedra that conform to the surface meshes. Every stage works on a simpler problem than the one above it, and conformity is enforced by construction because a lower-dimensional mesh is the fixed boundary of the next. A sizing field (a function from position to desired edge length) steers every stage, and a final quality loop (swap, collapse, split, smooth) repairs what the first pass left behind. Delaunay triangulation is the workhorse of the core because its empty-circumsphere rule gives good shapes in 2D and a canonical answer; it needs robust geometric predicates and a step that recovers the prescribed boundary.

Parts:

• Geometry model and kernel — solids as boundary representations with curved faces; intersection and Boolean operations for building the model — B-rep topology (no page yet), NURBS (no page yet), Surface-surface intersection (no page yet), Solid Boolean operations (no page yet), Newton root finding, Interval arithmetic (no page yet)
• Robust arithmetic — decisions about orientation and circumspheres that never lie — Exact geometric predicates (no page yet), Floating-point error
• Spatial queries — nearest geometry, point location, size lookups, overlap tests — Bounding volume hierarchy, k-d tree, R-tree, Quadtree and octree (no page yet), Space-filling curves (no page yet)
• Triangulation core — incremental Delaunay insertion, cavity retriangulation, dual diagrams — Delaunay triangulation (no page yet), Voronoi diagrams (no page yet), d-ary heap, PRNGs (PCG, xoshiro)
• Surface and polygon meshing — triangles on curved faces and flat polygons with holes — Polygon tessellation (no page yet), Half-edge mesh (no page yet)
• Mesh storage — cell-face-node connectivity for solvers — Unstructured mesh representation (no page yet), Sparse matrices (CSR)
• Partitioning and ordering — balanced parts and bandwidth-reducing numbering for the solver — Graph partitioning, Fill-reducing ordering

 CAD file (STEP/BREP) ─► geometry model ──────────────┐   size specification
                         (vertices, curves, surfaces,  │   (points, curvature,
                          volumes; NURBS)               │   distance fields)
                               │                        ▼
                               │               ┌──────────────────┐
                               │   queries     │ sizing field h(x)│◄─ k-d tree / R-tree /
                               ▼ ◄──────────── └────────┬─────────┘   octree lookups
        1D: edges along curves  ───────────────────────┤
        2D: triangles in (u,v), mapped to 3D            │ Delaunay insertion
        3D: tetrahedra, boundary recovery               │ + robust predicates
                               │                        │
                               ▼                        │
                    quality loop: flip, split, collapse, smooth
                               │
                               ▼
              mesh data structure ─► renumber ─► partition ─► solver files

How the parts fit together

Geometry model and kernel

The input is a hierarchy: vertices bound curves, curves bound surfaces, surfaces bound volumes (B-rep topology, no page yet). Curves and surfaces carry geometry, usually as splines (NURBS, no page yet), so the mesher can ask two questions of any entity: evaluate a point from parameters, and find the closest point on the entity to a given point. The second is a root-finding problem solved with Newton root finding, started from a nearby guess; a Bounding volume hierarchy over the faces gives the guess and discards the faces that cannot be near. Building a model from primitives (a plate with holes, a pipe joining a tank) needs Solid Boolean operations (no page yet), which in turn need Surface-surface intersection (no page yet) to find where two curved faces meet; Interval arithmetic (no page yet) gives rigorous bounds for isolating those intersection curves. Real CAD files have gaps, overlaps and tolerance mismatches, so a mesher normally includes a repair or tolerance-merge step; the geometry kernel (Gmsh uses OpenCASCADE for this) owns that work, and the mesher consumes a clean topology.

The contract with the meshing stages: every entity can be sampled by parameter, queried for its nearest point and its curvature, and has a well-defined boundary made of lower-dimensional entities.

Sizing field

A mesh needs to say "edges of length 0.5 mm here, 20 mm there". Gmsh takes the minimum of several sources: sizes given at geometric points and interpolated through the model, a size derived from curvature (a stated number of elements per 2π of turning angle), and general size fields such as boxes, distance-to-entity fields with a transition law, and a background mesh that stores sizes at nodes (Gmsh manual, "Specifying mesh element sizes"). Nearest-entity lookups for distance fields use spatial indices: a k-d tree for point sets, an R-tree for boxes of curves and faces, and a Quadtree and octree (no page yet) when the size field itself is stored hierarchically. For an adaptive loop in which a solver's error estimate produces the next size field, the field is defined on the previous mesh and must be interpolated at the points of the new one.

Triangulation core

Delaunay triangulation (no page yet) is built by incremental insertion (Bowyer 1981, Watson 1981): find the simplex that contains the new point, collect every simplex whose circumsphere contains it (the cavity, a connected star-shaped region), delete them, and connect the point to the boundary of the cavity. The empty-circumsphere test and the point-in-simplex test are determinants of coordinates; their sign decides the combinatorics, so a wrong sign from rounding can corrupt the structure. That is why Exact geometric predicates (no page yet) are part of the core: the standard approach (Shewchuk 1997) first evaluates the determinant in floating point with a rigorous error bound and falls back to higher precision only when the sign is in doubt, so the common case costs little (Floating-point error explains why the filter is needed). The dual of the triangulation is the Voronoi diagram (no page yet), whose vertices are the circumcenters used as candidate points for refinement.

Insertion order matters for speed. Random order gives good expected behaviour for point location; sorting points along a Space-filling curve (no page yet) gives cache locality and short walks from one insertion to the next; the BRIO scheme (biased randomized insertion order, Amenta, Choi and Rote 2003) combines both and is also what makes the parallel version feasible. Randomness uses a small generator (PRNGs (PCG, xoshiro)). Refinement processes the worst cells first through a priority queue (d-ary heap) keyed by quality.

Surface and polygon meshing

For a curved face, the usual route is to mesh in the parameter plane. The 1D meshes of the bounding curves become a polygon in (u,v); a 2D Delaunay mesh is built there, and the points are mapped to 3D by evaluating the surface. Because the parametrization stretches lengths, the sizes must be converted to a metric in parameter space (from the surface's first fundamental form); ignoring this gives skinny triangles in 3D. Gmsh describes its 2D approach concretely: an initial Delaunay mesh containing all the points of the 1D mesh is built by divide-and-conquer, missing edges are recovered by edge swaps, and then one of several algorithms improves it: "MeshAdapt" (local edge swaps, splits and collapses, the most robust on very complex curved surfaces), "Delaunay" (insert at the circumcenter of the element with the largest dimensionless circumradius) or "Frontal-Delaunay" (placing points from an advancing front, best element quality). For flat polygons with holes, as in rendering and some solver pre-processing, Polygon tessellation (no page yet) is a lighter-weight alternative. Editing surface meshes (flips, collapses) is easiest on a Half-edge mesh (no page yet), which gives constant-time access to the neighbours of an edge.

Volume meshing

Three families exist and generators combine them:

• Delaunay with boundary recovery. Mesh the convex hull of the surface vertices, then force the prescribed surface triangles to appear as faces of the tetrahedra (boundary recovery by flips and local point insertion, as in TetGen, Si 2015), delete the tetrahedra outside, and insert interior points to meet the size field. Gmsh's 3D "Delaunay" does exactly this using the TetGen boundary-recovery algorithm; "HXT" is a parallel reimplementation (Marot, Pellerin and Remacle 2018, "One machine, one minute, three billion tetrahedra").
• Advancing front (Gmsh's "Frontal" uses the Netgen algorithm): start from the boundary triangles and grow tetrahedra inward; good control of element size and boundary layers, harder to make robust at the end where fronts meet.
• Octree-based: a recursive subdivision of a bounding box refined to the size field, cut cells fitted to the boundary. Fast and robust, and the natural base for adaptive refinement of hierarchical forests (p4est), at the cost of weaker boundary conformity.

Delaunay in 3D has a known flaw: it can create slivers, tetrahedra with four vertices almost on a circle, which satisfy the Delaunay rule and have good circumradius-to-edge ratios but almost zero volume. Slivers are removed by a post-processing pass of flips, vertex smoothing and perturbation, which is why the quality loop is not optional.

Mesh storage

Solvers need cell-face-node connectivity in compact arrays: for each cell its nodes, for each face its two neighbouring cells, and boundary tags (Unstructured mesh representation, no page yet). The same adjacency is the sparsity pattern of the stiffness matrix, stored in Sparse matrices (CSR) form. The mesher's internal structure (with neighbour pointers for fast cavity search) is converted to this solver-facing form on output.

Partitioning and ordering

For parallel solvers the mesh graph (cells as vertices, shared faces as edges) is split into balanced parts with few cut edges by Graph partitioning (multilevel schemes such as METIS: Karypis and Kumar 1998). For the linear solver on each part, nodes are renumbered with a Fill-reducing ordering to shrink the bandwidth or the fill of a sparse factorization. Both operations are pure graph algorithms on the connectivity from the previous section.

Data and control flow

Mesh a CAD solid (Gmsh-style pipeline):

1. Import the CAD model, merge near-coincident entities, build the topology tree (vertices, curves, surfaces, volumes).
2. Evaluate the sizing field at the model's points and from curvature; combine with user fields (minimum rule).
3. Mesh each curve: place vertices along it so the integral of 1/h along the arc matches the number of edges.
4. Mesh each surface: map the curve meshes into parameter space, build the initial Delaunay mesh, recover boundary edges, refine to the size field, map back to 3D. Time goes to point insertion and predicates; for large surfaces to the cavity searches.
5. Mesh each volume: tetrahedralize the surface-mesh vertices, recover the surface triangles, remove outside cells, insert interior points to meet sizes.
6. Run the quality loop (flips, smoothing, optional sliver removal) until the worst-element quality is above a threshold or iterations run out.
7. Number nodes and cells, assign physical-group tags, write the mesh; optionally partition.

Insert one point into a Delaunay mesh (Bowyer-Watson):

1. Locate a containing simplex by walking from the last inserted point's neighbourhood (orientation tests).
2. Starting from it, flood to neighbours whose circumsphere contains the point (insphere tests); this set is the cavity.
3. Delete the cavity's simplices, create one new simplex per cavity boundary face with the new point, and set the neighbour links.
4. Cost is proportional to the cavity size; with a good insertion order it is a small constant on average.

An adaptive solve-remesh loop:

1. Solve on the current mesh and compute a per-cell error estimate.
2. Turn the error into a new sizing field (smaller where the error is high).
3. Remesh from scratch with the new field, or refine and coarsen locally (octree forest) while keeping the topology.
4. Transfer the solution to the new mesh; repeat until the estimate falls below tolerance.

Scaling

In the order you would pull the levers:

1. Predicate cost. Use filtered arithmetic so exact evaluation is rare; the filter makes this nearly free in practice.
2. Point location and cavity search. Good insertion order (BRIO or a space-filling-curve sort), walking with a locality hint, compact cell storage with neighbour indices.
3. Single-thread throughput. Cache-friendly layout of vertices and cells, avoid allocation per cell, process the quality queue in batches.
4. Shared-memory parallelism. Partition space, insert points in different regions concurrently with conflict checks on cavities (HXT's approach, with BRIO rounds), then merge the borders.
5. Memory. Billions of cells need compact indices and out-of-core or distributed storage; generate the mesh already partitioned, with each rank meshing its own block, or use an octree forest distributed across ranks.
6. Fewer cells. The cheapest speed-up is a sizing field that is only as fine as the physics needs; graded transitions beat uniform fine meshes, and anisotropic cells (stretched in the boundary layer) cut counts again.

Failure modes

Real implementations

Sources

• C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Meth. Engng 79(11), 2009. doi; Gmsh reference manual (algorithms, size fields). manual
• C. Marot, J. Pellerin, J.-F. Remacle, One machine, one minute, three billion tetrahedra, Int. J. Numer. Meth. Engng, 2018. doi
• A. Bowyer, Computing Dirichlet tessellations, The Computer Journal 24(2), 1981. doi; D. F. Watson, Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes, The Computer Journal 24(2), 1981. doi
• J. R. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete Comput. Geom. 18, 1997. doi
• J. R. Shewchuk, Delaunay refinement algorithms for triangular mesh generation, Comput. Geom. 22, 2002. doi
• H. Si, TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator, ACM Trans. Math. Softw. 41(2), 2015. doi
• N. Amenta, S. Choi, G. Rote, Incremental constructions con BRIO, SoCG 2003. doi
• G. Karypis, V. Kumar, A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs, SIAM J. Sci. Comput. 20(1), 1998. doi
• C. Burstedde, L. C. Wilcox, O. Ghattas, p4est: Scalable Algorithms for Parallel Adaptive Mesh Refinement on Forests of Octrees, SIAM J. Sci. Comput. 33(3), 2011. doi