Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2846-g6fb608615f 2025-03-15 04:10:00+00:00
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
19
20#ifdef DEAL_II_WITH_CGAL
21// Functions needed by the CGAL mesh generation utilities are inside
22# include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
23# include <CGAL/Labeled_mesh_domain_3.h>
24# include <CGAL/Mesh_triangulation_3.h>
26#endif
27
28
30
31// work around the problem that doxygen for some reason lists all template
32// specializations in this file
33#ifndef DOXYGEN
34
35namespace GridGenerator
36{
37 template <int dim>
38 void
40 const Function<3> &dealii_implicit_function,
42 const Point<3> &interior_point,
43 const double &outer_ball_radius)
44 {
45# ifdef DEAL_II_WITH_CGAL
46 Assert(dealii_implicit_function.n_components == 1,
48 "The implicit function must have exactly one component."));
49 Assert(dealii_implicit_function.value(interior_point) < 0,
51 "The implicit function must be negative at the interior point."));
52 Assert(outer_ball_radius > 0,
53 ExcMessage("The outer ball radius must be positive."));
54 Assert(tria.n_active_cells() == 0,
55 ExcMessage("The triangulation must be empty."));
56
57 if constexpr (dim == 3)
58 {
59 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
60 using NumberType = K::FT;
61 using Point_3 = K::Point_3;
62 using Sphere_3 = K::Sphere_3;
63
64 using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
65 using Tr =
66 CGAL::Mesh_triangulation_3<Mesh_domain,
67 CGAL::Default,
69 using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr, int, int>;
70 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
71
72
73 auto cgal_implicit_function = [&](const Point_3 &p) {
74 return NumberType(
75 dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
76 };
77
78 Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
79 cgal_implicit_function,
80 K::Sphere_3(
81 Point_3(interior_point[0], interior_point[1], interior_point[2]),
82 outer_ball_radius * outer_ball_radius));
83
84 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
85 CGAL::parameters::facet_angle = data.facet_angle,
86 CGAL::parameters::facet_distance =
87 data.facet_distance,
88 CGAL::parameters::cell_radius_edge_ratio =
89 data.cell_radius_edge_ratio,
90 CGAL::parameters::cell_size = data.cell_size);
91
92 auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
93 CGALWrappers::cgal_triangulation_to_dealii_triangulation(
94 cgal_triangulation, tria);
95 }
96 else if constexpr (dim == 2)
97 {
98 // default triangulation for Surface_mesher
99 using Tr = CGAL::Surface_mesh_default_triangulation_3;
100 using C2t3 = CGAL::Complex_2_in_triangulation_3<Tr>;
101 using GT = Tr::Geom_traits;
102 using Sphere_3 = GT::Sphere_3;
103 using Point_3 = GT::Point_3;
104 using FT = GT::FT;
105 using Function = FT (*)(Point_3);
106 using Surface_3 = CGAL::Implicit_surface_3<GT, Function>;
107 using Surface_mesh = CGAL::Surface_mesh<Point_3>;
108
109
110 auto cgal_implicit_function = [&](const Point_3 &p) {
111 return FT(
112 dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
113 };
114
115 Surface_3 surface(cgal_implicit_function,
116 Sphere_3(Point_3(interior_point[0],
117 interior_point[1],
118 interior_point[2]),
119 outer_ball_radius * outer_ball_radius));
120
121 Tr tr;
122 C2t3 c2t3(tr);
123 Surface_mesh mesh;
124
125 CGAL::Surface_mesh_default_criteria_3<Tr> criteria(data.angular_bound,
126 data.radius_bound,
127 data.distance_bound);
128 CGAL::make_surface_mesh(c2t3,
129 surface,
130 criteria,
131 CGAL::Non_manifold_tag());
132 CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, mesh);
133 CGALWrappers::cgal_surface_mesh_to_dealii_triangulation(mesh, tria);
134 }
135 else
136 {
137 Assert(false, ExcImpossibleInDim(dim));
138 }
139
140# else
141
142 (void)tria;
143 (void)dealii_implicit_function;
144 (void)data;
145 (void)interior_point;
146 (void)outer_ball_radius;
147 AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));
148
149# endif
150 }
151
152
153
154 void
156 Triangulation<3> &vol_tria,
158 {
159# ifdef DEAL_II_WITH_CGAL
160 Assert(
161 surface_tria.n_cells() > 0,
163 "The input triangulation cannot be empty when calling this function."));
164 Assert(
165 vol_tria.n_cells() == 0,
167 "The output triangulation must be empty when calling this function."));
168 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
169 using Point_3 = K::Point_3;
170
171 using Mesh_domain =
172 CGAL::Polyhedral_mesh_domain_with_features_3<K,
173 CGAL::Surface_mesh<Point_3>>;
174 using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,
175 CGAL::Default,
177 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
178 using C3t3 =
179 CGAL::Mesh_complex_3_in_triangulation_3<Tr,
180 Mesh_domain::Corner_index,
181 Mesh_domain::Curve_index>;
182
183 CGAL::Surface_mesh<Point_3> mesh;
184 // This function "fills" the missing arrow of the following diagram.
185 // Tria<2,3> Tria<3>
186 // | ^
187 // | |
188 // | |
189 // | |
190 // V |
191 // CGAL::Surface_mesh -----------> CGAL::C3t3
193 CGAL::Polygon_mesh_processing::triangulate_faces(mesh);
194 CGAL::Polygon_mesh_processing::stitch_borders(mesh);
195 Mesh_domain domain(mesh);
196 domain.detect_features();
197 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
198 CGAL::parameters::facet_angle = data.facet_angle,
199 CGAL::parameters::facet_distance =
200 data.facet_distance,
201 CGAL::parameters::cell_radius_edge_ratio =
202 data.cell_radius_edge_ratio,
203 CGAL::parameters::cell_size = data.cell_size);
204 const auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
205 CGALWrappers::cgal_triangulation_to_dealii_triangulation(cgal_triangulation,
206 vol_tria);
207
208# else
209
210 (void)surface_tria;
211 (void)vol_tria;
212 (void)data;
213 AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));
214
215# endif
216 }
217} // namespace GridGenerator
218
219// explicit instantiations
220# include "grid/grid_generator_cgal.inst"
221
222#endif // DOXYGEN
223
const unsigned int n_components
Definition function.h:164
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:113
unsigned int n_active_cells() const
unsigned int n_cells() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
Definition mpi.cc:746
void dealii_tria_to_cgal_surface_mesh(const ::Triangulation< dim, spacedim > &triangulation, CGAL::Surface_mesh< CGALPointType > &mesh)
CGAL::Surface_mesh< K_exact::Point_3 > Surface_mesh
CGAL::Sequential_tag ConcurrencyTag
Definition utilities.h:82
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
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)