Delaunay Renement: Mesh Generation

Download as ps, pdf, or txt
Download as ps, pdf, or txt
You are on page 1of 4

Meeting 5

February 2, 1999

Delaunay Re nement
mesh generation, triangle quality, Delaunay re nement,
algorithm, local feature size, constants, invariants, analysis.

Topics:

Mesh generation. The general objective in mesh

would be the largest angle and the aspect ratio. We


argue that a good lower bound for the smallest angle
implies good bounds for the other two expressions of
triangle quality. The largest angle is at most  2,
so if the smallest angle is bounded away from zero then
the largest angle is bounded away from . The aspect
ratio is the length of the longest edge, which we assume
is ac, divided by the distance of b from ac. We have
1  ka ck  2 :
sin 
kb ack sin 
In words, if the smallest angle is bounded away from
zero then the aspect ratio can be neither too small not
too large.
The goal is to construct K so its smallest angle is no
less than some constant, , and the number of triangles
in K is at most some constant times the minimum. We
see from the example in Figure 18 that a small angle
between two input edges cannot possibly be resolved.
A reasonable way to deal with this diculty is to accept sharp input features as unavoidable and to isolate
them so they cause no deterioration of the triangulation
nearby. In these notes we assume there are no sharp input features and in particular that all input angles are
at least 2 .

generation is to decompose a geometric space into elements. The elements are restricted in type and shape,
and the number of elements should not be too big. We
discuss a concrete version of the 2-dimensional mesh
generation problem.
A polygonal region in the plane, possibly with
holes and with constraining edges and vertices inside the region.
Output. A triangulation of the region whose edges
cover all input edges and whose vertices cover all
input vertices.

Input.

The graph of input vertices and edges is denoted by


G, and the output triangulation is denoted by K . It
is convenient to enclose G in a bounding box and to
triangulate everything inside that box. A triangulation
of the input region is obtained by taking a subset of the
triangles. Figure 18 shows input and output of a mesh
generation problem.

Delaunay re nement. We construct K as the De-

launay triangulation of a set of points that includes all


input points. Other points are added one by one to
resolve input edges that are not covered and triangles
that have too small an angle.
1. Suppose ab is a segment of an edge in G that is
not covered by edges of the current Delaunay triangulation. This can only be because some of the
vertices lie inside the diameter circle of ab, see Figure 19. We say these vertices encroach upon ab,
and we use function Split1 to add the midpoint of
ab and to repair the Delaunay triangulations with
a series of ips.

Figure 18: The solid vertices and edges de ne the input

graph, and together with the hollow vertices and dotted


edges they de ne the output triangulation.

Triangle quality. The quality of a triangle abc is ex-

pressed by its smallest angle, . Two alternative choices

12

The choice of B implies that no circumcenter x will ever


lie outside the box. Initially, all 12 triangles next to the
box boundary are non-obtuse so their circumcenters lie
inside B . The algorithm maintains the non-obtuseness
of triangles next to input edges. Since the circumcircles
of Delaunay triangles are empty, this implies that all
circumcenters lie inside B .
It remains to prove that the algorithm halts. Assume for the time being that it does and that n is the
number of vertices in the nal triangulation. The most
expensive step of the algorithm is ipping edges to repair the Delaunay triangulation. The expected linear
upper bound on the number of ips in the randomized
incremental construction does not apply because points
are not added in a random order. The number of ips
per new vertex is bounded by the degree, which is less
than n. The total amount of time to run the algorithm
is therefore in O(n2 ).

Figure 19: Vertex p encroaches upon segment ab. After


adding the midpoint we have too smaller dotted diameter
circles both contained in the diameter circle of ab.

2. Suppose a triangle abc in the current Delaunay triangulation K is skinny, that is, it has an angle less
than . We use function Split2 to add the circumcenter as a new vertex, see Figure 20. Since
the circumcircle of abc is no longer empty, it is
guaranteed to be removed by one of the ips used
to repair the Delaunay triangulation.

Local feature size. We can understand the algorithm better if we relate its actions to the local feature
size de ned as a map f : R ! R. For a point x 2 R ,
2

f (x) is the smallest radius r such that the closed disk


with center x and radius r

x
2

(i) contains two vertices of G,

(ii) intersects one edge of G and contains one vertex of


G that is not endpoint of that edge, or

Figure 20: The angle \axb is twice the angle \acb.

(iii) intersects two vertex disjoint edges of G,

Algorithm. The rst priority of the algorithm is to

see Figure 21. Because of (i) we have f (p)  ka bk


for all vertices a; b 2 G. The local feature size is con-

cover all input edges and the second is to resolve skinny


triangles. Before starting the algorithm we draw the
boundary of a 3-by-3 packing of minimum bounding
rectangles, and we place G inside the center rectangle.
The box is bounded by 12 edges, see Figure 18 where
the box is too small but has the right combinatorics.
Initially, K is the Delaunay triangulation of the input
points, which includes the 12 vertices of B .

loop
while 9

Figure 21: In each case the radius of the circle is the local

encroached upon segment ab do


Split1 (ab)
endwhile;
if no skinny triangle left then exit endif;
let abc 2 K be skinny and x its circumcenter;
x encroaches upon segments s1 ; s2 ; : : :; sk ;
if k  1 then Split1 (si ) for all i
else Split2 (abc)

feature size at x.

tinuous. To see this assume f (x) < f (y) c, for some


positive constant c, but kx yk < c. Then the disk of
x is contained in the interior of the disk of y, which contradicts the de nition of f (y). This argument proves a
Lipschitz condition, which is stronger than continuity.

endif
forever.

Claim 8. jf (x) f (y)j

13

 kx yk.

Constants. The analysis of the algorithm uses two


carefully chosen positive constants C1 and C2 such that
p
1 + 2C2  C1  C2 2sin 1 :
(1)
The constraints that correspond to the two inequalities are bounded by lines, and we have a solution i
the slope ofp the rst line is greater than that of the
second: 1= 2 > 2 sin . Figure 22 illustrates the two
constraints for < 20:7. In this case the two lines
intersect at a point in the positive quadrant. The coordinates of that point are the smallest constants C1 and
C2 that satisfy (1).

Since 1  C2  C1 we have f (a)  L  C1 in all three


cases. Let r = kx ak be the radius of the circumcircle
of abc. Using the Lipschitz condition and sin  = 2Lr
from Figure 20 we get
f (x)  f (a) + r
 L  C1 + r
 2r  sin   C1 + r:
Since  < and C2  1 + 2C1 sin we get
f (x)
f (x)
;

r 
1 + 2C1 sin
C2
as required. We use a similar argument to prove (A),
where x is the midpoint of a segment ab. Let p be a vertex that encroaches upon ab, see Figure 19. Consider
rst the case where p lies on a segment that shares no
endpoint with ab. Then f (x)  r by condition (iii) of
the de nition of local feature size. Consider second the
case where the splitting of ab is triggered by rejecting
the addition of a circumcenter. Let p be this circumcenter and let r0 be the radius of its circle. Since
p p lies
inside the diameter circle of ab we have r0  2r. Using the Lipschitz condition and the argument for (B)
we get
f (x)  f (p) + r
 p
r0  C2 + r
 2r  C2 + r:
p
Using C1  1 + 2C2 we get
f (px)
r 
 fC(x) ;
1 + 2C2
1
as required.

C1
1
1

C2

Figure 22: For < 20:7 the shaded region of pairs


(C1 ; C2) that satisfy (1) is non-empty.

Invariants. We show that with the above constants


the distance between new and old points cannot be
much smaller than the local feature size at the new
point.
Invariant 9. Let p and x be two vertices such that x
is added after p. If x is added by

(A) Split1 then kx pk  f (x)=C1,


(B) Split2 then kx pk  f (x)=C2.
Proof. We use induction to prove (B), where x is the
circumcenter of a skinny triangle abc. Let  < at c be
the smallest angle, see Figure 20. Assume that either a
and b both belong to G or that a is added after b. We
distinguish three cases depending on how a became to
be a vertex. Let L be the length of ab.

Analysis. Invariant 9 guarantees that vertices added

to the triangulation cannot get arbitrarily close to preceding vertices. We show that this implies that they
cannot get close to succeeding vertices either.
Lemma 10. ka bk 

a is a vertex of G. Then b is also a vertex of


G and f (a)  L.

Case 1.

f (a)
1+C1

for all vertices a; b.

Proof. If b precedes a then

ka bk  fC(a)  1f+(aC) :
1
1
Otherwise, we have kb ak  f (b)=C1 and therefore
f (a)  f (b) + ka bk  ka bk  (1 + C1):

a is added as circumcenter of a circle with


radius r0. Before the addition of a the circumcircle
is empty and hence r0  L. By induction we have
f (a)  r0  C2 and therefore f (a)  L  C2.

Case 2.

a is added as the midpoint of a segment. Then


f (a)  L  C1 by the same argument as in Case 2.

Case 3.

14

A simple packing argument implies that the algorithm halts after adding a nite number of vertices. Let
K be the nal triangulation. We relate the number of
triangles to the integral of 1=f 2 (x). Recall that B is
the bounding box used in the construction of K .

[2]

Guaranteed-quality triangular meshes. Report TR-98-983, Comput. Sci. Dept., Cornell Univ.,
Ithaca, New York, 1989.
[3] J. Ruppert. A new and simple algorithm for quality 2-dimensional mesh generation. Report UCB/CSD
92/694, Comput. Sci. Div., Univ. California, Berkeley,
California, 1992.
[4] J. Ruppert. A Delaunay re nement algorithm for quality 2-dimensional mesh generation. J. Algorithms 18
(1995), 548{585.

Lemma 11. The number of vertices in K is at most

some constant times


Z

1 dx:

2
B f (x)

Proof. For each vertex a of K let Da be the disk with


center a and radius ra = f (a)=(2 + 2C1). By Lemma
10, the disks are pairwise disjoint. At least one quarter
of each disk lies inside B . Therefore
Z
XZ
1
dx
dx
 4
2 (x)
2 (x)
f
f
B
a Da
2
X
 41  (f (a)ra+ r )2
a
a
X
 14  (3 + 2C )2 ;
1
a

because (f (a) + ra )2 =ra2  1 + (f (a)=ra )2 = 3 + 2C1.


This is a constant times the number of vertices.
It can also be shown that every triangulation solving
the mesh generation
problem for G has at least some
R
constant times dx=f 2 (x) triangles. This nally shows
that the size of K is at most some constant times the
minimum.

Bibliographic notes. The algorithm in these notes

is due to Jim Ruppert. It achieves best results if the


skinny triangles are removed in order of non-decreasing
smallest angle. The details left out of the journal publication [4] can be found in the technical report [3]. Bern,
Eppstein and Gilbert [1] show that the same technical result (constant minimum angle and proportional to
minimum size) can also be achieved using quad trees.
However, the constants are considerably worse than the
ones obtained using Delaunay re nement. Another predecessor of Ruppert's algorithm is the Delaunay re nement algorithm by Paul Chew [2] that achieves a
constant minimum angle but with uniform density over
the entire region. What distinguishes Ruppert's from
Chew's algorithm is the idea of local feature size and
that it relates to the density of the triangulation.
[1]

M. Bern, D. Eppstein and J. Gilbert. Provably


good mesh generation. J. Comput. Syst. Sci. 48 (1994),
384{409.

15

P. Chew.

You might also like