A Survey of Unstructured Mesh Generation Technology
2. Tri/Tetrahedral Meshing
TET/TRI
METHODS
Triangle and tetrahedral meshing are by far the most common forms of unstructured mesh generation. Most techniques currently in use can fit into one of three main categories:
Although there is certainly a difference in complexity when moving from 2D to 3D, the algorithms discussed are for the most part applicable for both triangle and tetrahedral mesh generation.2.1 Octree
The Octree technique was primarily developed in the 1980s by Mark Shephard's [5][6] group at Rensselaer. With this method, cubes containing the geometric model are recursively subdivided until the desired resolution is reached. Figure 1 shows the equivalent two-dimensional quadtree decomposition of a model. Irregular cells are then created where cubes intersect the surface, often requiring a significant number of surface intersection calculations. Tetrahedra are generated from both the irregular cells on the boundary and the internal regular cells. The Octree technique does not match a pre-defined surface mesh, as an advancing front or Delaunay mesh might, rather surface facets are formed wherever the internal octree structure intersects the boundary. The resulting mesh also will change as the orientation of the cubes in the octree structure is changed and can also require. To ensure element sizes do not change too dramatically, a maximum difference in octree subdivision level between adjacent cubes can be limited to one. Smoothing and cleanup operations can also be employed to improve element shapes.
Figure 1. Quadtree decomposition of a simple 2D object
From the survey, only four of the 38 codes generating tetrahedral
meshes reported using some form of octree technique. SCOREC
[7] at Rensselaer
develops a set of mesh generation tools called MEGA that utilizes the
Octree technique that is available through their partners program. A
public domain octree mesh generator called QMG
[8] is available from Steve
Vivasis at Cornell.
2.2 Delaunay
By far the most popular of the triangle and tetrahedral meshing
techniques are those utilizing the Delaunay
[9] criterion. The
Delaunay criterion, sometimes called the "empty sphere" property
simply stated, says that any node must not be contained within
the circumsphere of any tetrahedra within the mesh. A circumsphere
can be defined as the sphere passing through all four vertices
of a tetrahedron. Figure 2 is a simple two-dimensional illustration
of the criterion. Since the circumcircles of the triangles in (a)
do not contain the other triangle's nodes, the empty circle property
is maintained. Although the Delaunay criterion has been known for
many years, it was not until the work of Charles Lawson
[10] and Dave
Watson[11]
that the criterion was utilized for developing algorithms
to triangulate a set of vertices. A simple public domain 3D Delaunay
triangulation program called Qhull is available from the University
of Minneapolis. The criterion was later used in developing meshing
algorithms by Timothy Baker[12]
at Princeton, Nigel Weatherill[13]
at Swansea, Paul-Louis George[14]
at INRIA among others.
Figure 2. Example of Delaunay criterion. (a) maintains the criterion while (b) does not.
The Delaunay criterion in itself, is not an algorithm for
generating a mesh. It merely provides the criteria for which
to connect a set of existing points in space. As such it is
necessary to provide a method for generating node locations
within the geometry. A typical approach is to first mesh the
boundary of the geometry to provide an initial set of nodes.
The boundary nodes are then triangulated according to the
Delaunay criterion. Nodes are then inserted incrementally
into the existing mesh, redefining the triangles or tetrahedra
locally as each new node is inserted to maintain the Delaunay
criterion. It is the method that is chosen for defining where
to locate the interior nodes that distinguishes one Delaunay
algorithm from another.
2.2.1 Point insertion
The simplest point insertion approach is to define nodes from
a regular grid of points covering the domain at a specified
nodal density. In order to provide for varying element sizes,
a user specified sizing function can also be defined and
nodes inserted until the underlying sizing function is
satisfied. Another approach is for nodes to be recursively
inserted at triangle or tetrahedral centroids. Weatherill and Hassan
[13]
propose a tetrahedral mesh generation scheme where nodes are inserted
at a tetrahedron's centroid provided the underlying sizing function
is not violated.
An alternate approach is to define new nodes at element circumcircle/sphere centers as proposed by Chew[15] and Ruppert[16]. When a specific order of insertion is followed, this technique is often referred to as "Guaranteed Quality" as triangles can be generated with a minimum bound on any angle in the mesh. Jonathon Shewchuk[17] at CMU has developed a 2D version of this algorithm and makes it available free of charge for research purposes.
Similar to the circumcircle point insertion method, another technique introduced by Rebay[18] is the so-called, Voronoi-segment point insertion method. A Voronoi segment can be defined as the line segment between the circumcircle centers of two adjacent triangles or tetrahedra. The new node is introduced at a point along the Voronoi segment in order to satisfy the best local size criteria. This method tends to generate very structured looking meshes with six triangles at every internal node.
Another method, introduced by Marcum[19] is an advancing front approach to node insertion. Nodes are inserted incrementally, but added from the boundary towards the interior. Each facet is examined to determine the ideal location for a new fourth node on the interior of the existing Delaunay mesh. The node is then inserted and local reconnection is performed. This method tends to generate elements well aligned with the boundary with a very structured appearance to the mesh. Dave Marcum provides both a 2D and 3D version of his mesh generators through the ERC [20] at Mississippi State.
One straightforward method used by INRIA[21]
in their mesh generator
GSH3D[22], is point
insertion along edges. A set of candidate
vertices is generated by marching along the existing internal
edges of the triangulation at a given spacing ratio. Nodes
are then inserted incrementally, discarding nodes that would
be too close to an existing neighbor. This process is continued
recursively until a background sizing function is satisfied.
A variety of other methods for point insertion have also been
proposed, but most have a similar flavor to those discussed above.
2.2.2 Boundary Constrained Triangulation
In many finite element applications, there is a requirement that
an existing surface triangulation be maintained. In most Delaunay
approaches, before internal nodes are generated, a three dimensional
tessellation of the nodes on the geometry surface is produced.
In this process, there is no guarantee that the surface
triangulation will be satisfied. In many implementations,
the approach is to tessellate the boundary nodes using a standard
Delaunay algorithm without regard for the surface facets. A
second step is then employed to force or recover the surface
triangulation. Of course by doing so, the triangulation may
no longer be strictly "Delaunay", hence the term "Boundary
Constrained Delaunay Triangulation".
In two dimensions the edge recovery is relatively straightforward. George[23] describes how the edges of a triangulation may be recovered by iteratively swapping triangle edges. The process is considerably more complex in three dimensions, since after recovering all edges in the surface triangulation, there is no guarantee that the surface facets themselves will be recovered. Additional facet recovery operations can be required to maintain the surface triangulation. While the two dimensional recovery process is guaranteed to produce a boundary conforming triangulation, there are cases[23] in three dimensions where a valid triangulation can not be defined without first inserting additional vertices. This fact increases the complexity of any three dimensional boundary recovery procedure. Two different methods presented in the literature for recovery of the boundary include George[14]and Weatherill [13]
In the first approach defined by George[14] and implemented in INRIA's GSH3D[22]software, edges are recovered by performing a series of tetrahedral transformations by swapping two adjacent tetrahedra for three, as shown in Figure 3. Where a swap cannot resolve the edge, nodes must sometimes be inserted. After edges have been recovered, in order to recover the face, additional transformations are performed, mostly characterized by swapping three adjacent tetrahedra at an edge for two. More complex transformations or additional nodes can be inserted during the face recovery phase if the transformations do not resolve the surface facet.
Figure 3. Tetrahedral transformation where two tets are swapped for three.
The second approach defined by Weatherill [13] also involves an edge recovery phase and a face recovery phase. The main difference with this approach is that rather than attempting to transform the tetrahedra to recover edges and faces, nodes are inserted directly into the triangulation wherever the surface edge or facet cuts non-conforming tetrahedra. This process temporarily adds additional nodes to the surface. Once the surface facets have been recovered, additional nodes that were inserted to facilitate the boundary recovery are deleted and the resulting local void retriangulated.
Another approach presented by Barry Joe[25],
is able to avoid the
boundary recovery problem altogether. Provided the geometry is convex,
Joe is able to define a boundary conforming tetrahedral mesh.
The emphasis in this method, rather than attempting to repair the
boundary of an arbitrary non-convex surface triangulation, is to
decompose the geometry into convex regions that can be separately
processed. An older unsupported public domain version of Barry Joe's
code, Geompack, is available from the University of Alberta
[26] via
anonymous ftp.
2.3 Advancing Front
Another very popular family of triangle and tetrahedral mesh
generation algorithms is the advancing front, or moving front
method. Two of the main contributors to this method are
Rainald Lohner
[27][28]
at George Mason University and S. H. Lo
[29][30]
at the University of Hong Kong. In this method, the tetrahedra
are built progressively inward from the triangulated surface.
An active front is maintained where new tetrahedra are formed.
Figure 4 is a simple two-dimensional example of the advancing
front, where triangles have been formed at the boundary. As the
algorithm progresses, the front will advance to fill the remainder
of the area with triangles. In three-dimensions, for each triangular
facet on the front, an ideal location for a new fourth node is
computed. Also determined are any existing nodes on the front
that may form a well-shaped tetrahedron with the facet. The
algorithm selects either the new fourth node or an existing node
to form the new tetrahedron based on which will form the best
tetrahedron. Also required are intersection checks to ensure
that tetrahedron do not overlap as opposing fronts advance
towards each other. A sizing function can also be defined in
this method to control element sizes. Lohner
[28] proposed using
a course Delaunay mesh of selected boundary nodes over which the
sizing function could be quickly interpolated. A version of S.
H. Lo's advancing front mesh generator is available with the ANSYS
[31]
suite of mesh generation tools.
Figure 4. Example of advancing front where one layer of triangles has been placed
A form of the advancing front method, sometimes called "advancing layers", is also used for generating boundary layers for CFD, Navier-Stokes applications. This method lends itself well to control of element sizes near the boundary. Pirzadeh[32] presents a method where the elements are stretched in the direction of the boundary, the expected direction of fluid flow. A public domain version of Pirzadeh's code, VGRID [33] is available from NASA, Langley.
Back to Introduction.
Continue to Hex/Quad Methods.