Generation of Triangular and Tetrahedral Meshes

The tree-based approach [6,9,15,28,38,36,44,46,57,68,94,128] uses a spatial decomposition procedure representing the domain as a set of disjoint cells stored in a hierarchic tree. The domain to be meshed is placed in a square or cubic cell that entirely encloses it. This cell is called the root of the tree. The structure of the tree is then defined in a recursive manner by subdividing the individual cells into four quadrants (quadtree in 2D) or eight octants (octree in 3D) until all cells are at the desired level, dictated by the geometry of the domain and by the mesh size specification. Since the tree structure may be highly irregular, a one-level difference is enforced. This ensures that cells sharing an edge do not differ in the size by a factor greater than two and that a gradual variation of the final mesh density is guaranteed. After that, each cell is classified to be inside, outside, or boundary with respect to the domain being discretized (some inside cells may be later reclassified as boundary cells if they are too close to the boundary). No further attention is paid to the outside cells. The inside cells are meshed using templates, predefined sets of elements filling the entire cell. Since the one-level difference rule has been applied, which effectively limits the number of midside and midface nodes to one per cell edge and face, respectively, the overall number of templates is significantly reduced. While the 16 templates in the 2D case can be easily created manually, the 6210 templates (4096 when only midside nodes are considered, 2114 additional templates when also midface nodes are considered) in the 3D case can hardly be implemented manually even if rotations and mirroring are taken into account. This problem has been resolved by manual implementation of six most frequently occurring templates that cover over 90 % of all cases and by application of a simple but fast procedure in the remaining cases. The discretization of boundary cells is usually more difficult because of requirements on a consistent and unambiguous way of treating the domain boundary. In the early versions [6,9], only a limited number of specific locations in a cell was reserved for the intersection of the boundary with the cell. This allowed to use templates also for the meshing of boundary cells, but a limit on the geometrical complexity of the domain being discretized was imposed. In the more recent implementations [15,68], where the use of real intersections has effectively eliminated the possibility of application of templates to boundary cells, an element removal concept has been used to mesh the boundary cells. Simultaneously, a small-feature elimination process (on the mesh level) has been introduced to remove small elements arising by cutting off small portions of boundary cells by the domain boundary. It should be mentioned that not all strategies use templates to complete the discretization of the tree structure. Some of them use the element removal concept [44] or Delaunay technique in the global or local sense [28,38,46,71,128]. Although the tree-based approach proved to be quite efficient and robust, it suffers from the limited flexibility of the mesh size control due to the virtual predetermination of the element size by the cell size, which is considered to be a significant deficiency in the framework of an adaptive analysis.

Another very attractive method for the generation of triangular and tetrahedral meshes is the advancing front technique [19,25,26,88,99,105,125,126,137,149,166], which is filling the empty (not yet meshed) space by introducing one new element at a time. Initially, a front, defined as an interface between the meshed and unmeshed parts of the domain, and consisting of facets (segments in 2D, triangles in 3D) forming the boundary of the domain, is established. The generation then continues on the basis of an element (triangle in 2D, tetrahedron in 3D) removal algorithm that consists in the selection of an appropriate facet from the front and an attempt to complete it to a proper element, either by picking a suitable existing node in the local neighbourhood of the facet or by creation of a new node. The original facet together with the already existing ones that were used to form the new element are removed from the front and the newly created facets are added to the front. This process is repeated until there are no remaining facets in the front. The algorithm respects very well the desired element size because of its high-quality runtime point placement strategy, which is in contrast with its early version [11] using nodes distributed in a global process prior to the discretization. Since this method is based on a hierarchical strategy (new elements are created from mesh entities of a lower dimension), the integrity of the boundary is guaranteed. The main difficulty in this approach appears to be the need to track the front as it develops and folds over itself, which requires an intensive computational effort to avoid undesirable intersections. Therefore the spatial localization, intersection checking, and front management are considered to be general bottlenecks of the advancing front technique, increasing its computational complexity. Various data structures and strategies [24,55,148] have been proposed to alleviate this negative aspect. It should be mentioned that the advancing front technique can be successfully employed also for the generation of anisotropic meshes (meshes with prestretched and preoriented elements) by simply using the metric given by the transformation between the isotropic and anisotropic spaces. However, extremely high stretches of the order of one thousand (and more), typically used in computational fluid dynamics (CFD), cause numerical difficulties due to the truncation and rounding errors. The advancing front technique is also suitable for the discretization of 3D curved surfaces, either in the real space [126,137] taking into account the surface curvature, or in a parametric space [105,125,166] (if available) taking into account the stretching and distortion induced by the mapping between the real and parametric spaces.

The most widely used approach for the generation of triangular and tetrahedral meshes is based on Delaunay triangulation [1,10,28,30,31,34,38,42,46,67,71,75,93,94,103,113,124,127,128,140,143,145,147], which is a method modifying an existing grid by introducing new nodes. In this technique, discretization of the domain boundary and an initial mesh (possibly consisting of just one element), satisfying the Delaunay property and surrounding the whole domain, have to be available. Each node of the boundary discretization is then inserted into the mesh. Since the node insertion typically violates the Delaunay property, a local mesh transformation has to be carried out after each node insertion to recover it. This transformation is usually performed on the basis of Bowyer's [3] or Watson's [4] algorithm using the ideas of Voronoi diagram. All elements the circumscribed sphere of which contains the node being inserted are removed from the dicretization and the created cavity is remeshed taking into account the inserted node. Further nodes are then inserted into the interior of the domain. These nodes are either created during the runtime [93,113,145] or gathered from a point distribution generated before the actual meshing [1,10,28,38,46,71,128,147]. The former approach provides usually a much better node distribution with respect to the desired mesh size variation. Again, the local transformation process must be accomplished after each internal node insertion. Note that the insertion of internal nodes can be also governed by alternative (usually combined) approaches not necessarily yielding a Delaunay triangulation [124,140]. Finally, the original boundary has to be recovered (by overriding the Delaunay property) to ensure compatibility between the model and the mesh, which is a non-trivial task, especially in 3D. This computationally intensive phase may be skipped if a constrained version of Delaunay triangulation is adopted in which the connections near the boundary are restricted to ensure that the boundary discretization remains intact all the time [30]. The popularity of Delaunay triangulations lies especially in their solid mathematical background, merely referring to a particular connectivity associated with a given set of points, and also in the fact that a valid mesh is maintained throughout the grid generation process. Thus the grid represented by a simple data structure serves as an efficient tool for searching and spatial localization. It should be noted that the quadtree and octree based techniques may represent also one possible realization of Delaunay triangulations under the assumption that templates satisfy the Delaunay property [28,38,46]. Similarly to the advancing front technique, Delaunay triangulation can be also employed for generation of anisotropic meshes using the appropriate metric reflecting the desired stretching and orientation of elements [42]. On the other hand, the use of Delaunay technique for triangulation of 3D curved surfaces is not as straightforward because the formulation of Delaunay property in 2.5D (2D elements in the 3D space) is not available. This is usually resolved by application of Delaunay property to mesh entities from the local neighbourhood of the inserted node projected to the tangential plane at that node [127,143].

*Daniel Rypl
2005-12-07*