CGAL 4.7 - 2D Conforming Triangulations and Meshes: User Manual

\( \newcommand{\E}{\mathrm{E}} \) \( \newcommand{\A}{\mathrm{A}} \) \( \newcommand{\R}{\mathrm{R}} \) \( \newcommand{\N}{\mathrm{N}} \) \( \newcommand{\Q}{\mathrm{Q}} \) \( \newcommand{\Z}{\mathrm{Z}} \) \( \def\ccSum #1#2#3{ \sum_{#1}^{#2}{#3} } \def\ccProd #1#2#3{ \sum_{#1}^{#2}{#3} }\)

CGAL 4.7 - 2D Conforming Triangulations and Meshes

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
User Manual

Author
  Laurent Rineau

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

# Conforming Triangulations

## Definitions

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 as 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). This edge 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. Conforming Gabriel triangulations are thus also conforming Delaunay triangulations.

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

## Building Conforming Triangulations

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

template<class CDT>
void make_conforming_Delaunay_2 (CDT& t)
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 (see Chapter 2D Triangulations).

The geometric traits of the constrained Delaunay triangulation used to instantiate the parameter `CDT` 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 into a conforming Gabriel triangulation by adding vertices. The user is advised to make a copy of the input triangulation in the case where the original triangulation has to be preserved for other computations

The algorithm used by `make_conforming_Delaunay_2()` and `make_conforming_Gabriel_2()` builds internal data structures 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. For additional control of the refinement algorithm, this class also provides separate functions to insert one Steiner point at a time.

## Example: Making a Triangulation Conforming Delaunay and Then Conforming Gabriel

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 Mesh_2/conforming.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_conformer_2.h>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
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;
}

See Figure 50.1

Figure 50.1 Initial triangulation and the corresponding Delaunay and Gabriel triangulation.

Figure 50.2 The corresponding conforming Delaunay triangulation.

# Meshes

## Definitions

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.

See Figure 50.3 for an example of a domain defined without using seed points, and a possible mesh of it.

Figure 50.3 A domain defined without seed points and the generated mesh.

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 Figure 50.4 for another domain defined with the same Pslg and two seed points used to define holes. In the corresponding mesh these two holes are triangulated but not meshed.

Figure 50.4 A domain with two seeds points defining holes and the generated mesh.

## Shape and Size Criteria

The shape criterion for 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{\frac{1}{2B}}\) on the minimum angle of the triangle and an upper bound of \( \pi - 2* \arcsin{\frac{1}{2B}}\) on the maximum angle. Unfortunately, the termination of the algorithm is guaranteed only if \( B \ge \sqrt{2}\), which corresponds to a lower bound of \( 20.7\) degrees over 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 vary 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 Meshing Algorithm

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. This method inserts new vertices to the triangulation, as far as possible from other vertices, 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 terminate 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 violate the criteria near small input angles. This is unavoidable since small angles formed by input segments cannot be suppressed. Furthermore, it has been shown ([1]), 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.

## Building Meshes

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

template<class CDT, class Criteria>
void refine_Delaunay_mesh_2 (CDT &t, const Criteria& criteria)

The template parameter `CDT` must be instantiated by a constrained Delaunay triangulation class. 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 two models for this concept:

- `Delaunay_mesh_criteria_2<CDT>`, that defines a shape criterion that bounds the minimum angle of triangles,
- `Delaunay_mesh_size_criteria_2<CDT>`, that adds to the previous criterion 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 rebuilds the internal data structure used for meshing at every call. In order to avoid rebuild the data structure at every call, the advanced user can use the class `Delaunay_mesher_2<CDT>`. This class provides also step by step functions. Those functions insert one vertex at a time.

Any object of type `Delaunay_mesher_2<CDT>` is constructed from a reference to a `CDT`, and 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.

Once the mesh is constructed, one can determine which faces of the triangulation are in the mesh domain using the `is_in_domain()` member function of the face type (see the concept `DelaunayMeshFaceBase_2`).

## Optimization of Meshes with Lloyd

The package also provides a global function that runs Lloyd optimization iterations on the mesh generated by Delaunay refinement. The goal of this mesh optimization is to improve the angles inside the mesh, and make them as close as possible to 60 degrees.

template< class CDT >
Mesh_optimization_return_code lloyd_optimize_mesh_2(CDT& cdt);

Note that this global function has more parameters (see details in reference pages) to tune the optimization process.

This optimization process alternates relocating vertices to the center of mass of their Voronoi cells, and updating the Delaunay connectivity of the triangulation. The center of mass is computed with respect to a sizing function that was designed to preserve the local density of points in the mesh generated by Delaunay refinement.

See Figure Figure 50.5 for a mesh generated by `refine_Delaunay_mesh_2()` and optimized with `lloyd_optimize_mesh_2()`. Figure Figure 50.6 shows the histogram of angles inside these meshes.

Figure 50.5 (Left) Mesh generated by `refine_Delaunay_mesh_2()` for a uniform sizing criterion. (Right) shows the same mesh after 100 iterations of Lloyd optimization.

Figure 50.6 Histograms of angles inside the mesh after Delaunay refinement, and after 10 and 100 iterations of Lloyd optimization. After Delaunay refinement, angles are in the interval [28.5; 121.9] degrees. After 10 iterations of Lloyd optimization, they are in [29.1; 110.8]. 100 iterations take them to [29.3; 109.9].

## Examples

### Example Using the Global Function

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 default ones provided by the criteria class `Delaunay_mesh_criteria_2<K>`. No seeds are given, meaning that the mesh domain covers the whole plane except the unbounded component.

File Mesh_2/mesh_global.cpp

#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>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
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;
}

### Example Using the Class Delaunay_mesher_2<CDT>

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 Mesh_2/mesh_class.cpp

#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>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
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;
}

### Example Using Seeds

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 default ones provided by the criteria class `Delaunay_mesh_criteria_2<K>`.

Once the mesh is constructed, the `is_in_domain()` member function of faces is used to count them.

File Mesh_2/mesh_with_seeds.cpp

#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>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
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;
std::cout << "Number of finite faces: " << cdt.number_of_faces() << std::endl;
int mesh_faces_counter = 0;
for(CDT::Finite_faces_iterator fit = cdt.finite_faces_begin();
fit != cdt.finite_faces_end(); ++fit)
{
if(fit->is_in_domain()) ++mesh_faces_counter;
}
std::cout << "Number of faces in the mesh domain: " << mesh_faces_counter << std::endl;
}

### Example Using the Lloyd optimizer

This example uses the global function `lloyd_optimize_mesh_2()`. The mesh is generated using the function `refine_Delaunay_mesh_2()` of `CGAL::Delaunay_mesher_2`, and is then optimized using `lloyd_optimize_mesh_2()`. The optimization will stop after 10 (set by `max_iteration_number`) iterations of alternating vertex relocations and Delaunay connectivity updates. More termination conditions can be used and are detailed in the Reference Manual.

File Mesh_2/mesh_optimization.cpp

#define CGAL_MESH_2_OPTIMIZER_VERBOSE
//#define CGAL_MESH_2_OPTIMIZERS_DEBUG
//#define CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES
#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_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/lloyd_optimize_mesh_2.h>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_mesh_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(-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,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..." << std::endl;
Mesher mesher(cdt);
mesher.set_criteria(Criteria(0.125, 0.05));
mesher.refine_mesh();
std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl;
std::cout << "Run Lloyd optimization...";
CGAL::lloyd_optimize_mesh_2(cdt,
CGAL::parameters::max_iteration_number = 10);
std::cout << " done." << std::endl;
std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl;
}