Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
grid_generator_cgal.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2022 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
18
19#ifdef DEAL_II_WITH_CGAL
20// Functions needed by the CGAL mesh generation utilities are inside
22#endif
23
24
26
27// work around the problem that doxygen for some reason lists all template
28// specializations in this file
29#ifndef DOXYGEN
30
31namespace GridGenerator
32{
33 template <int dim>
34 void
36 const Function<3> &dealii_implicit_function,
38 const Point<3> &interior_point,
39 const double &outer_ball_radius)
40 {
41# ifdef DEAL_II_WITH_CGAL
42 Assert(dealii_implicit_function.n_components == 1,
44 "The implicit function must have exactly one component."));
45 Assert(dealii_implicit_function.value(interior_point) < 0,
47 "The implicit function must be negative at the interior point."));
48 Assert(outer_ball_radius > 0,
49 ExcMessage("The outer ball radius must be positive."));
50 Assert(tria.n_active_cells() == 0,
51 ExcMessage("The triangulation must be empty."));
52
53 if constexpr (dim == 3)
54 {
55 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
56 using NumberType = K::FT;
57 using Point_3 = K::Point_3;
58 using Sphere_3 = K::Sphere_3;
59
60 using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
61 using Tr =
62 CGAL::Mesh_triangulation_3<Mesh_domain,
63 CGAL::Default,
65 using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr, int, int>;
66 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
67
68
69 auto cgal_implicit_function = [&](const Point_3 &p) {
70 return NumberType(
71 dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
72 };
73
74 Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
75 cgal_implicit_function,
76 K::Sphere_3(
77 Point_3(interior_point[0], interior_point[1], interior_point[2]),
78 outer_ball_radius * outer_ball_radius));
79
80 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
81 CGAL::parameters::facet_angle = data.facet_angle,
82 CGAL::parameters::facet_distance =
83 data.facet_distance,
84 CGAL::parameters::cell_radius_edge_ratio =
86 CGAL::parameters::cell_size = data.cell_size);
87
88 auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
89 CGALWrappers::cgal_triangulation_to_dealii_triangulation(
90 cgal_triangulation, tria);
91 }
92 else if constexpr (dim == 2)
93 {
94 // default triangulation for Surface_mesher
95 using Tr = CGAL::Surface_mesh_default_triangulation_3;
96 using C2t3 = CGAL::Complex_2_in_triangulation_3<Tr>;
97 using GT = Tr::Geom_traits;
98 using Sphere_3 = GT::Sphere_3;
99 using Point_3 = GT::Point_3;
100 using FT = GT::FT;
101 using Function = FT (*)(Point_3);
102 using Surface_3 = CGAL::Implicit_surface_3<GT, Function>;
103 using Surface_mesh = CGAL::Surface_mesh<Point_3>;
104
105
106 auto cgal_implicit_function = [&](const Point_3 &p) {
107 return FT(
108 dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
109 };
110
111 Surface_3 surface(cgal_implicit_function,
112 Sphere_3(Point_3(interior_point[0],
113 interior_point[1],
114 interior_point[2]),
115 outer_ball_radius * outer_ball_radius));
116
117 Tr tr;
118 C2t3 c2t3(tr);
119 Surface_mesh mesh;
120
121 CGAL::Surface_mesh_default_criteria_3<Tr> criteria(data.angular_bound,
122 data.radius_bound,
123 data.distance_bound);
124 CGAL::make_surface_mesh(c2t3,
125 surface,
126 criteria,
127 CGAL::Non_manifold_tag());
128 CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, mesh);
129 CGALWrappers::cgal_surface_mesh_to_dealii_triangulation(mesh, tria);
130 }
131 else
132 {
133 Assert(false, ExcImpossibleInDim(dim));
134 }
135
136# else
137
138 (void)tria;
139 (void)dealii_implicit_function;
140 (void)data;
141 (void)interior_point;
142 (void)outer_ball_radius;
143 AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));
144
145# endif
146 }
147
148
149
150 void
152 Triangulation<3> &vol_tria,
154 {
155# ifdef DEAL_II_WITH_CGAL
156 Assert(
157 surface_tria.n_cells() > 0,
159 "The input triangulation cannot be empty when calling this function."));
160 Assert(
161 vol_tria.n_cells() == 0,
163 "The output triangulation must be empty when calling this function."));
164 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
165 using Point_3 = K::Point_3;
166
167 using Mesh_domain =
168 CGAL::Polyhedral_mesh_domain_with_features_3<K,
169 CGAL::Surface_mesh<Point_3>>;
170 using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,
171 CGAL::Default,
173 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
174 using C3t3 =
175 CGAL::Mesh_complex_3_in_triangulation_3<Tr,
176 Mesh_domain::Corner_index,
177 Mesh_domain::Curve_index>;
178
179 CGAL::Surface_mesh<Point_3> mesh;
180 // This function "fills" the missing arrow of the following diagram.
181 // Tria<2,3> Tria<3>
182 // | ^
183 // | |
184 // | |
185 // | |
186 // V |
187 // CGAL::Surface_mesh -----------> CGAL::C3t3
189 CGAL::Polygon_mesh_processing::triangulate_faces(mesh);
190 CGAL::Polygon_mesh_processing::stitch_borders(mesh);
191 Mesh_domain domain(mesh);
192 domain.detect_features();
193 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
194 CGAL::parameters::facet_angle = data.facet_angle,
195 CGAL::parameters::facet_distance =
196 data.facet_distance,
197 CGAL::parameters::cell_radius_edge_ratio =
199 CGAL::parameters::cell_size = data.cell_size);
200 const auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
201 CGALWrappers::cgal_triangulation_to_dealii_triangulation(cgal_triangulation,
202 vol_tria);
203
204# else
205
206 (void)surface_tria;
207 (void)vol_tria;
208 (void)data;
209 AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));
210
211# endif
212 }
213} // namespace GridGenerator
214
215// explicit instantiations
216# include "grid_generator_cgal.inst"
217
218#endif // DOXYGEN
219
const unsigned int n_components
Definition function.h:163
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:111
unsigned int n_active_cells() const
unsigned int n_cells() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void dealii_tria_to_cgal_surface_mesh(const ::Triangulation< dim, spacedim > &triangulation, CGAL::Surface_mesh< CGALPointType > &mesh)
CGAL::Sequential_tag ConcurrencyTag
Definition utilities.h:82
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
CGAL::Surface_mesh< K_exact::Point_3 > Surface_mesh
void surface_mesh_to_volumetric_mesh(const Triangulation< 2, 3 > &surface_tria, Triangulation< 3 > &vol_tria, const CGALWrappers::AdditionalData< 3 > &data=CGALWrappers::AdditionalData< 3 >{})
void implicit_function(Triangulation< dim, 3 > &tria, const Function< 3 > &implicit_function, const CGALWrappers::AdditionalData< dim > &data=CGALWrappers::AdditionalData< dim >{}, const Point< 3 > &interior_point=Point< 3 >(), const double &outer_ball_radius=1.0)