 # A 2D Advancing-Front Delaunay Mesh Refinement Algorithm

I present a generalization of Chew's first algorithm for Delaunay mesh refinement. In his algorithm, Chew splits the line segments of the input planar straight line graph (PSLG) into shorter subsegments whose lengths are nearly identical. The constrained Delaunay triangulation of the subsegments is refined based on the length of the radii of the circumcircles of the triangles. This algorithm produces a uniform mesh, whose minimum angle can be at most π/6. My algorithm generates both truly Delaunay and constrained Delaunay size-optimal meshes. In my algorithm, I split the line segments of the input PSLG such that their lengths are asymptotically proportional to the local feature size (LFS) by solving ordinary differential equations (ODEs) that map points from a closed 1D interval to points on the input line segments in the PSLG. I then refine the Delaunay triangulation (truly or constrained) of the PSLG by inserting off-center Steiner vertices of "skinny" triangles while prioritizing such triangles with shortest edges first. As in Chew's algorithm, I show that the Steiner vertices do not encroach upon any subsegment of the PSLG. The off-center insertion algorithm places Steiner vertices in an advancing front manner such that we obtain a size-optimal Delaunay mesh (truly or constrained) if the desired minimum angle is less than π/6. In addition, even in the presence of a small angle ϕ < π/2 in the PSLG, the bound on the minimum angle "across" the small angle tends to arctan((sinϕ)/(2-cos(ϕ)) as the PSLG is progressively refined. Also, the bound on the maximum angle across any small input angle tends to π/2 + ϕ/2 as the PSLG is progressively refined.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Delaunay mesh refinement techniques are commonly used to generate high-quality meshes in two or higher dimensions. These meshes are usually used to solve partial differential equations (PDEs) with the finite volume method (FVM) or the finite element method (FEM)

. The FVM is typically used to solve fluid-flow problems by defining a control volume surrounding each vertex in a mesh and measuring the flux entering and exiting the control volume. Delaunay meshes are extensively used in the FVM because it is easy to define a control volume around a vertex using a Delaunay mesh and its corresponding Voronoi diagram.

The FEM, on the other hand, may be used with any mesh. For instance, in order to solve isotropic elliptic PDEs, the FEM requires meshes whose elements are “regularly” shaped, i.e., the length of their edges should be nearly identical . In 2D triangular meshes, this regularity requirement translates to triangles having angles as close to 60 degrees as possible. In 2D, given a set of vertices, Delaunay triangulation maximizes the minimum angle among all possible triangulations over the set of vertices. As this property aligns with the goal of mesh generation, Delaunay meshes are also used in the FEM to solve isotropic elliptic PDEs. In this paper, I focus only on 2D meshes.

A challenge in generating Delaunay meshes is to make the meshes conform to both internal and external boundaries, which can be accomplished using either truly Delaunay meshes or constrained Delaunay meshes . A truly Delaunay mesh does not permit any mesh vertex inside the circumcircle of a triangular element. With truly Delaunay meshes, the Delaunay triangulation of mesh vertices recovers the boundaries automatically. In the context of solving PDEs, truly Delaunay meshes are necessary when the FVM is being used, and the mesh might be rather large. For other applications (such as the FEM), a constrained Delaunay meshes might be sufficient.

A constrained Delaunay mesh permits mesh vertices inside the circumcenter of a triangular element if the mesh vertices are not visible (explained in Fig. 1) from any point inside the triangular element. This flexibility allows the generation of smaller high-quality meshes when compared to truly Delaunay meshes. Fig. 1 shows some examples of truly Delaunay triangulation and constrained Delaunay triangulation. Note that it is easy to show that if no vertices are allowed inside the diametral circles of subsegments of an input segment, a Delaunay triangulation of vertices will yield a conforming Delaunay mesh that is truly Delaunay.

A typical Delaunay mesh refinement algorithm starts from an initial Delaunay triangulation of the input geometric domain. The input geometric domain is also called a planar straight line graph (PSLG) since the input is planar and consists of vertices and line segments that may be thought of as an embedded graph. The Delaunay refinement algorithms progressively add vertices to the mesh and retriangulate the vertices such that poorly shaped triangles are eliminated. In Section 2, I will discuss some of the techniques developed to obtain Delaunay meshes. Most algorithms liberally add vertices (also called Steiner vertices) within the input domain. When the techniques, however, attempt to add a Steiner vertex close to the boundaries or outside the domain, it is called an encroachment of the domain. Most algorithms handle this case by splitting the relevant segments in the input PSLG, deleting a few other vertices (only in some algorithms), and retriangulating the remaining vertices. These technique produce well-graded, high-quality meshes.

My technique is a generalization of a technique by Chew . Chew splits the input segments in the PSLG into subsegments whose lengths are nearly identical. The split PSLG segments are triangulated and refined by adding Steiner vertices, which are circumcenters of triangles whose radius of the circumcircle is larger than the length of the shortest subsegment in the split PSLG. The technique never attempts to add a vertex outside the domain. As a result of the nearly uniform splitting of the PSLG segments and the refinement technique, this algorithm produces uniform meshes.

In my technique, I generate well-graded meshes by refining the PSLG such that the lengths of the split segments are asymptotically proportional to the local feature size (formally defined in Section 3) at the end points of the split segments. Such asymptotically proportional splits are achieved by solving an ordinary differential equation, whose solution maps points from a reference line segment to a line segment in the PSLG. I then refine poor-quality triangles by adding their circumcenter or an off-center Steiner vertex [8, 28, 27]. I prioritize poor-quality triangles with shortest edges first. More details are provided in Section 4. In Section 5, I show that my adaptive splitting also ensures that the algorithm never attempts to insert a vertex outside the domain. I have separately analyzed the algorithm for generating truly and constrained Delaunay meshes.

My technique improves the upper bound on the minimum angle in a size-optimal mesh to 30 degrees (except “across” small angles). Moreover, even in the presence of small angles, the technique improves the lower bound on the maximum angle to , where is the smallest angle in the input PSLG. These results hold for both truly and constrained Delaunay meshes, but the constants associated with size optimality (defined in Section 3) are different for the two types of meshes. As expected, truly Delaunay meshes are larger than constrained Delaunay meshes.

## 2 Related Work

An extensive body of literature has focused on generating 2D Delaunay meshes with a specified minimum angle. The first subsection below reviews prior algorithms in which the PSLG does not contain any angle smaller than (or , depending on the algorithm). The second subsection reviews techniques to deal with PSLGs with small angles and the limitations of any algorithm used to mesh such PSLGs.

In most algorithms, a triangle is considered to be of poor quality or “skinny” if its minimum angle is less than some user-defined threshold; it is an input to the algorithm. The ratio of the length of the radius of the circumcircle of a triangle and its shortest length is equal to , where is the minimum angle in the triangle (see Fig. 2). A user may equivalently provide the radius-edge ratio to define a skinny triangle.

### 2.1 Delaunay Mesh Refinement

The first Delaunay mesh refinement algorithm was developed by Chew  to obtain constrained Delaunay meshes. Chew first splits the segments of the input PSLG such that the length of each subsegment is between and , where is small enough that such a division in possible111In his report , Chew provides more details about how to find and how to split the PSLG.. The split PSLG is then triangulated using the constrained Delaunay triangulation algorithm . Then triangles whose radius of the circumcircle is less than are “split” by adding its circumcenter. Upon retriangulation, the addition of the circumcenter eliminates the skinny triangle because the new vertex (the circumcenter) is inside the skinny triangle, which Delaunay triangulation does not allow. Since Chew always adds vertices at a distance of at least from other vertices (Delaunay triangulation ensures that there are no vertices inside the circumcircles of all triangles), we eventually run out of space to add more points. In addition, Chew demonstrated that no triangle will be formed such that the radius of its circumcircle is less than and its circumcircle is outside the domain. Also, since the radius of the circumcircle is less than and the shortest edge length is at least , all the angles in the mesh are greater than or equal to . This algorithm results in a uniform mesh. In this paper, I generalize this algorithm by developing a technique to split PSLG segments in a size-optimal fashion rather than uniformly.

Ruppert [20, 21] developed a similar Delaunay refinement algorithm that splits the PSLG on the fly to generate truly Delaunay meshes. Instead of splitting a triangle based on the length of the radius of the circumcircle, the radius-edge ratio became the criterion. A triangle is split only when the ratio is greater than 1 (it should be greater than for reasons explained below), which corresponds to the minimum angle being . Clearly, if the triangle is split when the ratio is less than , the algorithm does not terminate as it places vertices progressively closer to each other. If the Steiner vertex, however, is outside the domain or inside the diametral circle of a PSLG subsegment, the PSLG subsegment is split at its midpoint. Because the midpoint may introduce short edges, the threshold of the radius-edge ratio to split a triangle has to be increased to , which corresponds to a minimum angle of about 20.7 degrees. Rand  showed that Ruppert’s algorithm (for truly Delaunay meshes) terminates for angles almost as large as 22.2 degrees. Ruppert showed that this algorithm provides a truly Delaunay, size-optimal mesh (elements’ edge lengths are graded based on their proximity to small features in the input PSLG) when the input angles in the PSLG are greater than . This threshold was later lowered to by Shewchuk [22, 24]. Shewchuk was also able to improve the bound on the minimum angle in a triangle to at the cost of size optimality. He split the PSLG segments such that their lengths were in specific ranges, and this restriction in the subsegment length resulted in the loss of the size-optimality guarantee.

Chew  independently devised a second constrained Delaunay refinement algorithm in which he incorporated the idea of refining the PSLG only when necessary. In addition, the technique also removes points within the diametral circle of a subsegment when the subsegment is to be refined. Shewchuk [22, 24] showed that this algorithm produces size-optimal constrained Delaunay meshes when the desired radius edge ratio is greater than , which translates to a minimum angle of about 26.57 degrees.

Shewchuk  improved upon both the techniques of Chew and Ruppert above by refining triangles near any PSLG segment differently from those away from the segments. Shewchuk’s technique ensured that the quality of triangle in the interior is better than the ones near the PSLG boundary segments. Also, he introduced the notion of “diametral lenses” to analyze Chew’s second algorithm and showed the bound of about 26.57 degrees mentioned above.

Miller, Pav, and Walkington [13, 15] showed that in a modified version of Ruppert’s technique (for truly Delaunay triangulation), at least three circumcenters have to be inserted between two refinements of a PSLG segment/subsegment. Thus, in the worst case and when the input is restricted (conditions on the angles between PSLG segments and their lengths), they were able to show that the algorithm terminates with a truly Delaunay, size-optimal mesh if the minimum desired angle is set to some value strictly less than degrees even when the input angle is as small as . Rand  extended this analysis to Chew’s second algorithm (for constrained Delaunay triangulation) and showed that it can produce a size-optimal mesh with a minimum desired angle is set to some value strictly less than such that , which corresponds to about 28.60 degrees.

Üngör and Erten [8, 28, 27]

developed a heuristic technique to place Steiner vertices at “off-center” points such that the shortest edge of skinny triangles subtend the desired minimum angle at the points. When coupled with the prioritization of skinny triangles with shortest edges first, they found that the resulting meshes have fewer elements and vertices. Their technique works for both truly and constrained Delaunay meshes. I use the heuristic technique in my algorithm, and it plays an important role in the proofs.

Foteinos et al. [3, 9] generalized Chew’s second algorithm to show that the insertion of Steiner vertices can happen at any point (not just the circumcenter or the off-center point) in a skinny triangle’s “selection circle”, which is a circle concentric with its circumcircle, but with a radius shorter by the length of the shortest edge of the skinny triangle. In addition, they showed that the PSLG may be split at any point (not just the midpoint) sufficiently far away from the end points. Such an algorithm would still produce a size-optimal constrained Delaunay mesh with quality guarantees. Hudson  generalized the selection circle to a larger selection region.

### 2.2 Managing Small Input Angles

If the PSLG has any small input angles, the results above do not hold for the whole mesh. One may assume that if the algorithms above ignore triangles at small input angles, it should be possible to construct a mesh with larger angles everywhere else. Unfortunately, Shewchuk [22, 24] showed that it is impossible to construct such a mesh. Researcher have since attempted to mitigate the effects of small input angles.

Shewchuk [22, 23, 24] developed a technique he called Terminator that splits skinny triangles near a small input angle only if the new edge that would be introduced in the mesh were long enough to ensure that the algorithm terminates. This technique (for truly Delaunay meshes) does not have a theoretical guarantee of size optimality. He showed that the minimum angle in the mesh may be as large as , where is the smallest angle in the input PSLG. He applied a similar technique for Chew’s algorithm (for constrained Delaunay meshes) and improved the bounds to .

Rand and Walkington  employed a collar-based strategy to protect regions around a PSLG vertex with small angles and prevented vertex insertion near a small angle. Pav and Walkington  employed a similar strategy. This strategy was inspired by Ruppert’s heuristic in his papers [20, 21].

Pav et al. [13, 15] redefined a poor-quality triangle as a skinny triangle whose end points of the shortest edge do not lie on two segments emanating from an input vertex with a small angle. They split all skinny triangles in the mesh that were not “across” a small angle. They ignored skinny triangles across a small input angle, but due to the segment splits beforehand, they showed that the angles in such skinny triangles were at least . This technique was inspired by Ruppert’s concentric shell splitting heuristic. In my algorithm, I split the PSLG before splitting triangles. Thus, I am also able to show almost identical bounds. In addition, the Pav et al. algorithm produces a mesh with a maximum angle of degrees. My algorithm improves this bound.

## 3 Background

The notations, concepts and algorithms presented in this section are used in the advancing-front algorithm.

### 3.1 A “Skinny” Triangle

In this paper, I will use to denote the desired minimum angle in a mesh and to denote the desired radius-edge ratio (see Fig. 2 for an explanation). A skinny triangle is any triangle whose minimum angle is less than or whose radius-edge ratio is greater than .

### 3.2 The Local Feature Size and Size Optimality

The local feature size, denoted as or simply , is the radius of the smallest circle centered at a point that intersects two nonadjacent features of the input PSLG. Note that the local feature size is dependent only on the input but not on the mesh generated. For any point on a line segment in the input PSLG, since is already on an input feature, the feature size at is the minimum of (a) the distance to the nearest feature from such that the feature is not adjacent to and (b) the distance to or , whichever is farthest from . The distance to or is also considered because and are not adjacent features, and both and are inside a disk centered at with radius equal to or , whichever is greater.

Ruppert  introduced this metric and showed that his algorithm produces meshes in which the length of the edges is greater than some fraction of the local feature size at its end points. In other words, the length of any edge in the final mesh produced by his algorithm is greater than , where is some constant. He called the meshes size-optimal meshes. Since then, many researchers have used this metric to show that their algorithms also produce meshes that are size optimal.

The local feature size is a Lipschitz continuous function, i.e., and because the nonadjacent features that are contained in a disk centered at with a radius are also contained in the disk centered at with a radius .

### 3.3 Nonhomogenous Ordinary Differential Equations

My algorithm involves the solution of ordinary differential equations (ODEs) of first and second order. There is an extensive body of literature addressing the solution of such equations analytically and/or numerically. Fortunately, it is possible to derive an analytical solution to the equations presented in Section 4. Below, I will go through the two equations we will see in Section 4.

First, let us consider an equation of the form , where and are constants. Multiplying both sides by , we get . Integrating both sides w.r.t. , we get , where is some constant. Thus, , where is determined using a boundary condition. If , then , where is some constant.

Second, let us consider an equation of the form , where is a constant. The solution to this equation is given by , where and are constants to be determined using two independent boundary conditions.

I use the solution of the differential equations discussed above to adaptively split all line segments of the input PSLG such that their lengths are asymptotically proportional to the local feature size. Given a line segment that is a subsegment of the PSLG line segment, I ensure that the length of is such that , where and are some constants. I will denote the ratio as .

### 3.4 The Off-Center Refinement Algorithm

In earlier Delaunay mesh refinement algorithms, Steiner vertices were added at the circumcenter of skinny triangles in order to produce a mesh with no skinny triangles. Erten and Üngör [27, 8] developed an algorithm in which the position of a Steiner vertex is an off-center point on the perpendicular bisector of the shortest edge of the skinny triangle. The point is chosen such that the angle subtended by the shortest edge at that point is equal to the desired angle specified in the input (see Fig. 3). If the off-center point is farther than the shortest edge than the circumcenter, the algorithm reverts to the circumcenter insertion technique. The algorithm processes poor-quality triangles with the shortest edges first, i.e., it considers the shortest edge in every skinny triangle, and processes the triangle with the shortest edge among those considered edges.

In this algorithm, if the length of the shortest edge of a skinny triangle is , the location of the off-center point is at a distance of from the vertices of the shortest edge, where and is the desired minimum angle in the mesh. Note that , where is the desired radius-edge ratio. Since there are no other vertices within the circumcircle of the skinny triangle (the Delaunay property), this algorithm always places Steiner vertices that are at least a distance of from all other vertices in the mesh and at most a distance of from both vertices of the shortest edge of the skinny triangle. If (and if all the Steiner vertices are placed inside the domain), it is easy to see that the algorithm terminates when it runs out of place to add more vertices in the domain. Also note that users could purposefully insert a Steiner vertex at a distance between and if they choose to slowly grow the size of elements as the vertices are placed away from the PSLG segments.

In my algorithm, I use the off-center Steiner vertex insertion algorithm described above (including the prioritization of shortest edges). I ensure that the adaptive splitting results in Steiner vertex insertion that does not encroach (defined below) upon any of the split segments.

### 3.5 Line Segment Encroachment

All the algorithms mentioned in this paper work without issues only if the Steiner vertices are always placed inside the domain. When a potential vertex is outside the domain, the vertex is considered to have encroached upon the domain. More formally, a Steiner vertex encroaches upon a PSLG subsegment if the skinny triangle that resulted in the Steiner vertex and the Steiner vertex lie on the opposite sides of the PSLG line segment. Note that this definition considers the internal boundaries in the PSLG as well. Since the triangulation obeys the Delaunay property, the circumcircle of a skinny triangle cannot enclose any vertex of a PSLG subsegment. For the circumcenter to encroach upon the subsegment, the circumcircle and the skinny triangle have to be on different sides of the PSLG subsegment. This implies that the skinny triangle has to be inside the diametral circle of the PSLG subsegment (see Fig. 4). In the algorithm, I insert the off-center vertex or the circumcenter, whichever is closer. Therefore, I will always insert a vertex that is at most the distance to the circumcenter. Thus, if the circumcenter does not encroach any PSLG subsegment, we can be guaranteed to insert a vertex that does not encroach upon a PSLG subsegment. If a Steiner vertex encroaches upon a PSLG line segment, prior algorithms impose corrective measures. As in Chew’s algorithm, I will show that no circumcenter of skinny triangles encroaches upon a PSLG line segment if the line segments are appropriately split.

### 3.6 Small Angle

I denote a small angle as . Given a size-optimal splitting of the PSLG such that for any subsegment of length and , I define a small angle as any angle .

### 3.7 Skinny Triangles “Across” a Small Angle

Miller, Pav, and Walkington [13, 15] developed an algorithm that provides guarantees on the mesh quality even in the presence of small angles in the PSLG. In their algorithm, they first adaptively split the line segments of the PSLG that form a small angle such that their lengths are in the powers of two (in some global scale). During the mesh refinement phase, Miller et al. ignore skinny triangles if the vertices of their shortest edges lie on adjacent line segments of the PSLG that form an angle . As in their algorithm, I too ignore skinny triangles “across” a small angle, but I will define the small angle as as in subsection above.

## 4 The Advancing Front Algorithm

The advancing front algorithm carries out the following three steps in succession to generate a size-optimal mesh:

1. The computation of the piecewise-smooth local feature size functions for the input line segments of the PSLG.

2. The splitting of the input line segments into subsegments whose lengths are asymptotically proportional to the local feature size.

3. The refinement of the truly or the constrained Delaunay triangulation of the PSLG until all skinny triangles are eliminated.

Each step is described in detail in the following subsections. After the input PSLG line segments are split, I will refer to each of the individual split segments as subsegments.

### 4.1 The Feature Size Function Computation

The algorithm requires the knowledge of the local feature size at every point on the input line segments of the PSLG. In order to compute the feature size, let the line segment of PSLG, , be parameterized to lie on the -axis from to , where is its length. I compute the piecewise smooth function that provides the feature size at any point on the line segment. I call this the feature size function of . As mentioned in the previous section, the local feature size at any point on is the distance to the nearest feature that is not adjacent to or the distance to the farthest end point of , whichever is smaller.

Clearly, is the lower envelope of many different functions, which plot the distance to a vertex or a line (nonadjacent features) in the PSLG from . Examples of such distance functions are shown in Fig. 5. The distance from a point on to a vertex in the PSLG is . The function is shown in Fig. 5(a). The distance between line segments (also denoted as ) to some other line segment , whose end points are and , is given by a linear function whose domain is from to , where and are points on such that and are perpendicular to as shown in Fig. 5(b). From to and from to , the distance to vertex and , respectively, defines the distance function. We should also consider the function in Fig. 5(c) that accounts for the farthest end points of to the list of functions over which we compute the lower envelope. The value of this function is at and , and it is at . To compute the feature size function of , we consider the distance function from all vertices and line segments of the PSLG (including end points of the line segments) and the farthest end points of . Note that the distance functions and the feature size function need to be computed only from to .

The lower envelope of the distance function can be computed using a sweep line algorithm that maintains a balanced binary search tree or a heap that orders the functions based on their value at the current location of the sweep line. As I focus on mesh generation in this paper, I direct the readers to a paper  that solves the problem of computing the lower envelope efficiently. The paper is written to compute the lower envelope of lines and line segments, but it can be easily adapted for nonlinear distance functions in our context.

### 4.2 The PSLG Segment Splitting

After the feature size function is computed for each line segment in the PSLG, the next step is to split the line segments such that the length of each subsegment is asymptotically proportional (see Section 3.3 for the definition of asymptotic proportionality) to the feature size at the end points of each subsegment. In order to achieve this goal, I construct a mapping function from a reference segment to the PSLG line segment (see Fig. 6) such that for each point on the reference segment, there is a corresponding point on the PSLG segment (and vice versa). When I split the reference segment evenly and correspondingly split the PSLG line segment at the mapped location, the mapping function ensures that the length of the subsegments in the PSLG line segment is asymptotically proportional to the local feature size. In this section, I will explain how the mapping function is computed by constructing and solving a differential equation. I will also explain how all the reference segments (there is one reference segment for every line segment in the PSLG) are split such that the corresponding splits in the PSLG line segments are size optimal.

#### 4.2.1 Deriving the Differential Equations

Let the mapping function for the line segment be denoted by , where is a point on the reference segment . Let the length of the reference segment be , which is yet to be determined. The mapping function should be designed such that and , where is the length of . Let be split into equal subsegments, which means that we split at , , , , and so on until . The length of each split in the reference segment is given by .

We want the mapping function to result in splits that are asymptotically proportional to the local feature size. Consider a vertex at on . Its corresponding point on is . The vertex adjacent to on is . Its corresponding point on is . The length of this subsegment on is . This length should be proportional to the feature size at , i.e., , where is the local feature size function computed above by constructing the lower envelope of the distance functions. Thus, . If is small enough, we know that . Thus, my intuition is to compute such that . In Section 5, I will show that this intuitive choice of results in asymptotically proportional splits of the PSLG. Note that the feature size function is always positive. Therefore, the mapping function is monotonically increasing as its derivative is also always positive.

#### 4.2.2 Solving the Differential Equations

Let us consider a line segment of the PSLG and the corresponding reference segment . Let the feature size function (obtained in Section 4.1) on that line segment have parts, i.e., there are pieces in the piecewise-smooth function. As seen in Section 4.1, has parts that are either linear or the square root of a quadratic function. Let us denote as . Our differential equation is . When a part of the feature size function is linear, the equation becomes for some and . The solution (see Section 3.3) to this equation is , where needs to be determined using a boundary condition. When the part of the feature size function is the square root of a quadratic function, the equation is of the form for some and . Squaring both sides, we get . Differentiating w.r.t. , we get . Dividing by on both sides, we get . The solution(see Section 3.3) to this equation is , where and need to be determined using two boundary conditions. As is piecewise smooth, is also piecewise smooth, and each part of is given by the solution above.

In the solution to the differential equation provided in Section 3.3, there are constants that need to be evaluated based on the boundary conditions. Let us consider the first part of the feature size function along . When the first part of is linear, we use the value of the mapping function at the initial point. In our case, the initial value at is , i.e., the starting point on maps to the starting point . On the other hand, when is of the form , the corresponding differential equation is of the form , which needs two initial value conditions. The first one is as above. The second condition is given by the local feature size at . Since is equal to the local feature size at , our second boundary condition is .

Using the boundary value conditions provided above, we can analytically compute the first part of the solution of the differential equation. Thus, we have computed the first part of the mapping function from the reference line segment to the PSLG line segment. Let us assume that the first part of the feature size function (on the actual PSLG segment ) starts at and ends at . The first part of the solution to the differential equation (on the reference segment ) starts at and ends at , where . Unfortunately, cannot be analytically computed in all cases, and hence, it needs to be numerically computed in a practical implementation.

Let us assume that the part of the feature size function starts at and ends at . For the second (and subsequent) parts of the solution of the differential equation, the boundary values are given by and for . If there are parts of the feature size function, the length of the reference segment is , where . Note that when the feature size is small, the length of the reference segment is large because the mapping function, whose derivative is proportional to the local feature size, grows slowly when the feature size is small. An example of the construction of the reference segment is provided in the proof of Lemma 5.1.

#### 4.2.3 Splitting the Line Segments

Now that I have computed the mapping functions from every reference segment to the PSLG segment , our task is to split the reference segments evenly so that the PSLG segments are split asymptotically proportional to local feature size. I first split the reference segment with the smallest length into parts, where is determined by the satisfaction of the lemmas in Section 5 (specific lemmas are mentioned in Section 4.4). Let the length of the shortest reference segment be . I then split the reference segment into parts, where

 ni=⌊n∗(t∗)it∗min⌋

and is the length of . Note that when the feature size is small, the length of the reference segment is large, and therefore, the number of splits is also large.

### 4.3 The Off-Center Vertex Insertion

After the PSLG line segments are split into subsegments, I use prior algorithms to obtain a mesh with the desired quality. I use Üngör et al.’s [8, 27] algorithm to place off-center vertices to eliminate skinny triangles from the mesh. As in their algorithm, I prioritize skinny triangles with shortest edges. To reiterate, I will consider the shortest edge in every skinny triangle and pick the triangle with the shortest edge among those considered edges. In the presence of small angles in the input PSLG, I use the algorithm by Pav et al. [13, 15] to decide which skinny triangles to ignore because no algorithm can eliminate all of them. I will show that the splits in Section 4.2 ensure that there is no vertex encroachment if the lemmas in the next section are satisfied.

In order to obtain a high-quality mesh, a truly Delaunay triangulation or a constrained Delaunay triangulation (whichever is desired) of the domain is constructed. Then a skinny triangle (if any) with the shortest edge is chosen. If the end points of the shortest edge of the skinny triangle belong to line segments of the PSLG that form a small angle (defined in Section 3.6), and if the shortest edge is shorter than a certain threshold (explained in Section 5.3), the skinny triangle is ignored. This skinny triangle is considered to be “across” a small angle. If not, its off-center point (see Section 3.4) and the circumcenter of the skinny triangle are considered for insertion into the mesh. Whichever point is closer to the shortest edge of the skinny triangle is inserted, and the domain is retriangulated. Another skinny triangle (if any) with the shortest edge is chosen to be eliminated. These steps are repeated until all skinny triangles are eliminated (except the ones across a small angle).

### 4.4 Satisfaction of Lemmas

For truly Delaunay meshes, in the absence of small angles, the PSLG segments should be split such that the conditions in Lemmas 5.6 and 5.10 are satisfied. In the presence of small angles, in addition to satisfying the conditions in these lemmas, the PSLG should be refined until the Delaunay triangulation of the vertices on the PSLG segments recover the PSLG segments. For constrained Delaunay meshes, the conditions in Lemmas 5.6,  5.12, and 5.14 should be satisfied. For both truly and constrained Delaunay meshes, in the presence of small angles, as we progressively refine the PSLG, the bounds on the minimum and the maximum angle improve.

## 5 An Analysis of the Algorithm

I will recap some of the notations from Section 3 because I use them extensively in the analysis of the algorithm. In my analysis, I will first show that the differential equation-based splitting of the input PSLG line segments will result in size-optimal subsegments such that , where is the length of a subsegment one of whose end points is , and are used to denote the local feature size function, and and are some constants. I will denote the ratio as . A skinny triangle is any triangle whose minimum angle is , where is provided as an input to the algorithm. I will then derive conditions such that no vertex encroaches upon a PSLG subsegment. The conditions are a function of , , , , and , where is the desired minimum radius-edge ratio of triangles in the mesh. I will also show that the algorithm terminates with a size-optimal mesh. In the bound associated with the size optimality, is the maximum distance (normalized to the length of the shortest edge of a skinny triangle) at which a Steiner vertex is placed from the vertices of the shortest edge of a skinny triangle. A small angle is any angle in the input that is . In the initial analysis, I will consider any angle as a small angle. In the appendix, I will show why can be smaller. Finally, I will show that even in the presence of skinny triangles, as long as we refine the PSLG sufficiently, it is possible to obtain truly or constrained Delaunay meshes such that the maximum angle is bounded.

### 5.1 Splitting the PSLG Segments

In the first few lemmas below, I will show that as the PSLG segments are progressively refined, the bounds and increase, but their ratio approaches . In addition, I will show that given an upper bound on or a lower bound on , it is possible to split the PSLG such that is bounded from above.

###### Lemma 5.1.

The length of the shortest reference segment .

###### Proof.

In order to obtain the shortest reference segment, the local feature size at any point on the PSLG line segment of length has to be as large as possible. The feature size cannot be arbitrarily large even if other features in the PSLG are very far away because the feature size is bounded from above by from to and by from to (see Fig. 5(c)), if is assumed to be on the x-axis and is at the origin (no loss of generality).

Let us now derive the mapping function to compute the length of the reference segment. The differential equation for the first piece of the mapping function is given by , where is the mapping function, which implies . The solution to this equation is given by , where is a constant. When , . Therefore, , which implies . Thus, the first piece of the mapping function is . This piece spans from to such that because the first piece of the differential equation spans from to , which implies , which implies . As the other half of the local feature size function is symmetric, the length of the reference line segment is at least . Thus, . ∎

The following lemma applies only to the PSLG segment with the shortest reference segment because the variables and pertain to the PSLG segment. If the two variables are replaced with their equivalent quantities for other segments, the lemma also holds for other segments.

###### Lemma 5.2.

If the PSLG segment with the shortest reference segment whose length is is split into subsegments, the bound on the ratio of the local feature size and length of a subsegment on at some vertex is given by , where is the length of a subsegment one of whose end points is , is the local feature size function, and

 A∗=n∗t∗min−1 and B∗=n∗t∗min+1.
###### Proof.

Let us denote the PSLG line segment under consideration as . In the algorithm, I split the corresponding reference segment into equal parts. The length of each part of the reference segment is . Let be the mapping function. Let us assume that is a vertex on such that it maps to vertex at on , i.e., . The length of one of the segments 222Note that the lemma also holds for the other segment, but I omit that case since the proofs are identical. at is given by because a vertex next to on is at . By the mean value theorem,

 lp=Mi(tp+h)−Mi(tp)=hM′i(tp+h0),

where is some constant. Since is the local feature size function (also denoted as ),

 lp=Mi(tp+h)−Mi(tp)=hF(Mi(tp+h0)). (1)

Note that corresponds to vertex on the PSLG segment, and let correspond to vertex on the PSLG segment. With these vertices, we can apply the property of Lipschitz functions in the next step. Since the local feature size function is a Lipschitz function (see Section 3.2),

 F(Mi(tp+h0))≤F(Mi(tp))+|(Mi(tp+h0)−Mi(tp))|

and

 F(Mi(tp+h0))≥F(Mi(tp))−|(Mi(tp+h0)−Mi(tp))|.

Since the mapping function is a monotonically increasing function,

 F(Mi(tp+h0))≤F(Mi(tp))+(Mi(tp+h0)−Mi(tp))

and

 F(Mi(tp+h0))≥F(Mi(tp))−(Mi(tp+h0)−Mi(tp)).

Substituting the inequalities above into Eq. 1, we get

 lp=hF(Mi(tp+h0))≤h(F(Mi(tp))+(Mi(tp+h0)−Mi(tp)))

and

 lp=hF(Mi(tp+h0))≥h(F(Mi(tp))−(Mi(tp+h0)−Mi(tp))).

Also since (because ), we get

 lp≤h(F(Mi(tp))+lp)

and

 lp≥h(F(Mi(tp))−lp).

Substituting with and rearranging,

 lp(1−h)≤hF(p)≤lp(1+h),

which implies

 (1−h)h≤F(p)lp≤(1+h)h.

Thus,

 A∗=(1−h)h=1−t∗min/n∗t∗min/n∗=n∗−t∗mint∗min=n∗t∗min−1 (2)

and

 B∗=(1+h)h=1+t∗min/n∗t∗min/n∗=n∗+t∗mint∗min=n∗t∗min+1. (3)

The following lemma applies to all PSLG segments. In the proof, I will substitute and (seen in the lemma above) with their equivalents, and , respectively, for any PSLG segment. I have explicitly mentioned about the substitution here so that there is no confusion about the lemmas.

###### Lemma 5.3.

If the PSLG segment with a reference segment of length is split into subsegments, the bound on the ratio of the local feature size and length of the subsegment at some vertex is given by , where is the length of a subsegment one of whose end points is , is the local feature size function, and

 A∗=n∗t∗min−12loge2−1 and B∗=n∗t∗min+1.
###### Proof.

In the proof of Lemma 5.2, Eq. 2 and Eq. 3 were obtained without any assumption about the length of the reference segment being the shortest. Therefore, the lower and upper bounds and for are given by

 A∗=nt∗i−1 and B=nt∗i+1.

Substituting in the above equation,

This equation can be rewritten as

 A∗=n∗t∗it∗min−ϵt∗i−1 and B∗=n∗t∗it∗min−ϵt∗i+1,

where . This equation is equivalent to

 A∗=n∗t∗min−ϵt∗i−1 and B∗=n∗t∗min−ϵt∗i+1.

Since and (by Lemma 5.1),

 A∗≥n∗t∗min−12loge2−1 and B∗≤n∗t∗min+1.

Thus,

 A∗≤LFS(p)lp≤B∗,

where

 A∗=n∗t∗min−12loge2−1 and B∗=n∗t∗min+1.

In the next two lemmas, I will show that it is possible to split the PSLG segments such that there is an upper bound on given a lower bound or an upper bound on the ratio , where is a vertex on the PSLG line segment, and is the length of a subsegment at .

###### Lemma 5.4.

Given a lower bound on , it is possible to split the PSLG line segments such that the upper bound , where is a vertex on the PSLG line segment, and is the length of a subsegment at .

###### Proof.

In Lemma 5.3, we showed that if the shortest reference segment is split into equal parts, the lower bounds on on any PSLG segment are given by

 A∗=nt∗min−12loge2−1,

which is linear in . From this equation, given a lower bound , it is possible to compute the minimum to obtain the lower bound by solving the linear equation. But this might not be an integer. Therefore, we choose the ceiling of , , i.e., we split the shortest reference segment into subsegments in order to obtain the lower bound on the ratio of and . Thus, our bound increases to

 A∗=n+ϵt∗min−12loge2−1,

where . Therefore,

 A∗≥nt∗min−12loge2−1. (4)

Similarly, for the upper bound (from Lemma 5.3),

 B∗=n+ϵt∗min+1=nt∗min+ϵt∗min+1.

Since (Lemma 5.1) and ,

 B∗

If we split the shortest reference segment into subsegments and is as small as it can be (Eq. 4) and is as large as it can be (Eq. 5), the difference between them is, at most, . Therefore, given a lower bound , we can split the PSLG line segments such that the upper bound is at most . ∎

###### Lemma 5.5.

Given an upper bound on the ratio of the upper and lower bound on , where is a vertex on the PSLG line segment and is the length of a subsegment at , it is possible to split the PSLG line segments such that there is an upper bound that is only a function of .

###### Proof.

In Lemma 5.3, we showed that if the shortest reference segment is split into equal parts, the lower bounds on on any PSLG segment are given by

 A∗=nt∗min−12loge2−1,

and the upper bound is given by

 B∗=nt∗min+1=A∗+12loge2+2.

The ratio is given by

 R=1+12A∗loge2+2A∗.

Clearly, it is possible to compute (as a function of only ) and the corresponding for which . The computed may not correspond to being an integer, so we take the ceiling . This operation can only decrease because as tends to infinity, tends to 1. As we saw in Lemma 5.4, the upper bound is bounded for a given lower bound . As we have computed the required lower bound for a given , it is possible to split the PSLG segments such that the upper bound is also bounded as a function of . ∎

### 5.2 Conditions for No Encroachment

In Lemmas 5.4 and 5.5, I have shown that the PSLG line segments can be split in a size-optimal manner if a minimum or a maximum is given. In the next set of lemmas, I will derive the condition on and (as a function of the minimum desired angle or the radius-edge ratio ) such that there is no encroachment of Steiner vertices upon PSLG subsegments. I have already shown that for any such condition, it is possible to obtain a size-optimal split. Further, I show that the final mesh is also size optimal. For now, I assume that there are no small angles in the PSLG. In the next subsection, I will analyze what happens when small angles are present in the PSLG.

#### 5.2.1 Truly Delaunay Refinement

I will first consider mesh refinement that gives us truly Delaunay meshes. In order to obtain truly Delaunay meshes, one should construct the Delaunay triangulation of the vertices inserted by our segment splitting algorithm and recover the PSLG, i.e., the Delaunay triangulation should automatically include all the segments in the PSLG. To achieve this, note that if the diametral circle of every subsegment is empty, the PSLG is recovered by the Delaunay triangulation of the vertices. To see why, consider the midpoint of any subsegment. Its nearest vertices are the two vertices forming the subsegment. Thus, the two vertices are neighbors in the Voronoi diagram, which is the dual of the Delaunay triangulation.

I will first provide the condition on and for which the PSLG is recovered. I will then provide the conditions on for which no Steiner vertices will be inserted in the diametral circle of the PSLG subsegments, which ensures that our algorithm terminates with a size-optimal, high-quality mesh.

###### Lemma 5.6.

If , the diametral circle of any PSLG subsegment (after the split) does not contain vertices from nonadjacent segments.

###### Proof.

If a subsegment is of length and , the feature size at and is at least . Thus, the nearest nonadjacent feature is at least away. The farthest point from or inside the diametral circle of is at a distance away. Thus, the diametral circle does not contain any vertices from nonadjacent features. ∎

For now, I will assume that a small angle in the PSLG is any angle less than . In the appendix, I will show why it is possible to lower the threshold for a small angle to . This small change has limited implications in the proofs in the rest of the paper.

###### Lemma 5.7.

If and the minimum angle , the PSLG is recovered by the Delaunay triangulation of the vertices on subsegments.

###### Proof.

Since vertices on nonadjacent segments are not inside the diametral circle of the PSLG subsegments, and since , the diametral circles of the PSLG subsegment are empty. Thus, the PSLG is recovered by the Delaunay triangulation. ∎

Before I proceed, I will define what I mean by the layer (layer of order ) of the advancing front of Steiner vertices. Let the shortest subsegment of the split PSLG be of length . The vertices on the shortest subsegment of the split PSLG are considered to be a part of the layer of the vertices. For other vertices on the PSLG, consider the shortest PSLG subsegment adjacent to the vertex. Let the length of the subsegment be . If , the vertex is considered to a part of the layer. Note that if there are small angles (see Fig. 7) in the PSLG, the shortest edge adjacent to a vertex on a PSLG subsegment may be a very short edge connecting it to a vertex on an adjacent PSLG segment. Assigning an order to such vertices can be confusing, which is why we assume no small angles are present in this subsection. In the next section, I will elaborate on how to assign their order.

When we insert a Steiner vertex (in the interior, not on a PSLG subsegment) into the mesh, we place it at a distance of from the nearest vertex in the mesh (at the time of insertion). If , I consider to be part of the layer of vertices.

I will define the parent of a Steiner vertex as one of the vertices of the shortest side of the skinny triangle . The vertex is either the off-center vertex of or its circumcenter. Between the two possible vertices of , pick the vertex that is on the lower-order layer. If we move from the Steiner vertex to its parent, then to its grandparent, and so on, we will reach a vertex on a subsegment of the PSLG. Let us call this vertex the ancestral vertex of the Steiner vertex. Si  analyzed Shewchuk’s 3D algorithm  through a similar sequence of vertices. As we move from a Steiner vertex to its ancestral vertex, I will show below that the order of vertices monotonically reduces. Note that we may skip layers as we move from a vertex to one of its children.

###### Lemma 5.8.

The order of a vertex is greater than that of its parent.

###### Proof.

If the length of the shortest segment of a skinny triangle is , one (or both) of the vertices on the shortest segment is at most on the layer . Why? Because the vertex of that was inserted later was inserted at a distance of at most from other vertices in the mesh. The Steiner vertex that is inserted to replace is at least at a distance of from all other vertices in the mesh (because we prioritize skinny triangles with shortest edges first). Thus, belongs to a layer whose order is at least , which is greater than the order of its parent. ∎

###### Lemma 5.9.

If a Steiner vertex belongs to layer and its ancestral vertex on the PSLG belongs to layer , the upper bound on the local feature size at is given by where is the length of the shortest subsegment in the split PSLG.

###### Proof.

Consider the path from to such that every vertex is preceded by its parent. The total length of this path is maximized when the path contains as many vertices as possible and when each edge (from one vertex to the next) has the maximum possible length. This maximization happens when Steiner vertices do not skip a layer (which translates to having as many vertices as possible in the path). Let the path be , , ,…,, . Let the length of the shortest PSLG subsegment at be . Note that . Since the layer order increases by at least , belongs to layer , belongs to layer , and so on. Also, because it belongs to layer , , and so on. Thus, the maximum length of the path is . The local feature size at is bounded by

 fb≤fa+la(α+α2+...αkb−ka),

where is the feature size at and is the feature size at . Since ,

 fb≤B∗la+la(α+α2+...αkb−ka).

Since (it belongs to layer ), where is the length of the shortest side of the PSLG,

 fb≤l0αka(B∗+α+α2+...αkb−ka).