2D Conforming Triangulations and Meshes

*Laurent Rineau*

This package implements Shewchuk's algorithm [She00] to construct conforming triangulations and 2D meshes. Conforming triangulations will be described in Section and meshes in Section .

A triangulation is a *Delaunay triangulation* if the circumscribing
circle of any facet of the triangulation contains no vertex in its
interior. A *constrained* Delaunay triangulation is a constrained
triangulation which is a much Delaunay as possible. The circumscribing
circle of any facet of a constrained Delaunay triangulation contains in its
interior no data point *visible* from the facet.

An edge is said to be a *Delaunay edge* if it is inscribed in an empty
circle (containing no data point in its interior). It is said to be a
*Gabriel edge* if its diametrical circle is empty.

A constrained Delaunay triangulation is said to be a *conforming
Delaunay triangulation* if every constrained edge is a Delaunay edge.
Because any edge in a constrained Delaunay triangulation is either a
Delaunay edge or a constrained edge, a conforming Delaunay triangulation is
in fact a Delaunay triangulation. The only difference is that some of the
edges are marked as constrained edges.

A constrained Delaunay triangulation is said to be a *conforming
Gabriel triangulation* if every constrained edge is a Gabriel edge. The
Gabriel property is stronger than the Delaunay property and each Gabriel
edge is a Delaunay edge. Thus conforming Gabriel triangulations are also
conforming Delaunay triangulations.

Any constrained Delaunay triangulation can be refined into a
conforming Delaunay triangulation or a conforming Gabriel
triangulation by adding vertices, called *Steiner vertices*, on
constrained edges until they are cut into subconstraints small enough
to be Delaunay or Gabriel edges.

Constrained Delaunay triangulations can be refined into
conforming triangulations
by two global functions:

*template<class CDT> void make_conforming_Delaunay_2 (CDT& t)* and

*template<class CDT> void make_conforming_Gabriel_2 (CDT& t)*.

In both cases, the template parameter *CDT* must be instantiated
by a constrained Delaunay triangulation class. Such a class must be a
model of the concept *ConstrainedDelaunayTriangulation_2*.

There are some requirements on the geometric traits of the constrained
Delaunay triangulation used to instantiate the parameter *CDT*.
It has to be a model of the concept
*ConformingDelaunayTriangulationTraits_2*.

The constrained Delaunay triangulation *t* is passed by reference
and is refined into a conforming Delaunay triangulation or a
conforming Gabriel triangulation by adding vertices, that is, the
triangulation is modified. If the user needs to keep the original
triangulation, he or she has to make a copy of it.

The algorithm used by *make_conforming_Delaunay_2* and
*make_conforming_Gabriel_2* builds internal data that would be
computed twice if the two functions are called consecutively on the same
triangulation. In order to avoid these data to be constructed twice, the
advanced user can use the class *Triangulation_conformer_2<CDT>* to
refine a constrained Delaunay triangulation into a conforming Delaunay
triangulation and then into a conforming Gabriel triangulation. That class
provides also step by step functions. Those functions insert one point at a
time.

This example inserts several segments into a constrained Delaunay triangulation, makes it conforming Delaunay, and then conforming Gabriel. At each step, the number of vertices of the triangulation is printed.

// file: examples/Mesh_2/conforming.C #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Constrained_Delaunay_triangulation_2.h> #include <CGAL/Triangulation_conformer_2.h> #include <iostream> struct K : public CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Constrained_Delaunay_triangulation_2<K> CDT; typedef CDT::Point Point; typedef CDT::Vertex_handle Vertex_handle; int main() { CDT cdt; // construct a constrained triangulation Vertex_handle va = cdt.insert(Point( 5., 5.)), vb = cdt.insert(Point(-5., 5.)), vc = cdt.insert(Point( 4., 3.)), vd = cdt.insert(Point( 5.,-5.)), ve = cdt.insert(Point( 6., 6.)), vf = cdt.insert(Point(-6., 6.)), vg = cdt.insert(Point(-6.,-6.)), vh = cdt.insert(Point( 6.,-6.)); cdt.insert_constraint(va,vb); cdt.insert_constraint(vb,vc); cdt.insert_constraint(vc,vd); cdt.insert_constraint(vd,va); cdt.insert_constraint(ve,vf); cdt.insert_constraint(vf,vg); cdt.insert_constraint(vg,vh); cdt.insert_constraint(vh,ve); std::cout << "Number of vertices before: " << cdt.number_of_vertices() << std::endl; // make it conforming Delaunay CGAL::make_conforming_Delaunay_2(cdt); std::cout << "Number of vertices after make_conforming_Delaunay_2: " << cdt.number_of_vertices() << std::endl; // then make it conforming Gabriel CGAL::make_conforming_Gabriel_2(cdt); std::cout << "Number of vertices after make_conforming_Gabriel_2: " << cdt.number_of_vertices() << std::endl; }

**Figure: **Initial triangulation.

**Figure: **The corresponding conforming Delaunay triangulation.

**Figure: **The corresponding conforming Gabriel triangulation.

A mesh is a partition of a given region into simplices whose shapes and sizes satisfy several criteria.

The domain is the region that the user wants to mesh. It has to be
a bounded region of the plane. The domain is defined by a *planar
straight line graph*, PSLG for short, which is a set of segments
such that two segments in the set are either disjoint or share an
endpoint. The segments of the PSLG are constraints that will be
represented by a union of edges in the mesh. The PSLG can also
contain isolated points that will appear as vertices of the mesh.

The segments of the PSLG are either segments of the boundary or internals constraints. The segments of the PSLG have to cover the boundary of the domain.

The PSLG divides the plane into several connected components. By default, the domain is the union of the bounded connected components. The user can override this default by providing a set of seed points. Either seed points mark components to be meshed or they mark components not to be meshed (holes).

See figures and for an example of a domain defined without using seed points, and a possible mesh of it. See figure for another domain defined with the same PSLG and two seed points. The two seed points define two holes in the domain. In the corresponding mesh (figure ), these two holes are triangulated but not meshed.

**Figure: **A domain defined without seed points.

**Figure: **A mesh of the domain defined without seed points.

**Figure: **A domain with two seeds points defining holes.

**Figure: **A mesh of the domain with seeds defining holes.

The shape criterion on triangles is a lower bound $$*B* on the ratio
between the circumradius and the shortest edge length. Such a bound
implies a lower bound of $$arcsin$$*(1)/(2B)* on the minimum angle
of the triangle and an upper bound of $$*- 2* *arcsin$$*(1)/(2B)*
on the maximum angle. Unfortunately, the termination of the algorithm
is guaranteed only if $$*B sqrt(2)* which corresponds to a lower
bound of $$*20.7* degrees on the angles.

The size criterion can be any criterion that tends to prefer small triangles. For example, the size criterion can be an upper bound on the length of longest edge of triangles, or an upper bound on the radius of the circumcircle. The size bound can be varying over the domain. For example, the size criterion could impose a small size for the triangles intersecting a given line.

Both types of criteria are defined in an object *criteria* passed as
parameter of the meshing functions.

The input to a meshing problem is a PSLG and a set of seeds describing the domain to be meshed, and a set of size and shape criteria. The algorithm implemented in this package starts with a constrained Delaunay triangulation of the input PSLG and produces a mesh using the Delaunay refinement method. That method inserts points into the triangulation, as far as possible from other points, and stops when the criteria are satisfied.

If all angles between incident segments of the input PSLG
are greater than $$*60* degrees and if the bound on the
circumradius/edge ratio is greater than $$*sqrt(2)*,
the algorithm is guaranteed to end up with a mesh
satisfying the size and shape criteria.

If some input angles are smaller than $$*60* degrees, the algorithm will
end up with a mesh in which some triangles near small input angles
violate the criteria. This is unavoidable since small angles formed
by input segments cannot be suppressed. Furthermore, it has been
proven ([She00]), that some domains with small input angles
cannot be meshed with angles even smaller than the small input angles.
Note that if the domain is a polygonal region, the resulting mesh will
satisfy size and shape criteria except for the small input angles.
In addition, the algorithm may succeed in producing meshes with a lower
angle bound greater than $$*20.7* degrees, but there is no such guarantee.

Meshes are obtained from
constrained Delaunay triangulations by calling the global function

*template<class CDT, class Criteria> void refine_Delaunay_mesh_2 (CDT &t, typename CDT::Geom_traits gt)*.

The template parameter *CDT* must be instantiated by a constrained
Delaunay triangulation class that is a model of the concept
*ConstrainedDelaunayTriangulation_2*. In order to override the domain,
a version of this function has two more arguments that define a sequence of
seed points.

The geometric traits class of *CDT* has to be a
model of the concept *DelaunayMeshTraits_2*. This concept
refines the concept *ConformingDelaunayTriangulationTraits_2*
adding the geometric predicates and constructors. The template parameter
*Criteria* must be a model of *MeshingCriteria_2*. This concept
defines criteria that the triangles have to satisfy.
CGAL provides several models for this concept such as:

*Delaunay_mesh_criteria_2<K>*, that defines a shape criterion that bounds the minimum angle of triangles,*Delaunay_mesh_size_criteria<K>*, that adds to the previous one a bound on the maximum edge length.

If the function *refine_Delaunay_mesh_2* is called several times on the
same triangulation with different criteria, the algorithm will rebuild used
internal data at every call. In order to avoid that, the advanced user can
use the class *Delaunay_mesher_2<CDT>*. That class provides also step
by step functions. Those functions insert one point at a time.

Any object of type *Delaunay_mesher_2<CDT>* is constructed from a
reference to a *CDT* and it has several member functions to define the
domain to be meshed and to mesh the *CDT*. See the example given below
and the reference manual for details. Note that the *CDT* should not be
externally modified during the life time of the *Delaunay_mesher_2<CDT>*
object.

The following example inserts several segments into a constrained
triangulation and then meshes it using the global function
*refine_Delaunay_mesh_2*. The size and shape criteria are the defaults
provided by the criteria class *Delaunay_mesh_criteria_2<K>*. No seeds are
given, meaning that the mesh domain covers the whole plane except for the
unbounded component.

// file: examples/Mesh_2/mesh_global.C #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Constrained_Delaunay_triangulation_2.h> #include <CGAL/Delaunay_mesher_2.h> #include <CGAL/Delaunay_mesh_face_base_2.h> #include <CGAL/Delaunay_mesh_size_criteria_2.h> #include <iostream> struct K : public CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Triangulation_vertex_base_2<K> Vb; typedef CGAL::Delaunay_mesh_face_base_2<K> Fb; typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds; typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT; typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria; typedef CDT::Vertex_handle Vertex_handle; typedef CDT::Point Point; int main() { CDT cdt; Vertex_handle va = cdt.insert(Point(-4,0)); Vertex_handle vb = cdt.insert(Point(0,-1)); Vertex_handle vc = cdt.insert(Point(4,0)); Vertex_handle vd = cdt.insert(Point(0,1)); cdt.insert(Point(2, 0.6)); cdt.insert_constraint(va, vb); cdt.insert_constraint(vb, vc); cdt.insert_constraint(vc, vd); cdt.insert_constraint(vd, va); std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl; std::cout << "Meshing the triangulation..." << std::endl; CGAL::refine_Delaunay_mesh_2(cdt, Criteria(0.125, 0.5)); std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl; }

This example uses the class *Delaunay_mesher_2<CDT>* and calls
the *refine_mesh()* member function twice changing the size and
shape criteria in between. In such a case, using twice the global
function *refine_Delaunay_mesh_2* would be less efficient,
because some internal structures needed by the algorithm would be
built twice.

// file: examples/Mesh_2/mesh_class.C #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Constrained_Delaunay_triangulation_2.h> #include <CGAL/Delaunay_mesher_2.h> #include <CGAL/Delaunay_mesh_face_base_2.h> #include <CGAL/Delaunay_mesh_size_criteria_2.h> #include <iostream> struct K : public CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Triangulation_vertex_base_2<K> Vb; typedef CGAL::Delaunay_mesh_face_base_2<K> Fb; typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds; typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT; typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria; typedef CGAL::Delaunay_mesher_2<CDT, Criteria> Mesher; typedef CDT::Vertex_handle Vertex_handle; typedef CDT::Point Point; int main() { CDT cdt; Vertex_handle va = cdt.insert(Point(-4,0)); Vertex_handle vb = cdt.insert(Point(0,-1)); Vertex_handle vc = cdt.insert(Point(4,0)); Vertex_handle vd = cdt.insert(Point(0,1)); cdt.insert(Point(2, 0.6)); cdt.insert_constraint(va, vb); cdt.insert_constraint(vb, vc); cdt.insert_constraint(vc, vd); cdt.insert_constraint(vd, va); std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl; std::cout << "Meshing the triangulation with default criterias..." << std::endl; Mesher mesher(cdt); mesher.refine_mesh(); std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl; std::cout << "Meshing with new criterias..." << std::endl; // 0.125 is the default shape bound. It corresponds to abound 20.6 degree. // 0.5 is the upper bound on the length of the longuest edge. // See reference manual for Delaunay_mesh_size_traits_2<K>. mesher.set_criteria(Criteria(0.125, 0.5)); mesher.refine_mesh(); std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl; }

This example uses the global function *refine_Delaunay_mesh_2* but
defines a domain by using one seed. The size and shape criteria are the
defaults provided by the criteria class *Delaunay_mesh_criteria_2<K>*.

// file: examples/Mesh_2/mesh_with_seeds.C #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Constrained_Delaunay_triangulation_2.h> #include <CGAL/Delaunay_mesher_2.h> #include <CGAL/Delaunay_mesh_face_base_2.h> #include <CGAL/Delaunay_mesh_size_criteria_2.h> #include <iostream> struct K : public CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Triangulation_vertex_base_2<K> Vb; typedef CGAL::Delaunay_mesh_face_base_2<K> Fb; typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds; typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT; typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria; typedef CDT::Vertex_handle Vertex_handle; typedef CDT::Point Point; int main() { CDT cdt; Vertex_handle va = cdt.insert(Point(2,0)); Vertex_handle vb = cdt.insert(Point(0,2)); Vertex_handle vc = cdt.insert(Point(-2,0)); Vertex_handle vd = cdt.insert(Point(0,-2)); cdt.insert_constraint(va, vb); cdt.insert_constraint(vb, vc); cdt.insert_constraint(vc, vd); cdt.insert_constraint(vd, va); va = cdt.insert(Point(3,3)); vb = cdt.insert(Point(-3,3)); vc = cdt.insert(Point(-3,-3)); vd = cdt.insert(Point(3,0-3)); cdt.insert_constraint(va, vb); cdt.insert_constraint(vb, vc); cdt.insert_constraint(vc, vd); cdt.insert_constraint(vd, va); std::list<Point> list_of_seeds; list_of_seeds.push_back(Point(0, 0)); std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl; std::cout << "Meshing the domain..." << std::endl; CGAL::refine_Delaunay_mesh_2(cdt, list_of_seeds.begin(), list_of_seeds.end(), Criteria()); std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl; }

Next chapter: 2D Conforming Triangulations and Meshes

The CGAL Project .
Tue, December 21, 2004 .