CGAL 6.2 - Tetrahedral Remeshing

Loading...
Searching...
No Matches
User Manual

Authors
  Jane Tournois, Noura Faraj, Jean-Marc Thiery, Tamy Boubekeur

# Multi-Material Isotropic Tetrahedral Remeshing

This package implements an algorithm for quality tetrahedral remeshing, introduced by Faraj et al in [1]. This practical iterative remeshing algorithm is designed to remesh multi-material tetrahedral meshes, by iteratively performing a sequence of elementary operations such as edge splits, edge collapses, edge flips, and vertex relocations following a Laplacian smoothing. The algorithm results in high-quality uniform isotropic meshes, with the desired mesh density, while preserving the input geometric curve and surface features.

## Remeshing Algorithm

Specific remeshing rules have been designed to satisfy the following criteria. First, the algorithm preserves the geometric complex topology, including multi-material surface patches and polyline features. Polyline features can be defined as intersections between more than two subdomains, or listed by the user. Second, it has been made possible to remesh only a selection of cells, instead of remeshing the whole domain, while preserving or remeshing the interface surfaces between the preserved and the remeshed tetrahedra.

All the local atomic operations that are performed by the algorithm preserve the input topology of the geometric complex.

The tetrahedral remeshing algorithm improves the quality of dihedral angles, while targeting the user-defined uniform sizing field and preserving the topology of the feature complex, as highlighted by Figure 66.1.

Experimental evidence shows that a higher number of remeshing iterations leads to a mesh with a improved fidelity to the sizing criterion, and higher quality dihedral angles.

Figure 66.1 Tetrahedral mesh, modified by our tetrahedral remeshing method. (Left) Before remeshing, dihedral angles were in the interval [1.3; 177.8]. (Right) After remeshing and keeping the same density, dihedral angles are in the interval [12,7; 157.7].

## Sizing Field

Tetrahedral isotropic remeshing algorithm is led by a target sizing field, that can be either uniform, adaptive, or custom.

The sizing field establishes the local target edge length for the remeshed tetrahedra. Two sizing fields are provided: a uniform and a mesh-adaptive sizing field. With `CGAL::Uniform_sizing_field`, all tetrahedra edges are targeted to have equal lengths, as depicted in Section Multi-Material Isotropic Tetrahedral Remeshing . With `CGAL::Adaptive_remeshing_sizing_field`, tetrahedra edge lengths remain locally close to the local size of simplices in the input mesh.

The difference between remeshing with uniform and adaptive sizing fields is depicted in Figure 66.2.

Figure 66.2 Tetrahedral mesh, modified by our tetrahedral remeshing method. (Left) Before remeshing; dihedral angles in the interval [0.1; 179.8], 33,936 vertices. (Middle) After remeshing with uniform sizing field, targeting the average edge length of input (left) mesh; dihedral angles in the interval [9.8; 165.5], 301,127 vertices. (Right) After remeshing with adaptive sizing field; dihedral angles in the interval [6.8; 166.2], 38,602 vertices.

It is also possible to use a custom sizing function, like for example the same sizing field as the one used for initial tetrahedral mesh generation using the CGAL 3D Mesh Generation package, and consecutive consistent isotropic remeshing, as shown by Figure 66.3. Note that the remeshed version is denser than the one generated by `CGAL::make_mesh_3()`. This happens because here, the sizing field is used both as a lower and an upper bound for edge lengths, whereas in mesh generation it is used only as an upper bound on simplex sizes.

Figure 66.3 Tetrahedral mesh, modified by our uniform tetrahedral remeshing method. (Left) Before remeshing, dihedral angles were in the interval `[0.2; 179.7]`. (Middle) After remeshing with a uniform target edge length equal to the average edge length of the input mesh. Dihedral angles are in the interval `[13.2; 157.3]`. (Right) After remeshing with the same sizing field as the one used for mesh generation. Dihedral angles are in the interval `[10.4; 162.3]`.

# API

The tetrahedral remeshing algorithm is implemented as a single free function `CGAL::tetrahedral_isotropic_remeshing()` that takes only two required parameters: the input triangulation, and the desired edge length or sizing field, which drives the remeshing process.

Named Parameters are used to deal with optional parameters. The page Named Parameters describes their usage.

# Examples

## Tetrahedral Remeshing Example

The following example shows the simplest use of the tetrahedral remeshing function. The only required parameter is a given target edge length that drives the remeshing process towards a high-quality tetrahedral mesh with improved dihedral angles, and a more uniform mesh, with edge lengths getting closer to the input parameter value.

By default, the cells with a non-zero `Subdomain_index` are selected for remeshing. The cells with `Subdomain_index` equal to zero are ignored by the remeshing algorithm.

Example: Tetrahedral_remeshing/tetrahedral_remeshing_example.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include "tetrahedral_remeshing_generate_input.h"
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;
int main(int argc, char* argv[])
{
const double target_edge_length = (argc > 1) ? atof(argv[1]) : 0.1;
const std::size_t nbv = (argc > 2) ? atoi(argv[2]) : 1000;
Remeshing_triangulation tr;
CGAL::Tetrahedral_remeshing::insert_random_points_in_cube(nbv, tr);
for (auto cell : tr.finite_cell_handles())
cell->set_subdomain_index(1);
CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length);
return EXIT_SUCCESS;
}
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3
The class Remeshing_triangulation_3 is a class template which provides a valid triangulation type tha...
Definition: Remeshing_triangulation_3.h:59
CGAL::tetrahedral_isotropic_remeshing
void tetrahedral_isotropic_remeshing(CGAL::Triangulation_3< Traits, TDS, SLDS > &tr, const SizingFunction &sizing, const NamedParameters &np=parameters::default_values())
remeshes a tetrahedral mesh.
Definition: tetrahedral_remeshing.h:187

## Tetrahedral Remeshing of a Selection

Optional BGL named parameters offer more precise control on the remeshing process.

The default behaviour is to remesh only the cells with a non-zero `Subdomain_index`. The use of the named parameter `cell_is_selected_map` enables to change this behavior, and select the cells to remesh with a custom selection criterion.

In this example, a triangulation with two subdomains (defined by indices stored in cells) is given as input, but only one (defined by the `Subdomain_index` 2) of its subdomains is remeshed.

Example: Tetrahedral_remeshing/tetrahedral_remeshing_of_one_subdomain.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/property_map.h>
#include "tetrahedral_remeshing_generate_input.h"
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Remeshing_triangulation = CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K>;
using Subdomain_index = Remeshing_triangulation::Tds::Cell::Subdomain_index;
template<typename Tr>
struct Cells_of_subdomain_pmap
{
const int m_subdomain;
public:
using key_type = typename Tr::Cell_handle;
using value_type = bool;
using reference = bool;
using category = boost::read_write_property_map_tag;
friend value_type get(const Cells_of_subdomain_pmap& map,
const key_type& c)
{
return (map.m_subdomain == c->subdomain_index());
}
friend void put(Cells_of_subdomain_pmap&,
const key_type&,
const value_type)
{
; //nothing to do : subdomain indices are updated in remeshing
}
};
Subdomain_index subdomain_on_side_of_plane(const K::Point_3& p,
const K::Plane_3& plane)
{
if (plane.has_on_positive_side(p))
return 1;
else
return 2;
}
K::Point_3 centroid(const Remeshing_triangulation::Cell_handle c)
{
return CGAL::centroid(c->vertex(0)->point(),
c->vertex(1)->point(),
c->vertex(2)->point(),
c->vertex(3)->point());
}
int main(int argc, char* argv[])
{
const double target_edge_length = (argc > 1) ? atof(argv[1]) : 0.05;
const std::size_t nbv = (argc > 2) ? atoi(argv[2]) : 1000;
const std::size_t nbv_on_plane = (argc > 3) ? atoi(argv[3]) : 100;
Remeshing_triangulation tr;
CGAL::Tetrahedral_remeshing::insert_random_points_in_cube(nbv, tr);
const K::Plane_3 plane(0, 0, 1, 0);
CGAL::Tetrahedral_remeshing::insert_points_on_plane(plane, nbv_on_plane, tr);
for (auto cell : tr.finite_cell_handles())
{
const K::Point_3 cc = centroid(cell);
cell->set_subdomain_index(subdomain_on_side_of_plane(cc, plane));
}
CGAL::tetrahedral_isotropic_remeshing(tr,
target_edge_length,
CGAL::parameters::cell_is_selected_map(
Cells_of_subdomain_pmap<Remeshing_triangulation>{2}));
std::ofstream os("out.mesh");
CGAL::IO::write_MEDIT(os, tr);
return EXIT_SUCCESS;
}
CGAL::IO::write_MEDIT
void write_MEDIT(std::ostream &os, const T3 &t3, const NamedParameters &np=parameters::default_values())
CGAL::centroid
CGAL::Point_2< Kernel > centroid(const CGAL::Point_2< Kernel > &p, const CGAL::Point_2< Kernel > &q, const CGAL::Point_2< Kernel > &r)

The following example is a variant of the previous one, where the input triangulation is converted from a `Mesh_complex_3_in_triangulation_3` (e.g., generated by the 3D Mesh Generation package).

Example: Tetrahedral_remeshing/mesh_and_remesh_one_subdomain.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Image_3.h>
#include <CGAL/Mesh_3/Detect_features_in_image.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/property_map.h>
#include <unordered_set>
#include <iostream>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Image_domain = CGAL::Labeled_mesh_domain_3<K>;
using Mesh_domain = CGAL::Mesh_domain_with_polyline_features_3<Image_domain>;
#ifdef CGAL_CONCURRENT_MESH_3
using Concurrency_tag = CGAL::Parallel_tag;
#else
using Concurrency_tag = CGAL::Sequential_tag;
#endif
// Triangulation
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
// Triangulation for Remeshing
using Triangulation_3 = CGAL::Triangulation_3<Tr::Geom_traits, Tr::Triangulation_data_structure>;
// Constraints for protecting 1D-features
using Vertex_handle = Triangulation_3::Vertex_handle;
using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
template <typename Tr>
struct Cells_of_subdomain_pmap
{
using Subdomain_index
= typename Tr::Triangulation_data_structure::Cell::Subdomain_index;
const Subdomain_index m_subdomain;
public:
using key_type = typename Tr::Cell_handle;
using value_type = bool;
using reference = bool;
using category = boost::read_write_property_map_tag;
friend value_type get(const Cells_of_subdomain_pmap& map, const key_type& c) {
return (map.m_subdomain == c->subdomain_index());
}
friend void put(Cells_of_subdomain_pmap&, const key_type&, const value_type) {
; // nothing to do : subdomain indices are updated in remeshing
}
};
// To avoid verbose function and named parameters call
namespace params = CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc>1) ? argv[1] : CGAL::data_file_path("images/420.inr");
CGAL::Image_3 image;
if(!image.read(fname)){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
Mesh_domain domain
= Mesh_domain::create_labeled_image_mesh_domain(image,
params::features_detector = CGAL::Mesh_3::Detect_features_in_image());
const CGAL::Bbox_3 bbox = domain.bbox();
double diag = CGAL::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
CGAL::square(bbox.ymax() - bbox.ymin()) +
CGAL::square(bbox.zmax() - bbox.zmin()));
double sizing_default = diag * 0.05;
Mesh_criteria criteria(params::edge_size = sizing_default,
params::facet_size = sizing_default,
params::facet_distance = sizing_default / 10
);
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::no_exude().no_perturb());
// Output
std::ofstream out0("out_meshed.mesh");
CGAL::IO::write_MEDIT(out0, c3t3.triangulation());
out0.close();
Constraints_set constraints;
Constraints_pmap constraints_pmap(constraints);
Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3),
params::edge_is_constrained_map(constraints_pmap));
const double remeshing_target_size = sizing_default * 0.5;
std::cout << "Remeshing with target edge length: " << remeshing_target_size << std::endl;
CGAL::tetrahedral_isotropic_remeshing(tr, remeshing_target_size,
params::cell_is_selected_map(Cells_of_subdomain_pmap<Triangulation_3>{4})
.number_of_iterations(5)
.edge_is_constrained_map(constraints_pmap)
.remesh_boundaries(false));
// Output
std::ofstream out("out_remeshed.mesh");
CGAL::IO::write_MEDIT(out, tr);
out.close();
return EXIT_SUCCESS;
}
CGAL::Bbox_3
CGAL::Bbox_3::ymin
double ymin() const
CGAL::Bbox_3::xmax
double xmax() const
CGAL::Bbox_3::zmin
double zmin() const
CGAL::Bbox_3::zmax
double zmax() const
CGAL::Bbox_3::ymax
double ymax() const
CGAL::Bbox_3::xmin
double xmin() const
CGAL::Image_3
CGAL::Image_3::read
bool read(const char *file)
CGAL::Labeled_mesh_domain_3
CGAL::Mesh_complex_3_in_triangulation_3
CGAL::Mesh_criteria_3
CGAL::Mesh_domain_with_polyline_features_3
CGAL::Triangulation_3
CGAL::Triangulation_3::Vertex_handle
Triangulation_data_structure::Vertex_handle Vertex_handle
unspecified_type
CGAL::square
NT square(const NT &x)
CGAL::sqrt
NT sqrt(const NT &x)
CGAL::parameters::no_perturb
unspecified_type no_perturb()
CGAL::convert_to_triangulation_3
CGAL::Triangulation_3< typename Tr::Geom_traits, typename Tr::Triangulation_data_structure > convert_to_triangulation_3(CGAL::Mesh_complex_3_in_triangulation_3< Tr, CornerIndex, CurveIndex > c3t3, const NamedParameters &np=parameters::default_values())
converts the triangulation contained in the input to a Triangulation_3.
Definition: tetrahedral_remeshing.h:357
CGAL::parameters
CGAL::Mesh_3::Detect_features_in_image
CGAL::Parallel_tag
CGAL::Sequential_tag

## Tetrahedral Remeshing With Polyline Features

Optional BGL named parameters offer more precise control on the remeshing process. In this example, a triangulation with polyline features that should be preserved - though resampled - during the remeshing process, is given as input. Preserving all surfaces exactly could also be achieved by setting the named parameter `remesh_boundaries` to `false`.

Example: Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/detect_features.h>
#include <CGAL/link_to_face_graph.h>
#include <CGAL/property_map.h>
#include <boost/iterator/function_output_iterator.hpp>
#include <CGAL/IO/File_medit.h>
#include "tetrahedral_remeshing_generate_input.h"
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Remeshing_triangulation = CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K>;
using Point = Remeshing_triangulation::Point;
using Vertex_handle = Remeshing_triangulation::Vertex_handle;
using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
using Surface_mesh = CGAL::Surface_mesh<K::Point_3>;
int main(int argc, char* argv[])
{
const double target_edge_length = (argc > 1) ? atof(argv[1]) : 0.02;
const int nb_iter = (argc > 2) ? atoi(argv[2]) : 1;
const int nbv = (argc > 3) ? atoi(argv[3]) : 500;
Remeshing_triangulation t3;
CGAL::Tetrahedral_remeshing::generate_input_cube(nbv, t3);
// detect sharp edges on the domain boundary and add them as constraints
Surface_mesh boundary;
auto sm_to_t3_vm = boundary.add_property_map<Surface_mesh::Vertex_index, Vertex_handle>("v:v2v").first;
auto v2v_out = boost::make_function_output_iterator([sm_to_t3_vm](const std::pair<Vertex_handle, Surface_mesh::Vertex_index>& p)
{put(sm_to_t3_vm,p.second, p.first);});
CGAL::link_to_face_graph(t3, t3.infinite_vertex(), boundary,
CGAL::parameters::vertex_to_vertex_output_iterator(v2v_out));
auto ecm = boundary.add_property_map<Surface_mesh::Edge_index, bool>("e:ecm", false).first;
CGAL::Polygon_mesh_processing::detect_sharp_edges(boundary, 60., ecm);
Constraints_set constraints;
for (Surface_mesh::Edge_index e : edges(boundary))
if (get(ecm,e))
constraints.emplace(CGAL::make_sorted_pair(get(sm_to_t3_vm,source(e, boundary)),
get(sm_to_t3_vm,target(e, boundary))));
// run remeshing on t3 with constraints
CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length,
CGAL::parameters::edge_is_constrained_map(Constraints_pmap(constraints))
.number_of_iterations(nb_iter));
std::ofstream out("cube_remeshed.mesh");
CGAL::IO::write_MEDIT(out, t3);
return EXIT_SUCCESS;
}
CGAL::Surface_mesh
CGAL::link_to_face_graph
boost::graph_traits< TriangleMesh >::vertex_descriptor link_to_face_graph(const Triangulation_3 &t, typename Triangulation_3::Vertex_handle vh, TriangleMesh &tm, const NamedParameters &np=parameters::default_values())

It is also possible to define the polyline features as the ones stored as complex edges in a `Mesh_complex_3_in_triangulation_3` (e.g., generated by the 3D Mesh Generation package).

Example: Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/property_map.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/IO/File_medit.h>
#include <string>
#include <unordered_set>
// Domain
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Polyhedron = CGAL::Mesh_polyhedron_3<K>::type;
using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K>;
#ifdef CGAL_CONCURRENT_MESH_3
using Concurrency_tag = CGAL::Parallel_if_available_tag;
#else
using Concurrency_tag = CGAL::Sequential_tag;
#endif
// Triangulation for Meshing
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<
Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>;
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
// Triangulation for Remeshing
using Triangulation_3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
using Vertex_handle = Triangulation_3::Vertex_handle;
using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
using Corners_set = std::unordered_set<Vertex_handle, boost::hash<Vertex_handle>>;
using Corners_pmap = CGAL::Boolean_property_map<Corners_set>;
// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/fandisk.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if (input.fail() || !CGAL::is_triangle_mesh(polyhedron)) {
std::cerr << "Error: Input invalid " << fname << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
// Get sharp features
domain.detect_features();
// Mesh criteria
Mesh_criteria criteria(edge_size = 0.025,
facet_angle = 25, facet_size = 0.05, facet_distance = 0.005,
cell_radius_edge_ratio = 3, cell_size = 0.05);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Property map of constraints
Constraints_set constraints;
Constraints_pmap constraints_pmap(constraints);
Corners_set corners;
Corners_pmap corners_pmap(corners);
Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3),
edge_is_constrained_map(constraints_pmap).
vertex_is_constrained_map(corners_pmap));
//note we use the move semantic, with std::move(c3t3),
// to avoid a copy of the triangulation by the function
// `CGAL::convert_to_triangulation_3()`
// After the call to this function, c3t3 is an empty and valid C3t3.
//It is possible to use : CGAL::convert_to_triangulation_3(c3t3),
// Then the triangulation is copied and duplicated, and c3t3 remains as is.
const double target_edge_length = 0.1;//coarsen the mesh
CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length,
number_of_iterations(5)
.smooth_constrained_edges(true)
.edge_is_constrained_map(constraints_pmap));
std::ofstream out("out_remeshed.mesh");
CGAL::IO::write_MEDIT(out, tr);
out.close();
return EXIT_SUCCESS;
}
CGAL::Polyhedral_mesh_domain_with_features_3
CGAL::is_triangle_mesh
bool is_triangle_mesh(const FaceGraph &g)
CGAL::Parallel_if_available_tag

## Tetrahedral Remeshing After Mesh Generation

The tetrahedral remeshing algorithm is designed as a postprocessing for mesh generation algorithms. The API allows to generate a tetrahedral mesh with the CGAL 3D Mesh Generation package, and further improve it with the tetrahedral remeshing algorithm. This example shows how to use tetrahedral mesh generation and remeshing in sequence, from a polyhedral domain with features.

Example: Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/property_map.h>
#include <CGAL/IO/File_medit.h>
#include <string>
#include <unordered_set>
// Domain
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Surface_mesh = CGAL::Surface_mesh<K::Point_3>;
using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K, Surface_mesh>;
#ifdef CGAL_CONCURRENT_MESH_3
using Concurrency_tag = CGAL::Parallel_if_available_tag;
#else
using Concurrency_tag = CGAL::Sequential_tag;
#endif
// Triangulation
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<
Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>;
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
// Triangulation for Remeshing
using Triangulation_3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
using Vertex_handle = Triangulation_3::Vertex_handle;
using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/fandisk.off");
const int nb_iter = (argc > 2) ? atoi(argv[2]) : 5;
std::ifstream input(fname);
Surface_mesh mesh;
std::string filename(fname);
input >> mesh;
if (input.fail()) {
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
if (!CGAL::is_triangle_mesh(mesh)) {
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(mesh);
// Get sharp features
domain.detect_features();
// Mesh criteria
const double size = 0.072;
Mesh_criteria criteria(edge_size = size,
facet_angle = 25,
facet_size = size,
facet_distance = 0.1 * size,
cell_radius_edge_ratio = 2,
cell_size = size);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude());
Constraints_set constraints;
Constraints_pmap constraints_pmap(constraints);
Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3),
edge_is_constrained_map(constraints_pmap));
// Remeshing
CGAL::tetrahedral_isotropic_remeshing(tr, size,
number_of_iterations(nb_iter)
.edge_is_constrained_map(constraints_pmap));
std::ofstream out("out_remeshed.mesh");
CGAL::IO::write_MEDIT(out, tr);
out.close();
return EXIT_SUCCESS;
}
CGAL::parameters::no_exude
unspecified_type no_exude()

## Tetrahedral Remeshing With Sizing Field

The API enables the user to run tetrahedral remeshing with a custom sizing field, that is used as a sizing criterion to trigger mesh refinement and decimation. The following examples show how to use this type of sizing field.

In the first example, we show how to use the same sizing field as the one used as a meshing criterion for tetrahedral mesh generation using the CGAL 3D Mesh Generation package. Example: Tetrahedral_remeshing/mesh_and_remesh_with_sizing.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_if_available_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT (Function)(const Point&);
typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
// Triangulation for Remeshing
using T3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
namespace params = CGAL::parameters;
// Sizing field
struct Spherical_sizing_field
{
typedef K::FT FT;
typedef K::Point_3 Point_3;
typedef Mesh_domain::Index Index;
FT operator()(const Point_3& p, const int, const Index&) const
{
FT sq_d_to_origin = CGAL::squared_distance(p, Point(CGAL::ORIGIN));
return CGAL::abs(CGAL::sqrt(sq_d_to_origin) - 0.5) / 5. + 0.025;
}
};
// Function
FT sphere_function (const Point& p)
{
return CGAL::squared_distance(p, Point(CGAL::ORIGIN)) - 1;
}
int main()
{
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain
(sphere_function, K::Sphere_3(CGAL::ORIGIN, K::FT(2))
);
// Mesh criteria
Spherical_sizing_field size;
Mesh_criteria criteria(params::facet_angle(30).facet_size(0.1).facet_distance(0.025).
cell_radius_edge_ratio(2).cell_size(size));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_exude().no_perturb());
T3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3));
//note we use the move semantic, with std::move(c3t3),
// to avoid a copy of the triangulation by the function
// `CGAL::convert_to_triangulation_3()`
// After the call to this function, c3t3 is an empty and valid C3t3.
//It is possible to use : CGAL::convert_to_triangulation_3(c3t3),
// Then the triangulation is copied and duplicated, and c3t3 remains as is.
std::cout << "Remeshing...";
std::cout.flush();
CGAL::tetrahedral_isotropic_remeshing(tr, size);
std::cout << "\rRemeshing done." << std::endl;
return EXIT_SUCCESS;
}
CGAL::Point_3
Index
CGAL::abs
NT abs(const NT &x)
CGAL::ORIGIN
const CGAL::Origin ORIGIN
CGAL::squared_distance
Kernel::FT squared_distance(Type1< Kernel > obj1, Type2< Kernel > obj2)

In the second example, we show how to use the provided adaptive sizing field, `CGAL::Adaptive_remeshing_sizing_field`, that automatically sets the target edge length to keep a mesh density consistent with the density of the input mesh. Example: Tetrahedral_remeshing/mesh_and_remesh_with_adaptive_sizing.cpp

Show / Hide
#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_3.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/Adaptive_remeshing_sizing_field.h>
#include <CGAL/IO/File_medit.h>
#include <string>
// Domain
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point = K::Point_3;
using FT = K::FT;
using Polyhedron = CGAL::Polyhedron_3<K>;
using Mesh_domain = CGAL::Polyhedral_mesh_domain_3<Polyhedron, K>;
// Triangulation for Meshing
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
// Triangulation for Remeshing
using T3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elk.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if (input.fail()) {
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
if (!CGAL::is_triangle_mesh(polyhedron)) {
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
std::cout << "Meshing...";
std::cout.flush();
// Mesh criteria
Mesh_criteria criteria(facet_angle = 25,
facet_distance = 0.2,
cell_size = 10.,
cell_radius_edge_ratio = 3.);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb().no_exude());
std::cout << "\rMeshing done." << std::endl;
std::ofstream out_mesh("out_meshing.mesh");
CGAL::IO::write_MEDIT(out_mesh, c3t3.triangulation());
out_mesh.close();
T3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3));
//note we use the move semantic, with std::move(c3t3),
// to avoid a copy of the triangulation by the function
// `CGAL::convert_to_triangulation_3()`
// After the call to this function, c3t3 is an empty and valid C3t3.
//It is possible to use : CGAL::convert_to_triangulation_3(c3t3),
// Then the triangulation is copied and duplicated, and c3t3 remains as is.
std::cout << "Remeshing...";
std::cout.flush();
CGAL::tetrahedral_isotropic_remeshing(tr,
CGAL::create_adaptive_remeshing_sizing_field(tr),
CGAL::parameters::number_of_iterations(5));
std::cout << "\rRemeshing done." << std::endl;
std::ofstream out("out_remeshing.mesh");
CGAL::IO::write_MEDIT(out, tr);
out.close();
return EXIT_SUCCESS;
}
CGAL::Polyhedral_mesh_domain_3
CGAL::Polyhedron_3

## Tetrahedral Remeshing from Any Tetrahedral Mesh

The following example shows how to read a mesh from a triangulation stored in a Medit file, perform tetrahedral remeshing, and save the output triangulation. The input triangulation should follow the validity requirements of a `CGAL::Triangulation_3` (valid connectivity, positive orientation of the cells, and coverage of the convex hull of the vertices).

Example: Tetrahedral_remeshing/tetrahedral_remeshing_from_mesh.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/IO/File_medit.h>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;
int main(int argc, char* argv[])
{
const char* filename = (argc > 1) ? argv[1] : "data/sphere.mesh";
const double target_edge_length = (argc > 2) ? atof(argv[2]) : 0.1;
Remeshing_triangulation tr;
std::ifstream is(filename, std::ios_base::in);
CGAL::IO::read_MEDIT(is, tr);
CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length);
std::ofstream os("after_remeshing.mesh");
CGAL::IO::write_MEDIT(os, tr);
return EXIT_SUCCESS;
}
CGAL::IO::read_MEDIT
bool read_MEDIT(std::istream &in, T3 &t3, const NamedParameters &np=parameters::default_values())

# Implementation History

This package implements the uniform version of the "Multi-Material Adaptive Volume Remesher" algorithm for quality tetrahedral remeshing, described by Noura Faraj et al. in [1].

A first version of the code was written by Noura Faraj, Jean-Marc Thiery, and Tamy Boubekeur. Jane Tournois worked on the finalization of the code, the API, and documentation.

It was initially published in CGAL-5.1.

The ability to use a sizing field as variable target edge length was introduced by Jane Tournois in CGAL-6.0.