Mesh Generation: Algorithms & Methods

The initial classification of geometric domains is based on their dimensionality: two-dimensional (2D) or three-dimensional (3D). This fundamental distinction significantly influences the choice of mesh generation techniques and algorithms. Two-dimensional domains are further categorized into several specific types, each presenting unique challenges and opportunities for meshing. Understanding this initial dimensional classification is a critical first step in the mesh generation process, determining the type of algorithms and data structures most suitable for the task. The complexity of the meshing process increases with the dimensionality of the domain and the complexity of the shapes involved.

Four types of planar (2D) domains are defined. A simple polygon encompasses both its boundary and interior. A polygon with holes is a simple polygon excluding the interiors of other simple polygons, resulting in a boundary with multiple connected components. Multiple domains are more complex, featuring internal boundaries; effectively, they can be any planar straight-line graph where the infinite face is bounded by a simple cycle. These are useful for modeling objects composed of various materials. Finally, curved domains allow for boundaries defined by algebraic curves like splines, again potentially incorporating holes and internal boundaries. The first three categories are collectively referred to as polygonal domains, highlighting the fundamental geometric shapes involved in their representation. These distinctions are crucial for selecting appropriate mesh generation algorithms that can effectively handle the unique topological characteristics of each type of domain, particularly concerning boundary conditions and internal structures.

The categorization of geometric domains is directly relevant to the modeling of real-world objects and materials. Multiple domains, for example, are expressly designed to represent objects made from more than one material. Each material might require a different meshing strategy or refinement level, due to variations in material properties or behaviors. The accurate representation of these domains is essential for the successful simulation of physical phenomena within the object. The complexity of the domain, in terms of holes, internal boundaries, and curved surfaces, poses significant challenges in terms of the efficiency and accuracy of the mesh generation algorithms. The choice of an appropriate algorithm will be influenced not just by the complexity of the geometry, but also by the need to maintain fidelity in representing the material properties and their spatial distribution within the modeled object.

The choice of mesh generation algorithm is directly influenced by the type of geometric domain. For instance, simple polygons might lend themselves to straightforward triangulation methods, while dealing with polygons containing holes or multiple domains requires more sophisticated algorithms capable of managing complex topological information. Similarly, curved domains necessitate techniques that can approximate curves effectively to produce a mesh that faithfully represents the curved boundaries and their interactions with adjacent elements. Algorithms tailored to specific domain types can improve mesh quality (measured through aspect ratios, angles, etc.) and computation efficiency. Selecting the right algorithm is critical to obtaining accurate and reliable simulation results, as mesh quality is directly linked to the accuracy and computational cost of the simulation itself. Different algorithms have varying computational complexities, and the optimal choice depends on the complexity of the geometry and the desired level of accuracy. The initial classification of the geometric domain provides a crucial foundation for choosing and implementing a suitable mesh generation algorithm.

The core of mesh generation lies in its application to solving partial differential equations, commonly using finite element (FEM) or finite difference (FDM) methods. FEM typically employs unstructured meshes, offering flexibility in handling complex geometries. In contrast, FDM often utilizes structured meshes due to their simplicity and efficiency in computation, particularly when using finite difference stencils. The choice between structured and unstructured meshes depends significantly on the specific problem and the desired balance between computational efficiency and geometric complexity. The inherent characteristics of each method and the suitability of various mesh types (e.g., triangular, quadrilateral, tetrahedral, hexahedral) to different types of problems must be carefully considered. Structured meshes, for example, while easier to implement, might require significantly more elements than unstructured meshes for the same problem, particularly when dealing with rapidly changing element sizes. The selection of the optimal mesh type and method is highly problem-dependent.

Domain decomposition presents a hybrid approach, dividing the problem domain into smaller, potentially overlapping subdomains. Subproblems are solved directly within these smaller domains, then iteratively combined to reach a global solution where neighboring subdomain solutions agree. This technique combines the advantages of direct solvers within smaller regions with the superior convergence properties of iterative solvers over the entire domain. It is particularly valuable for complex geometries where a global structured mesh might be impractical or computationally expensive. This approach offers a compromise, leveraging the strengths of structured meshes for subdomains while reducing the overall burden on the mesh generator. The partitioning of the domain allows for the use of efficient structured mesh generation techniques for each subproblem, simplifying the generation process while still allowing for an accurate and robust solution of the overall problem, mitigating some of the computational difficulties associated with larger, complex domains.

Solving the sparse linear system arising from either FEM or FDM is computationally the most demanding step in the mesh generation process. Direct factorization methods, such as sparse Cholesky or LU factorization, are robust but can be slow, especially for large systems. Their computational cost can be amortized if the factorization is reused for multiple right-hand sides. Iterative methods, although potentially faster, often require preconditioning to improve convergence, and their efficacy is heavily influenced by the mesh generation process and discrete formulation choices. A popular preconditioning technique uses an incomplete factorization, balancing computational cost against accuracy. The choice between direct and iterative methods involves a trade-off between speed and robustness, influenced by factors such as the size of the system, the mesh structure, and the desired accuracy of the solution. Higher-order basis functions in FEM can lead to smaller but denser linear systems, impacting the computational cost and the effectiveness of different solution methods.

The mesh generation step significantly impacts the effectiveness of solution methods. For instance, using higher-order basis functions in FEM allows for coarser meshes, resulting in smaller, denser linear systems that may be solved more efficiently by direct methods. Conversely, the finite difference method often relies on structured meshes for better accuracy approximation with orthogonal edges. The characteristics of the mesh—such as element shape, aspect ratio, and structure—directly influence the condition number of the linear system, affecting both the speed and stability of the solver. For example, poorly shaped elements can lead to slow convergence or even instability in iterative methods. Therefore, optimizing the mesh generation process for the chosen numerical method is crucial for overall efficiency and solution accuracy. Careful consideration must be given to the interplay between the chosen numerical method (FEM or FDM) and the mesh properties to minimize computational time and maximize the accuracy of the final solution.

Mesh quality is intrinsically linked to the aspect ratio of mesh elements. Large aspect ratios, where elements are excessively elongated or flattened, can lead to unacceptable interpolation errors, compromising the accuracy of numerical solutions. Early research demonstrated convergence as triangular elements shrink, provided angles remain bounded away from 0°. Later analysis relaxed this constraint, showing convergence if angles are bounded away from 180°. This concept extends to 3D, where dihedral angles must be bounded away from 0°. Bank and Smith's triangle quality measure, defined as the triangle's area divided by the sum of squared edge lengths, offers a practical metric, slightly favoring sharp over obtuse triangles. The influence of aspect ratio on interpolation error is a crucial factor determining the accuracy of the numerical solution. Large aspect ratios can lead to significant inaccuracies, while well-shaped elements with controlled aspect ratios are essential for obtaining reliable results.

While typically undesirable, high aspect ratio elements can be beneficial for anisotropic solutions, where the second derivative of the solution varies greatly with direction. Properly aligned high-aspect-ratio elements provide a highly efficient mesh in these cases. The ideal aspect ratio of a triangle is the square root of the ratio of the largest to smallest eigenvalues of the Hessian. For triangular meshes, the presence of large angles is more problematic than small angles when aspect ratios exceed the ideal. The alignment and aspect ratio of elements must be carefully considered in relation to the characteristics of the problem's solution. The optimal shape depends on the direction of the highest variation in the solution's derivatives, necessitating optimized element alignment for best results. The choice of element shape is not just about avoiding poor quality, but also about leveraging element properties to best represent the physical phenomenon being modeled.

Element shape influences the properties of the resulting linear system, beyond just its condition number. A triangular mesh with well-shaped elements yields a symmetric M-matrix—positive definite with negative off-diagonal entries—for a finite element formulation involving a Laplacian operator. M-matrices satisfy a discrete maximum principle, preventing numerical oscillations in the solution. Fluid flow problems, especially Navier-Stokes simulations, are strongly anisotropic; ideal aspect ratios in aerodynamics can reach 10,000. Quadrilateral and hexahedral meshes offer accuracy advantages over triangular and tetrahedral meshes for control-volume formulations of such problems, allowing element faces in boundary layers to be nearly parallel or orthogonal to the surface. The mathematical properties of the linear system arising from the mesh discretization are directly linked to the element shapes. Well-shaped elements lead to favorable properties in the linear system, improving the stability and efficiency of numerical solution techniques. In contrast, poorly shaped elements can negatively affect the solution's stability and lead to inaccuracies.

Tetrahedra can be poorly shaped in more ways than triangles. In 2D, problems arise from angles near 0° or 180°, but in 3D, classifications consider both dihedral and solid angles, resulting in five types of poorly shaped tetrahedra: needles, wedges, slivers, caps, and spires. A sliver, for instance, can have all face angles bounded away from 0° and 180°, yet still have arbitrarily small solid angles and volume. The same domain can be tetrahedralized with varying numbers of tetrahedra, leading to complexities in mesh optimization. Edge flipping, a technique used to improve mesh quality, has a more complicated generalization to 3D, making mesh improvement and refinement significantly more challenging in three dimensions than in two. The more complex geometric relationships in 3D require more sophisticated algorithms and techniques to achieve and maintain the desired mesh quality and avoid inaccuracies. The identification and removal of these poorly shaped elements are critical for achieving a robust and accurate numerical solution.

Structured meshes offer significant advantages in terms of simplicity and efficiency. They require considerably less memory than unstructured meshes with the same number of elements because neighbor connectivity is implicitly defined through array storage. This implicit definition allows for faster access to neighboring cells during computations, leading to optimized code and efficient execution, particularly on vector machines. The simpler data structure and inherent regularity of structured meshes translate into faster computation times and reduced memory requirements, making them highly attractive for certain applications. However, the inherent regularity can limit their applicability to complex geometries, as it can become difficult or impossible to create a structured mesh for domains with intricate shapes. This limitation often necessitates the use of unstructured meshes or hybrid approaches.

A conformal mapping, which preserves angles, provides a powerful approach for generating structured meshes. The Riemann mapping theorem guarantees the existence of a conformal mapping between any topological disk and another, such as a unit disk or a square. Mapping a region onto a square grid induces a structured mesh where element angles tend towards 90° with increasing discretization refinement. While offering advantages like near-orthogonal elements, true conformal mappings are not widely used in mesh generation due to the computational cost of obtaining such mappings and the lack of local control over point spacing. While this method offers highly desirable mesh characteristics such as near-orthogonal elements, the computational cost associated with calculating true conformal mappings and the limitations on local control over point spacing have restricted its widespread use. Other methods have been developed to approximate the benefits of conformal mappings, offering a balance between computational cost and mesh quality.

Delaunay triangulation is a cornerstone of unstructured mesh generation. Several approaches exist for generating unstructured meshes based on Delaunay triangulation. One involves combining vertices from multiple structured meshes. Another adds Steiner points layer by layer, working inward from the domain boundary – a technique similar to advancing front mesh generation. Random point scattering with filtering represents a third approach, where points are scattered randomly according to some distribution, followed by filtering to remove closely spaced points and improve the overall quality. Deterministic methods achieve a similar effect to random sampling with filtering by applying birth and death rules that depend on the density of neighboring points. These diverse methods highlight the versatility of Delaunay triangulation for meshing complex shapes, offering various trade-offs between control over element shape, ease of fitting complex geometries, and computational cost.

The generation of a conforming Delaunay triangulation, where the triangulation conforms to the boundary of the domain, is a crucial step in many 2D meshing algorithms. Achieving conformity may require adding Steiner points to the boundary. Efficient algorithms like Saalfeld's algorithm add uniformly spaced Steiner points (except near endpoints) along edges, ensuring that edges of the domain are represented as unions of edges in the Delaunay triangulation. Subsequent steps improve the mesh by adding Steiner points based on angle criteria, aiming to produce a mesh with fewer, rounder triangles. Algorithms like Ruppert's algorithm combine constrained Delaunay triangulation with Steiner point insertion to achieve a balance between mesh quality and efficiency, minimizing the total number of triangles while maintaining a specified minimum angle. These approaches ensure that the resulting mesh accurately captures the domain's boundary and has acceptable element shapes, paving the way for accurate numerical solutions.

Three-dimensional mesh generation often employs Delaunay triangulation, similar to the 2D case. Point placement methods, including combining structured meshes, advancing front techniques, and random scattering with filtering, are used. A crucial aspect is ensuring that the Delaunay triangulation conforms to the domain boundary, requiring sufficient points on the boundary. While the 3D conforming Delaunay triangulation problem is manageable for most practical domains, published solutions remain less common than in 2D. The same challenges of improper point spacing at junctures between fronts or patches persist from the 2D case. Additionally, even well-spaced point sets may contain sliver tetrahedra in their Delaunay triangulation because a sliver's circumsphere is not unusually large compared to its edge lengths. Therefore, post-processing steps are frequently included to detect and remove these problematic elements. The ANSYS Inc. advancing front mesh generator is mentioned as an example of a system that places elements directly, rather than using a point placement followed by Delaunay triangulation approach.

An algorithm by Marcum and Weatherill combines advancing front methods with Delaunay triangulation and iterative improvement. It starts with a coarse mesh and adds Steiner points using advancing front, subdividing tetrahedra to maintain triangulation. The mesh is then improved using Delaunay flips followed by minmax-solid-angle flips, a strategy found to be more effective than either type of flip alone. Pure advancing front methods iteratively choose a face from the initial boundary faces and build a tetrahedron over it. New vertices, when needed, are strategically placed to minimize collisions and to account for aspect ratio considerations. The algorithm also prioritizes filling clefts to ensure overall mesh integrity, showcasing how mesh generation strategies in 3D demand a delicate balance between topology and geometry. The need for careful consideration of element shapes and the use of iterative improvement techniques such as flipping and smoothing are highlighted in this section.

Octrees, the 3D generalization of quadtrees, provide a hierarchical approach to mesh generation. An initial bounding cube is recursively subdivided into eight congruent cubes until each minimal cube intersects the domain in a simple way. A balance condition prevents adjacent cubes from having vastly different sizes. Refinement in 3D is discussed through the generalization of bisection to tetrahedra. Bisecting a tetrahedron across an edge creates two child tetrahedra. Maintaining mesh validity requires consistent bisection of shared faces between adjacent tetrahedra. A bisection algorithm by Bänsch generates a finite number of similarity classes by mapping tetrahedra to a canonical tetrahedron for which longest-edge bisection has a finite number of similarity classes. The algorithm ensures that the resulting tetrahedra belong to a limited set of similarity classes, reducing complexity. Consistent bisection of shared faces between adjacent tetrahedra is essential to prevent the creation of invalid meshes, emphasizing the topological constraints in 3D mesh refinement.

Three-dimensional mesh improvement often employs edge flipping, generalized from the 2D case. Flipping is performed first based on Delaunay's empty sphere criterion and then by minimizing the maximum solid angle. Laplacian smoothing, another common technique, is also applicable in 3D, but it doesn't guarantee improved element quality and may require explicit checks to prevent element inversion. The discussion highlights the challenges in 3D mesh refinement, noting that fewer proven techniques exist compared to 2D. Bänsch's bisection algorithm is presented as an example of a refinement method that generates a finite number of similarity classes, maintaining a controlled set of tetrahedral shapes throughout the refinement process. The iterative nature of these mesh improvement techniques (such as edge flipping and Laplacian smoothing), combined with the topological constraints present in 3D, makes the generation of high-quality meshes significantly more complex in three dimensions than in two.