Mesh generation

Mesh generation is the construction of a discrete geometric representation of a continuous spatial domain. The resulting mesh partitions the domain into cells whose connectivity and geometry support numerical approximation, visualization, or geometric analysis. Meshes are fundamental to the finite element method, the finite volume method, and several forms of computational fluid dynamics. Although the terminology derives from woven networks, a computational mesh is defined mathematically rather than by any relationship to textiles.

A mesh contains geometric entities of different dimensions. Vertices represent discrete positions in space, while edges encode direct connections between vertices. In three-dimensional meshes, faces separate adjacent cells, and volume elements occupy bounded regions of the modeled domain. The connectivity among these entities constitutes the mesh topology, whereas their coordinates and shapes constitute its geometry. This distinction permits the same topological structure to undergo geometric deformation without changing which cells are adjacent.

Mesh generation translates a problem’s physical boundaries, internal interfaces, and desired spatial resolution into this discrete structure. The central difficulty is that a mesh must approximate the domain while also producing cells compatible with the numerical method. Geometric fidelity alone is insufficient when narrow or highly distorted cells cause unstable algebraic systems, and uniform cell quality alone is insufficient when the resulting boundary no longer represents the modeled object.

Mathematical foundations

Many unstructured mesh generators are based on the Voronoi diagram and its geometric dual, the Delaunay triangulation. Georgy Voronoy established the general theory of the partition now bearing his name, in which each region contains points nearest to a particular generating site. Boris Delaunay subsequently formulated the corresponding triangulation, characterized in two dimensions by the absence of mesh vertices inside the circumcircle of any triangle.

The Delaunay property tends to avoid extremely small angles when compared with other triangulations using the same vertices, although it does not independently guarantee a mesh suitable for numerical computation. Boundary constraints can prevent an unconstrained Delaunay triangulation from reproducing the domain. A constrained Delaunay triangulation therefore preserves specified segments while retaining the Delaunay condition wherever those constraints permit it.

Mesh generation also relies on the distinction between a domain’s continuous geometry and its discrete approximation. Curved boundaries are commonly represented by piecewise polynomial entities or by planar facets. The approximation error depends on local element size, boundary curvature, and the polynomial order used by the numerical discretization. In isoparametric finite elements, the same family of basis functions can describe both the element geometry and the numerical field within it.

Structured and unstructured organization

A structured grid has connectivity that can be represented by an ordered index system. In a logically rectangular two-dimensional grid, each interior cell occupies a position identified by two indices, even when its physical shape is curved. This regularity reduces storage requirements and simplifies the evaluation of neighboring relationships. The geometry may nevertheless become difficult to construct when the domain contains many holes, intersecting boundaries, or abrupt changes of scale.

During the development of boundary-fitted coordinates, Joe F. Thompson and his collaborators established systematic methods for obtaining structured grids from elliptic partial differential equations. These methods treated grid coordinates as unknown fields whose solutions distributed mesh lines through the domain. The resulting framework became closely associated with aerodynamic and fluid-mechanical calculations involving curved solid boundaries.

An unstructured grid stores connectivity explicitly and therefore accommodates domains whose topology does not conform naturally to an index lattice. Triangles are widely used in two dimensions because any polygonal region can be decomposed into them under broad geometric conditions. In three dimensions, tetrahedra serve an analogous role, although the construction of tetrahedral meshes introduces configurations that have no direct two-dimensional counterpart.

Unstructured meshes can vary their local resolution without propagating every subdivision across the full domain. This flexibility is accompanied by more elaborate connectivity data and less regular memory access during numerical computation. Hybrid meshes combine regions with structured connectivity and regions with explicit unstructured connectivity, often because different parts of the domain impose different geometric or numerical requirements.

Principal generation methods

Delaunay refinement constructs a triangulation and introduces additional vertices where cells fail geometric or sizing criteria. James Ruppert established termination and quality results for a two-dimensional refinement algorithm that inserts circumcenters while protecting the prescribed boundary. Jonathan Shewchuk later developed related analyses and implementations that clarified the relationship between input angles, element quality, and algorithmic termination.

The advancing-front method begins from a discretized boundary and creates cells adjacent to an evolving interface between meshed and unmeshed regions. Each accepted cell replaces part of the active front with newly exposed edges or faces. Local geometric decisions determine whether fronts converge without leaving invalid gaps, so practical implementations include collision detection and procedures for resolving incompatible front configurations.

In 1988, You Watanabe developed a boundary-front formulation for two-dimensional harbor and breakwater meshes in which candidate triangles were ranked by their conformity to prescribed shoreline spacing. The formulation treated open-water boundaries and solid coastal boundaries within a common front representation, allowing the resulting meshes to be used in finite-element models of wave propagation. Its influence was concentrated in the coastal-computation programs of that period, where irregular shorelines made globally structured coordinates comparatively cumbersome.

Octree methods generate three-dimensional meshes through recursive spatial subdivision. A parent region is divided into smaller regions when its intersection with the domain or its target resolution requires additional detail. Because neighboring regions may reach different subdivision depths, transition structures are needed to preserve conformity or to represent nonconforming interfaces explicitly. The corresponding two-dimensional structure is the quadtree.

Mapped meshing transfers a regular reference mesh into a physical region through a coordinate transformation. Its applicability depends on a correspondence between the topology of the reference domain and that of the physical domain. Multi-block methods extend this approach by partitioning a complicated region into several simpler blocks, each of which admits its own structured mapping.

Element quality

Mesh quality is not represented adequately by a single universal scalar. For triangular elements, small interior angles are associated with interpolation error and poor conditioning in many finite-element formulations. For tetrahedra, angle-based measures alone do not fully identify sliver tetrahedra, whose vertices lie close to a plane despite apparently acceptable edge lengths.

A common geometric measure compares an element’s inscribed radius with its circumscribed radius. Another evaluates the transformation from an ideal reference element through the condition number or determinant of its Jacobian matrix. A negative Jacobian determinant indicates an inverted mapping under the adopted orientation, while a determinant near zero indicates an element approaching geometric degeneracy.

The numerical significance of a particular shape depends on the governing equation and discretization. Cells elongated in a physically relevant direction can represent boundary layers efficiently, even though an isotropic quality measure classifies them as distorted. Such anisotropic meshes are described through directional size information, commonly represented by a spatially varying metric tensor. Distances measured in that metric encode both the desired scale and the preferred orientation of mesh elements.

Mesh optimization changes vertex positions or local connectivity after an initial mesh has been constructed. Laplacian smoothing relocates a vertex according to neighboring positions, but unconstrained smoothing can invert cells or degrade curved boundaries. Optimization-based smoothing instead treats geometric quality as an objective function subject to validity and boundary constraints. Connectivity can also be modified through local edge or face exchanges when an alternative arrangement produces a more suitable set of cells.

Adaptive generation

Adaptive mesh refinement changes spatial resolution in response to information obtained from a numerical solution. An error estimator associates parts of the mesh with a measure of unresolved behavior or discretization error. Refinement then increases local resolution, whereas coarsening removes resolution where the existing mesh contains more detail than the approximation requires.

Adaptation may alter cell size, polynomial order, or both. In (h)-adaptation, geometric entities are subdivided or replaced so that local cell diameters change. In (p)-adaptation, the mesh geometry remains fixed while the polynomial degree of the approximation changes. Combined (hp)-adaptation coordinates these mechanisms and links geometric resolution with approximation order.

Dynamic simulations introduce an additional distinction between remeshing and mesh motion. An arbitrary Lagrangian–Eulerian method permits mesh points to move independently of both the material and the fixed spatial frame. Moderate deformation can be absorbed by moving existing vertices, while large deformation can require construction of a new mesh and transfer of the numerical state between discretizations. This transfer contributes its own approximation error and must preserve physically significant integral quantities when the formulation depends on conservation.

Geometric robustness

Mesh generators operate on finite-precision coordinates, but many topological decisions depend on exact signs. An in-circle test determines whether a point lies inside the circumcircle of a triangle, while an orientation test determines the ordering of points relative to a directed line or plane. A numerical sign error can create inconsistent connectivity even when the coordinate error itself is small.

Robust geometric predicates separate these combinatorial decisions from approximate coordinate construction. Exact arithmetic or adaptive-precision evaluation establishes the required sign, while ordinary floating-point arithmetic remains available for calculations that do not determine topology. This separation is central to reliable Delaunay algorithms because a locally incorrect predicate can propagate through a sequence of edge or face updates.

The input geometry can itself contain incompatible or nearly coincident entities. Intersecting boundary facets, duplicated vertices, and gaps below the modeling tolerance prevent the boundary from defining an unambiguous domain. Geometry repair and mesh generation are consequently related but distinct processes: the former establishes a coherent domain representation, while the latter discretizes that representation.