Quality Meshing for Polyhedra With Small Angles
Quality Meshing for Polyhedra with Small Angles
∗
SiuWing Cheng
†
Tamal K. Dey
‡
Edgar A. Ramos
§
Tathagata Ray
‡
April 11, 2005
Abstract
We present an algorithm to compute a Delaunay mesh conforming to a polyhedron possibly with small input angles. The radiusedge ratio of most output tetrahedra are boundedby a constant, except possibly those that are provably close to small angles. Furthermore,the mesh is not unnecessarily dense in the sense that the edge lengths are at least a constantfraction of the local feature sizes at the edge endpoints. This algorithm is simple to implement as it eliminates most of the computation of local feature sizes and explicit protectivezones. Our experimental results validate that few skinny tetrahedra remain and they lieclose to small acute input angles.
1 Introduction
The need for meshing a polyhedral domain with well shaped tetrahedra arises as an importantproblem in ﬁnite element methods. Mainly there are two approaches known for the problem,
octree based reﬁnement
[1, 13] and the
Delaunay reﬁnement
[7, 9, 11, 12, 16, 17]. Among thesetwo, often the latter is preferred for its directional independence and better quality meshing ingeneral.The Delaunay reﬁnement techniques developed with the work of Chew [6] and Ruppert [16]who showed how to mesh a polygonal domain in 2D with quality triangles. Shewchuk extendedthe technique to polyhedral domains in three dimensions [17]. This extension has two shortcomings. First, it guarantees bounded radiusedge ratio but not bounded aspect ratio. Theradiusedge ratio of a tetrahedron is the ratio of its circumradius to its smallest edge length.Bounded radiusedge ratio eliminates all kinds of poor tetrahedra but not
slivers
. Second, theDelaunay reﬁnement works only with input domains where no angle is smaller than
π
2
. The
sliver exudation
method of Cheng et al. [4] addressed the sliver problem though for unboundeddomains. Chew [7] and Li and Teng [11] proposed nondeterministic point placement to eliminate slivers for bounded domains. Cheng and Dey [3] combined the sliver exudation methodwith the Delaunay reﬁnement to give a deterministic algorithm for meshing bounded domainswith tetrahedra having bounded aspect ratio.
∗
Research is supported by RGC, Hong Kong, China, under grants HKUST6088/99E, HKUST6190/02E,Army Research Oﬃce, USA under grant DAAD190210347 and NSF, USA under grants DMS0310642 andCCR0430735.
†
Department of Computer Science Science, HKUST, Hong Kong. Email:
scheng@cs.ust.hk
‡
Department of Computer Science and Engineering, Ohio State University, Ohio, USA. Email:
tamaldey@cis.ohiostate.edu
and
rayt@cis.ohiostate.edu
§
Department of Computer Science, University of Illinois at UrbanaChampaign, UrbanaChampaign, USA.Email:
eramosn@cs.uiuc.edu
1
However, these works did not address the problem of small angles, speciﬁcally input anglesless than
π
2
. Shewchuk [18] addressed the problem of small angles with constrained Delaunaytriangulation though without any guarantees about the shape quality of the tetrahedra. Thequestion of handling small angles with Delaunay reﬁnement was left open. Murphy, Mount andGable [14] and CohenSteiner, de Verdi`ere and Yvinec [8] presented algorithms that can triangulate any polyhedral domain with Delaunay tetrahedra though without any quality guarantee.Recently, Cheng and Poon [5] proposed a Delaunay meshing algorithm that can handle smallangles with quality guarantees. The output mesh is graded, i.e., for each mesh vertex
v
, itsnearest neighbor distance is Ω(
f
(
v
)) where
f
(
v
) is the local feature size at
v
. The radiusedgeratio of all tetrahedra are bounded by a constant. The constant depends on the smallest input angle in the vicinity of input edges, but it is independent of the input domain elsewhere.Although this algorithm is a signiﬁcant theoretical progress for dealing with small angles, itspractical viability is doubtful. The main reason is that it constructs a protecting region (aunion of balls) around each input edge explicitly, which is complicated and time intensive. Thisconstruction involves computing a feature related distance for a number of points on the edges,intersections among spheres and input facets, and intersections among spheres. Moreover, thesubsequent reﬁnement has to deal with curved edges and spherical surfaces.In this paper we present a new algorithm for meshing polyhedra which has similar guaranteesas that of [5], but it is much simpler to implement. Our algorithm produces a graded Delaunaymesh. The radiusedge ratio of most output tetrahedra are bounded by a constant. For anyremaining tetrahedron of poor shape, all of its vertices are provably close to some small inputangle. In this new algorithm, only input vertices and edges at acute input angles are protected.More importantly, we get rid of the explicit construction of any protecting region, and limit thecomputation of local feature sizes at the input vertices only. Another novelty is that we employthe splitting of nonDelaunay triangles in recovering the input boundary. Earlier works, instead,relied on splitting nonGabriel triangles, that is, the triangles whose smallest circumscribingballs are not empty of vertices. Our experience is that the latter approach may cause manymore points to be inserted. The splitting of nonDelaunay triangles is also essential as we donot have an explicit protecting region. We have implemented our algorithm and we presentexperimental results to demonstrate that our algorithm indeed handles polyhedra with smallangles. The results validate our claim on the quality of the output tetrahedra and show thatfew skinny tetrahedra remain.
2 Preliminaries
Basic deﬁnitions.
We need the following deﬁnitions most of which have been introduced inearlier works.
Skinny tetrahedra.
Let
R
,
ℓ
be the circumradius and shortest edge length respectively of atetrahedron
τ
. Let
ρ
(
τ
) =
R/ℓ
. We say
τ
is
ρ
0

skinny
or simply
skinny
if
ρ
(
τ
)
> ρ
0
.
Input domain.
A vertex
v
is a point in
R
3
; its boundary
∂v
is
v
itself. An
edge
e
is a line segment joining two vertices, say
v
1
and
v
2
. The boundary
∂e
is
{
v
1
,v
2
}
. A
facet
is a subset of a planein
R
3
bounded by a collection of simple polygonal cycles, i.e., for a facet
F
,
∂F
is a collection of edges and vertices forming polygonal cycles. A
piecewise linear complex
PLC is a collection of
elements
which are vertices, edges and facets that form a complex in the following sense. Theintersection of any two elements is either empty or is a collection of lower dimensional elements2
on their boundaries. We assume that the input polyhedron is a subset of
R
3
bounded by a2
manifold
which is the underlying space of a PLC.Two adjacent facets
F
1
and
F
2
may or may not be coplanar. The intersection
F
1
∩
F
2
mayconsist of several connected components. A component may be a single vertex, a single edge,or a polygonal chain of edges. In the case of a single vertex, we call it an
isolated vertex
. If
F
1
and
F
2
are noncoplanar, the components in
F
1
∩
F
2
are collinear and they lie along theintersection line of the support planes of
F
1
and
F
2
.The requirement that any two elements can intersect only at their boundaries disallowsspecial cases such as isolated vertex or a dangling edge in the middle of a facet. We believethat, after adding extra steps to deal with special cases, our algorithm can be extended tohandle any polyhedral domain as long as the union of the elements in the input PLC form a2manifold without boundary.To construct the mesh, we bound this domain in a large enough rectangular box, mesh theinside of the box conforming to the input polyhedron, and only keep the tetrahedra coveringthe interior of the input polyhedron. The introduction of a bounding box has been usedbefore [10, 16]. Let
B
denote the boundary of the bounding box. We use
P
to denote the PLCof the input polyhedron together with
B
. For convenience, sometimes we abuse the notation
P
to denote the underlying space of this input PLC.
Incident and adjacent elements.
Let
∂F
denote the boundary of an element
F
∈ P
. We saythat an element
F
∈ P
is
incident
to an element
G
∈ P
if either
F
⊆
∂G
or
G
⊆
∂F
. Twoelements of
P
are
adjacent
if their intersection is nonempty.
Local feature size.
The local feature size for
P
is a function
f
:
R
3
→
R
where
f
(
x
) is the radiusof the smallest ball centered at
x
intersecting two nonadjacent elements of
P
.
Input angles.
We deﬁne diﬀerent types of angle at the vertices, edges and facets of
P
. Thereare two types of input angles at each vertex
u
of
P
. For any two incident edges of
u
, we measurethe angle between them. We call such angles
edgeedge angles
. For any edge
uv
and a facet
F
incident to
u
such that
uv
is neither incident on
F
nor coplanar with
F
, the angle between
uv
and
F
is min
{
∠
puv
:
p
∈
F,p
=
u
}
. We call such angles
edgefacet angles
. Note that
uv
and
F
are required to be noncoplanar because the angle between
uv
and
F
can be zero otherwise.This happens when the ray emitting from
u
through
v
intersects
F
at another point. At anedge of
P
, we measure the internal and external dihedral angles at the edge.
2
F
1
F
Figure 1: There are two edges and one isolated vertex in
F
1
∩
F
2
.The remaining angles involve two adjacent noncoplanar facets
F
1
and
F
2
. Let
H
1
and
H
2
3
be the support planes of
F
1
and
F
2
, respectively. Notice that
F
i
may cross the line
H
1
∩
H
2
and the intersection of
F
1
and
F
2
may be disconnected; see Figure 1. The planes
H
1
and
H
2
divide
R
3
into four closed regions. We call a region a
wedge
of
F
1
and
F
2
if it contains pointson both
F
1
and
F
2
other than those on the line
H
1
∩
H
2
. The dihedral angle subtended bythe wedge at the line
H
1
∩
H
2
is called its
wedge angle
. Notice that wedges and wedge angleare well deﬁned only when
F
1
and
F
2
are noncoplanar. It is essential to distinguish betweenwedges and nonwedges among the four closed regions induced by
H
1
and
H
2
. Suppose thatwe do not and consider the dihedral angles of all four regions as input angles instead. Theneven a very round convex polyhedron with almost equilateral triangular facets might containinput angles arbitrarily close to zero (when the internal dihedral angle of two adjacent facetsis arbitrarily close to
π
). This is against the intuition that only the vertex angles of the facetsmatter in this case.Throughout this paper, we use
α
m
to denote the minimum input angle in
P
. Notice that
α
m
is strictly positive according to our deﬁnitions. We will need the following lemma aboutthe angles between a line segment and a plane that it is attached to.
Lemma 2.1
Let
ab
be a line segment. Let
H
be a plane that contains
b
. Let
a
′
be the orthogonal projection of
a
onto
H
. Then for any point
x
∈
H
other than
b
,
cos
∠
abx
= cos
∠
aba
′
·
cos
∠
a
′
bx
.Proof.
Let
i
,
j
, and
k
be three orthonormal vectors with srcin at
b
such that
i
lies along
a
′
b
,
j
is normal to
H
and
j
·
ba
≥
0, and
k
·
bx
≥
0. Thencos
∠
abx
=(
a
−
b
)
a
−
b
·
(
x
−
b
)
x
−
b
= (
i
cos
∠
aba
′
+
j
sin
∠
aba
′
)
·
(
i
cos
∠
a
′
bx
+
k
sin
∠
a
′
bx
)= cos
∠
aba
′
·
cos
∠
a
′
bx.
Lemma 2.1 implies that if we rotate a ray
bx
around
b
in
H
starting from the direction of
a
′
−
b
, the angle between
ba
and
bx
increases continuously until the ray reaches the directionof
b
−
a
′
. After that the angle decreases continuously until it returns to the direction of
a
′
−
b
.See Figure 2.
aba’ xc
Figure 2: As the ray
bx
rotates a full circle in clockwise order starting from
ba
′
,
∠
abx
strictlyincreases as
x
goes from
a
′
to
c
, and then strictly decreases as it goes from
c
to
a
′
.
Sharp features.
We call an edge
sharp
if the internal or external dihedral angle at the edgeis acute. Sharp edges can be identiﬁed by taking the dot product of the outward normals of their incident facets. If the dot product is negative, the edge is sharp.4
We call a vertex
u
sharp
if an edgeedge angle or an edgefacet angle at
u
is acute, or if
u
is an endpoint of a sharp edge. Two edges
uv
and
uw
form an acute angle if the dot product(
v
−
u
)
·
(
w
−
u
) is positive. Checking if an edge
uv
forms an acute angle with an incident facet
F
of
u
requires a bit more work. We ﬁrst check that
uv
is not incident to
F
and
uv
is notcoplanar with
F
. Then the angle between
uv
and
F
is acute if and only if
∠
vuw
is acute forsome vertex
w
of
F
.We call a wedge
sharp
if its wedge angle is acute. Our algorithm only needs to identify sharpvertices and sharp edges. That is, our algorithm only works with edgeedge angles, edgefacetangles, and dihedral angles. The wedges, wedge angles, and the minimum input angle
α
m
areused in the analysis only.Notice that, according to our assumption no angle of the bounding box
B
is less than
π
2
though other input angles in
P
may be so.
3 Algorithm
The algorithm has three distinct phases, an initialization phase called
Initialize
, a conformingphase called
Conform
and a reﬁnement phase called
Refine
. In
Initialize
we protect thesharp vertices and compute the initial Delaunay triangulation. In
Conform
, we ﬁrst reﬁneedges and facets so that the 3D Delaunay triangulation conforms to
P
. Then, we determinecertain protecting balls around sharp edges. Some points are disallowed to be inserted in the
Refine
phase in these protecting balls and the vertex balls around sharp vertices. Noticethat there is no explicit construction of a protecting region. After
Conform
,
Refine
startssplitting skinny tetrahedra which may trigger more splitting of the edges and facets. Thealgorithm maintains a vertex set
V
which is updated along with its Delaunay triangulationDel
V
as more vertices are generated.At any generic step of the algorithm, the vertices on any edge of
P
divide it into
subsegments
. The subsegments inside the vertex balls are protected by the vertex balls, so we arenot concerned about them. We call a subsegment
sharp
if it lies outside the vertex balls andon a sharp edge. It is
nonsharp
otherwise. We deﬁne a
circumball
of a subsegment to be aﬁnite ball with the subsegment endpoints on its boundary. The
diametric ball
of a subsegmentis its smallest circumball. A point
p
encroaches
a subsegment if
p
is not an endpoint of thesubsegment and lies on or inside its diametric ball. See Figure 3(a). This deﬁnition of subsegment encroachment is stricter than the usual one as it includes the case where
p
lies on the ballboundary. This helps to simplify proofs later.
pq pq
(
a
) (
b
)Figure 3: In Figure (a),
p
and
q
encroach the subsegment. In Figure (b),
p
encroaches thesubfacet but
q
does not.5