Snow Globe: An Advancing-Front 3D Delaunay Mesh Refinement Algorithm ††thanks: The work of the author was supported in part by the NIH/NIGMS Center for Integrative Biomedical Computing grant 2P41 RR0112553-12 and a grant from ExxonMobil. The author would also like to thank Christine Pickett, an editor at the University of Utah, for proofreading and finding numerous typos in the paper.
One of the objectives of a Delaunay mesh refinement algorithm is to produce meshes with tetrahedral elements having a bounded aspect ratio, which is the ratio between the radius of the circumscribing and inscribing spheres. The refinement is carried out by inserting additional Steiner vertices inside the circumsphere of a poor-quality tetrahedron (to remove the tetrahedron) at a sufficient distance from existing vertices to guarantee the termination and size optimality of the algorithm. This technique eliminates tetrahedra whose ratio of the radius of the circumscribing sphere and the shortest side, the radius-edge ratio, is large. Slivers, almost-planar tetrahedra, which have a small radius-edge ratio, but a large aspect ratio, are avoided by small random perturbations of the Steiner vertices to improve the aspect ratio. Additionally, geometric constraints, called “petals”, have been shown to produce smaller high-quality meshes in 2D Delaunay refinement algorithms. In this paper, we develop a deterministic nondifferentiable optimization routine to place the Steiner vertex inside geometrical constraints that we call “snow globes” for 3D Delaunay refinement. We explore why the geometrical constraints and an ordering on processing of poor-quality tetrahedra result in smaller meshes. The stricter analysis provides an improved constant associated with the size optimality of the generated meshes. Aided by the analysis, we present a modified algorithm to handle boundary encroachment. The final algorithm behaves like an advancing-front algorithms that are commonly used for quadrilateral and hexahedral meshing.
Delaunay refinement is usually carried out by eliminating one poor-quality tetrahedron at a time. Jumpani and Üngör  extended the work to eliminate tetrahedra associated with a single short edge in the mesh by simultaneous insertion of multiple vertices around the edge. In this paper, we extend the algorithm to simultaneously eliminate tetrahedra all over the mesh by multiple Steiner vertex insertions. In addition, we replace the exisiting randomized algorithm to avoid slivers, which are almost-planar tetrahedra having small dihedral angles and a large aspect ratio, with a combinatorial algorithm for the optimization of nondifferentiable objective functions.
Delaunay refinement is carried out by placing additional Steiner vertices inside the circumspheres of poor-quality elements to obtain a higher-quality mesh [1, 2, 3, 10, 13, 14, 6]. The radius-edge ratio, the ratio between the circumscribing circle and the shortest edge of a tetrahedra, may be used to measure the quality of an element. If the ratio is large, the element is “split” by adding its circumcenter to the mesh (and re-tetrahedralizing the domain). If the ratio is greater than 1, the new vertex introduces edges that are at least as large as the existing shortest edge in the tetrahedron. Thus, the termination of the algorithm can be proved by the packing argument. To account for the elements near the boundary (where constraints may prevent the insertion of new vertices), the ratio has to be greater than . If a Steiner vertex results in a large sliver, the vertex is added to the mesh because a subsequent insertion of a vertex at the circumcenter of the sliver eliminates the sliver and introduces a long enough edge that still guarantees the termination of the algorithm. For small slivers, however, the new vertex is perturbed enough such that it introduces a long-enough edge and its location results in a better-quality tetrahedron. An analysis of the quality of tetrahedra produced in such adversarial cases provides the theoretical bounds on the mesh quality [7, 9, 8].
The perturbation of the vertex location may be formulated as a nondifferentiable optimization problem. In our recent work , we developed a combinatorial algorithm for such problems. We describe this algorithm and some background on mesh generation in Section 2. Some 2D off-center vertex insertion algorithms for Delaunay refinement place vertices to the extent possible from other vertices, but within certain geometric constraints (called “petals”) [17, 18, 4]. These algorithms place only one vertex at a time. We describe these algorithms in Section 3. In our algorithm in Section 4, we avoid constructing small slivers as far as possible by placing vertices using our combinatorial optimization algorithm. We also provide an iterative framework so that multiple vertices can be added simultaneously to further improve the mesh quality. The 2D geometric constraints (petals) get translated to “snow globes” in our 3D algorithm. In Section 5, we carry out geometric analyses to show that our algorithm terminates and produces size-optimal meshes. For single-vertex insertion, the proofs were provided in earlier work, but we carry out a stricter analysis in which we account for the fact that we process poor-quality tetrahedron with shorter edges first. Thus, we obtain an improved constant of proportionality associated with the size optimality of the mesh. This analysis is an extension of Si’s work . This analysis also helped us develop a modified algorithm to handle boundary facets and edges. For multiple-vertex insertion, we prove that small perturbations in vertex positions due to our combinatorial optimization algorithm still satisfy the conditions necessary to prove the termination and size optimality. Possible future research directions are discussed in Section 6.
We present some notations and background on mesh quality metrics, Delaunay refinement, and optimization of nondifferentiable functions.
2.1 Mesh Quality
We begin by defining the following quality metrics and other terms commonly used in mesh generation algorithms.
Radius-Edge Ratio: The radius-edge ratio, , of a simplicial element is the ratio of the radius of its circumsphere and its shortest edge.
Volume-Edge Ratio: The volume-edge ratio, , of a simplicial element is the ratio of its volume and the cube of its shortest edge.
Almost-Good Mesh: A mesh is -almost-good if the radius-edge ratio of all its elements is less than a constant .
Sliver: A sliver is a tetrahedron whose radius-edge ratio is less than a constant , and the volume-edge ratio is less than a constant .
Local Feature Size: The local feature size at a point is the smallest disk/sphere that can be placed at such that at least two nonincident features are enclosed by the disk/sphere.
Given a triangle , the region of potential locations of an off-plane vertex such that it forms a sliver is inside the region shown in Fig. 1, called the “forbidden” region. If is the shortest side of , the forbidden region is bounded by spheres whose radius is and two planes parallel to the plane of at approriate distances such that the volume-edge ratio is below . Li  showed that there are at most a finite number of small111If the circumradius of a sliver is below a certain threshold, it is called a small sliver. slivers that may be formed when a vertex is being inserted near the circumcenter to refine an almost-good mesh.
The size optimality of an algorithm is typically shown using the local feature size, at any point in the domain. The size-optimality of a mesh is proved by showing that the length of each edge has a lower bound proportional to the local feature size. This lower bound shows that the mesh obtained by an algorithm is, at most, a constant times the size of the smallest possible mesh of the desired quality.
2.2 Nondifferentiable Optimization
We provide a brief description of the first-order necessary conditions for the existence of a local optimum of nondifferentiable function at a point. Consider a composite function that is defined as the minimum over a set of differentiable functions at any point , . Clearly, is a nondifferentiable function, and we wish to find that maximizes . A vector is defined as a subgradient of at if an neighborhood such that . If only one function , defines the value of at , the gradient of at is also the subgradient of at . If multiple functions222They are also called the active set functions. define the value of at a point, any convex combination of the gradients of the active set functions is also a subgradient of at . The set of all subgradients is called the subdifferential, , of at .
For unconstrained optimization, a necessary condition for a local minimum is , i.e., some convex combinations of the gradients of the active set functions should vanish. For constrained optimization, the first-order necessary condition includes the gradients of the active constraints in the convex combination. Formally, let be a set of constraints, , that the location has to satisfy. At a local minimum, the Kahn-Karush-Tucker (KKT) conditions have to be satisfied, i.e., , , , and .
For a -dimensional function, if subgradients have a vanishing convex combination, there exists a set of subgradients (among the subgradients) that also have a vanishing convex combination. Thus, it is sufficient to consider all possible combinations of or fewer subgradients to verify if a convex combination vanishes333If the origin is inside the convex hull of the vectors that represent the subgradients, a convex combination vanishes.. In our algorithm, we consider all possible combinations of or fewer constituent functions, , and constraints, . For each combination, we find locations where the constituent functions are equal and intersect all the constraints. The locations are, typically, determined by solving a system of polynomial equations (or nonlinear equations) with variables ( dimensions and the value of the constituent function), which has a finite number of solutions. We consider all solutions for each combination and return the optimal solution. See  for an example involving the maximization of the minimum angle.
3 Related Work
In this section, we highlight prior research in Delaunay refinement.
3.1 Delaunay Refinement
In order to ensure that the algorithm terminates, the vertex has to be inserted at least a finite distance away from all other existing vertices. If the radius-edge ratio is , the circumcenter is at least a distance from all other vertices of the mesh, where is the length of the shortest side of the triangle. In fact, any point inside the circle/sphere with the circumcenter as the center and having a radius of (for ) is, at least, a distance from other vertices. This sphere/circle is called a “picking” region. In our algorithm, we choose a parameter such that the radius of the picking region is . This ensure that every new edge is longer than the current shortest edge in the mesh. It was shown in  that inserting any point in the picking region is sufficient to obtain a size-optimal mesh (with the minimum angle greater than for 2D refinement). In 3D, a larger and a smaller results in a bigger sphere where there is more room for the Steiner vertex to avoid creating slivers, but a larger also implies a poorer-quality tetrahedron and a smaller may a bigger mesh with shorter edges. Thus, the radius-edge ratio and the volume-edge ratio is chosen so that the aspect ratio (or the dihedral angle) of the resulting tetrahedra should be maximized. Li  analyzes this trade-off.
The algorithms above do not specify the exact location where the Steiner vertex should be placed so that the resulting mesh has fewer vertices and elements than the naive circumcenter insertion. In 2D, the answer to that question was provided by Üngör and Erten [17, 18, 4]. In the first paper , Üng̈or considers the user-defined minimum angle specification and carries out a targeted placement of vertices, which results in a triangulation whose minimum angle is the user-specified value. The vertex is placed on the perpendicular bisector of the shortest edge of a poor-quality triangle such that the angle subtended by the edge on the point is equal to the required minimum angle. If the point of intersection is too close to other vertices, the algorithm reverts to the circumcenter insertion technique. This algorithm produces meshes with nearly 50% fewer vertices.
In Erten and Üngör’s paper , they place the Steiner vertex within a “petal” (see Fig. 2), which is the locus of the points at which the shortest edge subtends the required minimum angle, at a location that maximizes the minimum distance from other vertices. Since this choice of the vertex location is within the constraints of the previous algorithm , it can be shown that the technique terminates and produces size optimal meshes with a minimum angle of . In Section 5, we show why placing vertices inside the petals associated with the shortest edges first results in meshes having few vertices and elements. This algorithm produces meshes with nearly 50% fewer vertices than the previous algorithm.
The algorithm in the first paper  was extended to 3D  by placing the vertices on the perpendicularly bisecting plane of the shortest side. Our algorithm extends the idea of the second algorithm to 3D by defining a “snow globe”, which is analogous to a petal, in which we place Steiner vertices as far away as possible from existing vertices while avoiding the forbidden regions defined in Section 2.
3.2 Simultaneous Multiple-Vertex Insertion
In our recent research , we extended Erten and Üngör’s algorithm  such that multiple vertices may be simultaneously placed by the algorithm. Their insertion algorithm maximizes the minimum distance to existing vertices in a greedy manner, i.e., each vertex is placed at such a location by considering existing vertices in the mesh, but not the vertices that may be placed later. In our extension, we consider simultaneously placing vertices inside some of the current petals. We begin with the greedy placement strategy as in the previous algorithm and iteratively maximize the minimum distance. Fig. 2 provides the motivation for this algorithm by providing a comparison of the greedy strategy with our strategy to place the vertices. It can be seen that our strategy provides lexicographically better angles than the greedy strategy.
Given a 2D input vertex set, one of the Voronoi vertices or one of the points of intersection of the Voronoi edges and the convex hull maximizes the minimum distance to any vertex of the input set if the location is constrained to be within the convex hull of the input vertex set.
In order to find the location that maximizes the minimum distance, consider the gradients of the distance function from each of the input vertices. The function is nondifferentialble at the Voronoi edges (and vertices) as distance functions from multiple vertices interact at these locations. As described in Section 2, a local maximum occurs when one, two, or three gradients have a vanishing convex combination. Since the distance function for a single vertex does not have a zero gradient, we consider the other two cases (see Fig. 3). When two gradients have a vanishing convex combination, they are antiparallel, which happens at the midpoint of two vertices. This point is a saddle point if it is inside the convex hull, but it is a local maxima if it is on the edge of the convex hull. It is easy to show that three gradients have a vanishing convex combination at all the Voronoi vertices. Since we have considered all possible local mamima, one of these locations maximizes the minimum distance. ∎
In Erten and Üngör’s algorithm, additional potential locations of a local maximum are considered on the boundary of a petal as KKT conditions may be satisfied at these locations. In our extension , we iteratively maximize the minimum distance for each new Steiner vertex (and constrain them to be within their respective petals) by allowing only one vertex to move in each step and maximizing its distance from all other vertices. We retriangulate the domain when we are satisfied with the vertex distribution. The iterative algorithm is similar to Lloyd’s algorithm to obtain centroidal Voronoi tessellation. The algorithm converges to a stationary point because we always maximize the minimum distance in each step. This algorithm can be extended to 3D in a straightforward manner.
4 Deterministic Algorithms for Sliver-Free Delaunay Refinement
The Snow Globe
For 2D refinement, petals were defined on the shortest side of a poor-quality triangle. For 3D refinement, we define the “smallest” facet as one of the two triangular faces, , that contains the shortest side, , of a poor-quality tetrahedron and the next shortest side that is adjacent to 444This reason for this choice becomes clear is Section 5.. If the circumradius of is , its radius-edge ratio is . If , the radius-edge ratio threshold for an almost-good mesh, a snow globe is defined by the sphere of radius circumscribing the facet as in Fig. 4. Note that there are two such spheres, and we consider the sphere with the larger part on the same side of as the fourth vertex of the tetrahedron. If the radius-edge ratio of all facets of a tetrahedron are too large, a snow globe may not be defined. In Section 5, we will elaborate of what needs to be done under such circumstances.
Additional Geometric Constraints
We know that the inserted Steiner point has to be at a distance from the vertices of facet . Thus, a vertex is constrained to be within the picking region described in Section 3. In addition, the forbidden regions shown in Fig. 1 are also added as constraints to avoid each potential small sliver.
Objective Function Formulation
For the distance maximization-based routine below, the objective function is simply the maximization of the minimum Eucleadian distance between the Steiner vertices and the existing vertices. We solve systems of four four-variable polynomials, where the minimum distance is itself another variable to be computed.
4.2 Single-Vertex Delaunay Refinement Algorithm
The distance maximization-based routine is quite simple. As in , we place our vertices as far as possible from existing vertices within the snow globe and the picking region, but outside the forbidden region. Our combinatorial optimization algorithm may be used to maximize the distance under these constraints. If an inserted vertex encroaches555A vertex encroaches upon a boundary facet or a boundary edge if it is inside the diametral “lemon”  of the facet or diametral “lens”  of the edge. A different definition is provided in  in which a vertex encroaches upon a boundary edge/facet if it is invisible from the interior of the triangle/tetrahedron due to the boundary edge/facet. Either definition may be used. upon a boundary facet or an edge, all vertices with the diametral circle/sphere are deleted and the midpoint/circumcenter of the edge/facet is added to the mesh. The readers are directed to prior research [14, 4, 5, 15] to learn more about handling such cases. It is important to add vertices in the snow globes associated with shortest edges first to get small meshes. In the next section, we provide a technique to enforce this even for boundary-encroaching Steiner vertices and boundary elements. A disadvantage of using the distance-based technique is that the definition of the sliver is essential to avoid creating the slivers; the definition is explicitly used as a constraint. The angle maximization-based routine provided in the appendix does not have this disadvantage, but it may produce larger, better-quality meshes.
4.3 Multiple-Vertex Delaunay Refinement Algorithm
For multiple vertices, we use an iterative technique for the distance maximization routine. We consider several snow globes and maximize the minimum distance one by one iteratively. After each such vertex movement, the forbidden regions are updated. We want to ensure that poor-quality tetrahedra with shortest sides are processed first. Therefore, we follow an incremental vertex placement strategy, where we place an addition vertex and redistribute them until we run out of the current set of snow globes or the minimum distance is less than the shortest edge of the tetrahedra responsible for the current or next snow globe. The algorithm is presented Algorithm 1. Note that the perimeter-based distance mentioned in the algorithm is the not just the picking region constraint, but also includes the distance from “candidate” Steiner vertices that are to be added in this step.
5 Proofs and Guarantees
The analysis of the algorithm in this section provides the reasoning for our choice of the ordering of poor-quality tetrahedra (and the corresponding snow globes) based on the shortest sides. The analysis also paves the way for handling vertex encroachment on the boundary facets and edges in a slightly different manner. The main difference between our analysis and that of Si  is that we consider the ordering and also constrain the new vertex to be within a certain distance dictated by the dimensions of the snow globes. This feature makes this technique an advancing front technique. Important: This analysis assumes that the length of the shortest edge adjacent to vertex , denoted by , is lsf or less for all vertices in the initial piecewise linear complex (PLC). We discuss the extension to a general case in the appendix. We also provide some additional analysis and discussion about our algorithm in the appendix.
For the single-vertex insertion algorithms, the proof of termination is identical to Foteinos et al.’s  work because we place the vertices within the geometric constraints defined by their algorithm.
The vertex insertion algorithm terminates for appropriate and values that are prescribed in the analysis by Li .
For , the maximization of the minimum distance inside the snow globe yields an edge length of at least the shortest side of the poor-quality tetrahedra. Since only a finite number of small slivers are possible, there exists some for a point inside the picking region and the snow globe, but outside the forbidden region (corresponding to the and ). Note that since we consider not only the picking region but also the intersection of the picking region with the snow globe, our bounds may be slightly different. We do not, however, analyze the bounds because it does not provide any additional insight into the problem, and the theoretical bounds are too small when compared with prior practical results. ∎
The multiple-vertex insertion algorithm terminates.
Every vertex movement in the distance-based algorithm increases the minimum distance between the vertex of interest and some other vertex. Thus, the minimum distance can only increase at every step. ∎
5.2 Size Optimality
We will prove the size optimality for the 2D single-vertex insertion algorithm. Based on the analysis, we make suitable modifications to construct a 3D size-optimal algorithm. Then, the size optimality of the 3D algorithm is self-evident. We assume that the desired radius-edge ratio is , and the desired volume-edge ratio is . The idea behind the analysis is to show that when we order the insertion for petals associated with smaller edges first, the vertices are inserted in an advancing front manner. Thus, new vertices are inserted relatively close to the existing set of good-quality triangles. The new poor-quality triangles have longer edges, and they will insert vertices slightly further away. In this process, the lengths of the new shortest edges of poor-quality elements progressively increase. This progressive increase is absent when poor-quality triangles are ordered in a different manner or when the vertices are inserted away from the petals. We will show this result assuming that boundary encroachments do not play a role. Then, we make suitable modifications to our algorithm when boundary edges/facets come into play.
There are constants and such that , where is the length of the shortest side of the poor-quality triangle in consideration and is the length of the shortest new edge adjacent to the current shortest side after the insertion of the vertex in the petal.
Since we enforce the picking region constraint, . The size of the petal is a function of the shortest side of the mesh. Thus, the Steiner vertex may be inserted within a certain distance from the vertices on the shortest edge. We will call as in the rest of the paper. If a snow globe exists, this lemma extends to 3D as well. ∎
As mentioned in the beginning of this section, we assume that is also the minimum local sizing field in the domain. For the sake of convenience, it is easy to view the vertex insertion procedure as taking place in stages. In the first stage, we insert vertices with the petals associated with poor-quality triangles of the shortest sides with lengths from to , where is the shortest edge in the mesh. In the second stage, the lengths of the shortest edges are from to (due to the lemma above). In the stage, the shortest sides have lengths from to .
At each stage of the algorithm, the locations of the inserted Steiner vertices have constraints on the local sizing field at their respective locations.
The local sizing field, lfs is a Lipschitz function, i.e., . Thus, in the first state, if a Steiner vertex is inserted at point , . In the stage, . ∎
We define a protected region around each vertex as follows. It is easy to see that no new Steiner vertex may be added inside the protected region.
Protected Region: A protected region is a circle/sphere centered at every vertex in the mesh having the radius of the shortest side of a poor-quality tetrahedron among all such tetrahedra at any stage of refinement.
The protected region grows exponentially after each stage.
The length of the shortest edge of a poor-quality tetrahedron in the first stage is . In the second stage, the length becomes . In the , the length becomes . Thus, the protected region grows exponentially. ∎
The single-vertex insertion algorithm produces size-optimal meshes.
In order to prove this statement, we will use Si’s  technique. Consider a vertex and the set of petals into which it was inserted666A vertex may be part of multiple petals by chance.. Let us define the parent of as the closest vertex among the vertices to which the petals belong. Let the edge associated with the petal be of length . We know that , where . Let the parent of be and the parent of be and so on. Let the corresponding lengths of the shortest sides be , , and so on. Let the corresponding ratio of the lengths be . If such a sequence ends at , where is the input vertex, the local sizing field at is bounded from above by . The length is equal to by definition. The ratio of the local sizing field to the length is given by
Since the ratio is bounded and the protected region grows, the vertex closest to can only be placed at a distance of at least . Thus, the ratio of the local sizing field at a vertex to the shortest edge edjecent to the vertex is given by
, and we obtain a size-optimal mesh. ∎
The constant of proportionality associated with our algorithm is better than the one derived by Si  as we have considered the ordering of the poor-quality elements. Note that the analysis was carried out for 2D refinement without considering boundary encroachment. Below, we explain the necessary modifications for the 3D refinement algorithm and for handling boundary encroachment.
5.2.1 3D Refinement
For the 3D Delaunay refinement algorithm, when the snow globe is present for a triangular facet, the edge length-based conditions above hold. When the snow globe is not present, we have to enforce the conditions. Consider a poor-quality tetrahedron that does not have a snow globe. Consider its facet with the shortest edge and the next shortest edge connected to it. Construct a petal associated with the triangle for some radius-edge ratio , and then construct the spindle torus costructed by rotating the petal about the shortest edge. Add a vertex inside the spindle torus and the parameter-based picking region, retetrahedralize the domain, and construct the snow globe on the new facet if necessary.
The intersection of the spindle torus and the picking region is non empty.
Consider the plane formed by the shortest edge of the tetrahedron and its circumcenter. The projection of the spindle torus on the plane is a petal, and the projections of the circumcircle and picking region on the plane are two concentric circle. For and , it is easy to show that the intersection of the petal and picking circle is nonempty. Thus, the intersection of the spindle torus and the picking region is nonempty. ∎
5.2.2 Handling Boundary Encroachment
Consider the queue that dictates the order of the petals/snow globe into which a Steiner vertex is inserted based on the shortest side associated with each petal/snow globe. If a potential Steiner vertex in a petal/snow globe encroaches upon a boundary edge or facet, the petal should be ordered in the queue as if its shortest side is either (a) its own shortest side or (b) half the length of the boundary edge or the radius of the circumcircle of boundary facet divided by , i.e., , whichever is shorter. Let us call this the effective length, . The idea is to get rid of poor-quality triangles with the shortest edges first. Since such petals or snow globes can bring about poor-quality triangles with shorter edges, the petals are appropriately ordered.
If a Steiner vertex does encroach upon the boundary edge/facet, and if there exists a point inside the petal/snow globe and the picking region, place the vertex at such a point. If not, first, we choose one of the vertices of the encroached boundary edge/facet that has the lowest local sizing field. From that vertex , we choose a point on the encroached boundary edge/facet such that , where and (see Fig. 5). If we are handling a facet, we choose on the line joining and the circumcenter of the facet. Second, add to the mesh, and delete all vertices within the circle/sphere formed with as the center and as the radius. Finally, retetrahedralize the domain, which creates new triangles/tetrahedra having the shortest side of length, at least, . Moreover, the vertex is added close to the “advancing front”. Thus, our size optimality results still hold when boundary edges/facets are considered in the above lemmas and theorem. Notice that when , the algorithm is equivalent to the prior Delaunay refinement algorithms.
5.2.3 Size Optimality for the Multiple-Vertex Insertion Algorithm
The multiple-vertex insertion algorithm terminates with size-optimal meshes.
In Algorithm 1, we simultaneously add multiple vertices only if the minimum distance is greater than times the shortest edge associated with the next snow globe. Thus, the constraints on the minimum and the maximum distance from the current “front” are still maintained. We enforce the ordering on the snow globe explicitly through the suitable condition in the while loop in Algorithm 1. Therefore, all of the lemmas above hold, and a size-optimal mesh is generated. ∎
6 Future Work
We have extended a 2D Delaunay refinement algorithm to 3D (petals to snow globes; deterministic maximization; multiple vertex insertion) and consolidated 3D algorithms (sliver-free refinement; handling boundary encroachment) to develop our advancing front algorithm. We have also extended the analysis by Si  to explain why our algorithm is likely to provide smaller high-quality meshes. We will extend this algorithm for constrained Delaunay refinement and implement it to study its performance in practice. Other future directions for research include adapting the algorithm for surface meshing algorithms and anisotropic meshing in 2D.
-  L. P. Chew. Constrained delaunay triangulations. In Proceedings of the Third Annual Symposium on Computational Geometry, SCG ’87, pages 215–222, New York, NY, USA, 1987. ACM.
-  L. P. Chew. Guaranteed-quality triangular meshes. Technical report, Department of Computer Science, Cornell University, 1989.
-  L. P. Chew. Guaranteed-quality mesh generation for curved surfaces. In Proceedings of the Ninth Annual Symposium on Computational Geometry, SCG ’93, pages 274–280, New York, NY, USA, 1993. ACM.
-  H. Erten and A. Üngör. Quality triangulations with locally optimal steiner points. SIAM Journal on Scientific Computing, 31(3):2103–2130, 2009.
-  P. A. Foteinos, A. N. Chernikov, and N. Chrisochoides. Fully generalized two-dimensional constrained delaunay mesh refinement: Revisited. SIAM Journal Scientific Computing, 32(5):2659–2686, 2010.
-  R. Jampani and A. Üngör. Construction of sparse well-spaced point sets for quality tetrahedralizations. In Proceedings of the 16th International Meshing Roundtable, October 14-17, 2007, Seattle, Washington, USA, Proceedings, pages 63–80, 2007.
-  X.-Y. Li. Sliver-free three dimensional delaunay mesh generation. PhD thesis, University of Illinois at Urbana-Champaign, 2000.
-  X.-Y. Li. Generating well-shaped d-dimensional delaunay meshes. Theoretical Computer Science, 296(1):145 – 165, 2003. Computing and Combinatorics.
-  X.-Y. Li and S. H. Teng. Sliver-free three dimensional delaunay mesh generation. In Symposium on Discrete Algorithms (SODA), pages 28–37. SIAM, 2001.
-  J. Ruppert. A delaunay refinement algorithm for quality 2-dimensional mesh generation. Journal of Algorithms, 18(3):548 – 585, 1995.
-  S. P. Sastry. Maximizing the minimum angle with the insertion of steiner vertices. In Canadian Conference on Computational Geometry. 2015. accepted. Reviewers may download the current draft at http://www.sci.utah.edu/~sastry/sastry_cccg.pdf for reference.
-  S. P. Sastry. Optimal insertion of multiple steiner vertices for two-dimensional delaunay mesh refinement. In International Meshing Roundtable. 2015. submitted. Reviewers may download the current draft at http://www.sci.utah.edu/~sastry/sastry_imr.pdf for reference.
-  J. R. Shewchuk. Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In M. C. Lin and D. Manocha, editors, Applied Computational Geometry: Towards Geometric Engineering, volume 1148 of Lecture Notes in Computer Science, pages 203–222. Springer-Verlag, 1996. From the First ACM Workshop on Applied Computational Geometry.
-  J. R. Shewchuk. Delaunay Refinement Mesh Generation. PhD thesis, Carnegie Mellon University, 1997.
-  J. R. Shewchuk and H. Si. Higher-quality tetrahedral mesh generation for domains with small angles by constrained delaunay refinement. In Proceedings of the Thirtieth Annual Symposium on Computational Geometry, SOCG’14, pages 290:290–290:299, New York, NY, USA, 2014. ACM.
-  H. Si. An analysis of shewchukâs delaunay refinement algorithm. In B. W. Clark, editor, Proceedings of the 18th International Meshing Roundtable, pages 499–518. Springer Berlin Heidelberg, 2009.
-  A. Üngör. Off-centers: A new type of steiner points for computing size-optimal quality-guaranteed delaunay triangulations. In M. Farach-Colton, editor, LATIN 2004: Theoretical Informatics, volume 2976 of Lecture Notes in Computer Science, pages 152–161. Springer Berlin Heidelberg, 2004.
-  A. Üngör. Off-centers: A new type of Steiner points for computing size-optimal quality-guaranteed delaunay triangulations. Computational Geometry, 42(2):109 – 118, 2009.
Appendix A Extension to a General Class of Piecewise Linear Complexes
In our ten-page abstract above, we assumed that the local sizing field at every input vertex is equal to (or smaller than) the length of the shortest edge in the input PLC. There are PLCs for which the condition does not hold. For instance, consider a line segment and a point close to it. If no vertices other than and are closer to , the local sizing field at is equal to its distance to , which may be smaller than and . We will first analyze why our algorithm may not function as an advancing-front algorithm in such instances and then provide a suitable preprocessing algorithm to achieve our objective of placing vertices in an advancing-front manner.
Let us assume that is a poor-quality triangle. If we were to carry out Steiner vertex placement as in our previous algorithm, would be split in an appropriate order depending on and . Further, if the Steiner vertex for encroaches upon segment , a new vertex would be added on , and all free vertices in the circle/sphere with as the center and having a radius (or ) would be deleted. Since is not a free vertex, it is not deleted even if it is inside the circle. Thus, we would create a triangle with a smaller edge length . This sabotages our idea of introducing progressively longer edges after each Steiner vertex insertion.
In order to tackle the issue above, we develop an algorithm to preprocess the PLC to add auxiliary Steiner points into the PLC such that the diametral/equatorial circles/spheres on the boundary edges and facets do not contain the input vertex. Thus, whenever a new vertex is introduced, the length of the shortest edge of a new triangle/tetrahedron can only be greater than the length of the shortest existing triangle/tetrahedron. Such a preprocessing technique would result in an advancing-front Steiner vertex insertion routine. Since additional vertices are added into the PLC, the constant associated with the size optimality will be compromised.
Let us first consider the 2D problem (see Fig. 6). If an input vertex lies inside the diametral circle on line segment , the perpendicular projection of on the line joining and lies on the line segment joining and . Let the point of projection be . Further, . If the point is added to the PLC, does not lie inside the diametral circle of edge or , but may not be the most optimal location of an auxiliary Steiner point because it may be arbitrarily close to or . Thus, we perturb the vertex to a location such that does not lie inside the new diametral circles, and its location maximizes the minimum distance from all other vertices.
In 2D refinement, may be placed such that the constant associated with the size optimality does not reduce, i.e., and does not encroach upon the line segment .
If encroaches upon the segment , and must lie inside the line segment . Thus, may be placed such that , and we obtain a triangulation where does not encroach on the boundary segment . If we process vertices in the increasing order of their local sizing field, no other input vertex will be closer than the distance to . Thus, in 2D refinement, the constant associated with the size optimality does not reduce. ∎
We recommend the same algorithm for 3D PLCs as well, but we do not yet have a bound on the size optimality in one of the cases described below. As above, we process the vertices in the order of their distance to the next nearest input features, i.e., the local sizing field at the vertex. Consider a point that is encroaching upon the triangular facet . Let the point of projection from onto the plane formed by be , and let us assume that . The new vertex may be added at if it is a sufficient distance from other vertices. Since is right above , no circumsphere of any triangle containing the vertex at can encroach upon it. Further, due to the Delaunay triangulation of the facet, does not lie not inside any of the circumcircles of other triangles. If is arbitrarily close to one of the vertices of the triangle, say, , we perturb to a new location on the line joining and . We will show that can be as large as before lies inside the equatorial spheres of the triangular facets. If is too close to some other vertices, some other point on the line segment may be chosen such that it maximizes the minimum distance from all other vertices.
If , no diametral sphere of triangular facets encroaches upon vertex .
Consider the unique plane perpendicular to the plane of and on the line joining points and . Consider the cross-section of a equatorial (hemi)sphere on the plane. This cross-section is similar to Fig. 6 (c) and (d) if diametral semicircles are added to it. By construction, . If is inside an equatorial sphere, the radius of the sphere is greater than or equal to because the center of the sphere lies on the plane of . The cross-section of the (hemi)sphere forms an arc smaller than or equal to a semicircle with the radius less than . Note that no arc on the plane may enclose vertices or because the triangular facets form the Delaunay triangulation. Thus, if , no equatorial sphere of triangular facets encroaches upon vertex . ∎
We still have not considered the following three cases: (a) when point lies outside all the triangles whose equatorial spheres enclose , (b) vertex lies close to an input line segment, and (c) point or lies very close to an input line segment. Case (a) cannot happen because the vertex would have to be inside an equatorial sphere of only one of the two such intersecting spheres, but the circle of their intersection projects onto the common edge of the two triangular facets. Thus, if the vertex is inside the equatorial sphere of one of the triangles, but its projection is outside, it is also inside the equatorial sphere of the triangle over which its projection falls. For case (b), we may follow the algorithm above used for 2D refinement. We will provide some nonrigorous analysis for case (c) below.
A natural solution for case (c) is to place the vertex at a point on the input line segment furthest away from existing vertices such that does not lie inside the resulting equatorial spheres after the insertion of the new vertex. What we do not yet know is how close to an existing vertex (relative to the local sizing field) such a point can lie. Consider the unique plane that intersects the encroached upon triangular facet on the input line segment and passes through the vertex . Since we know that or is closer to the input line segment than the distance , we know that the dihedral angle between our new plane and the plane of the triangular facet is less than 45 degrees. Thus, the cross-section of the equatorial (hemi)sphere (of a triangular facet incident on the input line segment) on the new plane forms a petal. If is inside the petal, the radius of the circle from which the petal is formed is constrained (but smaller than ). The constant associated with the constraint results in a proportional change in the constant associated with the size optimality of the mesh.
We have not yet considered other equatorial spheres that may enclose after is inserted into the mesh. If is still inside an equatorial sphere, it cannot be the equatorial sphere associated with a triangle that is incident on the input segment considered above (because of our choice of ). From the reasoning for case (a) above, since does not lie inside the equatorial sphere on the triangular facet in which its its projection lies, it does not lie in the equatorial sphere of any other triangular facet on the same plane. The analysis is nonrigorous, and we leave the rigorous analysis as future work.
Appendix B Angle Maximization-Based Vertex Placement
One of the disadvantages of the distance maximization-based routine provided in the paper is that the definition of sliver is necessary to place the vertex. Here, we present a slightly different algorithm where we attempt to improve the dihedral angles between the facets of a tetrahedra by maximizing the weighted distance from the potential new tetrahedra that may be formed as a result of the insertion of the new vertex.
For the angle maximization-based routine, we maximize the minimum weighted distance from the plane of those triangles that may form small slivers. The weights are proportional to the area of the triangles. Since the volume of the potential sliver is proportional to the area of the triangle and the distance of the Steiner vertex from the plane of the triangle, the increase in the weighted distance also increases the dihedral angle as well as the volume-edge ratio. As mentioned in Section 2, there may be, at most, a finite number of such facets. Li  showed that a small perturbation of the Steiner vertex may change the connectivity in the resulting Delaunay tetrahedralization. Thus, it is important to consider all potential triangular facets that may form small slivers.
For single-vertex insertion, however, this technique can place points that are too close to the existing vertices. Hence, the additional geometrical constraint in the form of the picking region is necessary. Our algorithm maximizes the weighted distance within the intersection of the picking region and the snow globe to insert the Steiner vertex. In the case of vertex encroachment, the same modification as above is applied. Although we expect better dihedral angles from this routine, it is likely to produce larger meshes.
For multiple-vertex insertion, instead of the picking region constraint, we place another constraint that does not allow vertices to get too close to each other. Specifically, we do not allow two vertices to be closer than , where is some feasible constant and is the larger of the shortest sides of the two poor-quality tetrahedra that the vertices are attempting to eliminate. This constraint is necessary to guarantee the termination and size optimality of the algorithm.
Appendix C Additional Analysis
The ordering of petals and snow globes is the reason why the Delaunay refinement algorithms provide smaller meshes. The algorithms attempt to get rid of poor-quality elements with short edges as quickly as possible.
Assuming boundary encroachment does not play a role, for petal- and snow globe-based refinement algorithms, poor-quality triangles or tetrahedra associated with the shortest edge in an existing poor-quality element can be eliminated with a finite number of Steiner vertex insertions.
For 2D refinement, with at most two additional Steiner vertices, poor-quality triangles associated with a short edge may be eliminated by placing the vertices on the either side of the shortest edge within the respective petals and the picking regions. For 3D refinement, because we eliminate slivers by avoiding the forbidden regions, at least one new tetrahedron is of good quality with bounded dihedral angles. Thus, after a finite number of Steiner vertex insertions, all poor-quality tetrahedra “around” the short edge are eliminated. ∎
If circumcenter insertion is used, there is no guarantee that the new triangle/tetrahedron is of acceptable quality, and the triangle/tetrahedron may have to be refined again. We prove the following lemma only for 2D refinement. The existence of slivers makes it hard to prove the lemma for 3D refinement.
Assuming boundary encroachment does not play a role, for circumcenter-based refinement algorithms, poor-quality triangles associated with the shortest edge in an existing poor-quality element can be eliminated with a Steiner vertex insertions, where is the diameter of the domain and is the length of the edge.
Let us assume the radius-edge ratio of the poor-quality element is . After the first vertex insertion, the quality of the new element of the shortest edge becomes . After the second vertex insertion, it becomes , and so on until vertices. We know that , . Since we assume boundary encroachment does not play a role, is bounded by . Thus, . ∎
The size of the mesh is , where is the diameter of the domain and is the minimum local sizing field in the domain.
Asymtotically, the size of the mesh does not change due to our algorithm. Thus, the proof is identical to the one provided by Si . The constant associated with the size optimality is different (as show in the ten-page abstract) due to the ordering. ∎