Reference documentation for deal.II version 9.5.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// Copyright (C) 1999 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
19
20#ifdef DEAL_II_WITH_CGAL
21// Functions needed by the CGAL mesh generation utilities are inside
23#endif
24
25
27
28// work around the problem that doxygen for some reason lists all template
29// specializations in this file
30#ifndef DOXYGEN
31
32namespace GridGenerator
33{
34 template <int dim>
35 void
37 const Function<3> & dealii_implicit_function,
39 const Point<3> & interior_point,
40 const double & outer_ball_radius)
41 {
42# ifdef DEAL_II_WITH_CGAL
43 Assert(dealii_implicit_function.n_components == 1,
45 "The implicit function must have exactly one component."));
46 Assert(dealii_implicit_function.value(interior_point) < 0,
48 "The implicit function must be negative at the interior point."));
49 Assert(outer_ball_radius > 0,
50 ExcMessage("The outer ball radius must be positive."));
52 ExcMessage("The triangulation must be empty."));
53
54 if constexpr (dim == 3)
55 {
56 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
57 using NumberType = K::FT;
58 using Point_3 = K::Point_3;
59 using Sphere_3 = K::Sphere_3;
60
61 using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
62 using Tr =
63 CGAL::Mesh_triangulation_3<Mesh_domain,
64 CGAL::Default,
66 using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr, int, int>;
67 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
68
69
70 auto cgal_implicit_function = [&](const Point_3 &p) {
71 return NumberType(
72 dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
73 };
74
75 Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
76 cgal_implicit_function,
77 K::Sphere_3(
78 Point_3(interior_point[0], interior_point[1], interior_point[2]),
79 outer_ball_radius * outer_ball_radius));
80
81 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
82 CGAL::parameters::facet_angle = data.facet_angle,
83 CGAL::parameters::facet_distance =
84 data.facet_distance,
85 CGAL::parameters::cell_radius_edge_ratio =
87 CGAL::parameters::cell_size = data.cell_size);
88
89 auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
90 CGALWrappers::cgal_triangulation_to_dealii_triangulation(
91 cgal_triangulation, tria);
92 }
93 else if constexpr (dim == 2)
94 {
95 // default triangulation for Surface_mesher
96 using Tr = CGAL::Surface_mesh_default_triangulation_3;
97 using C2t3 = CGAL::Complex_2_in_triangulation_3<Tr>;
98 using GT = Tr::Geom_traits;
99 using Sphere_3 = GT::Sphere_3;
100 using Point_3 = GT::Point_3;
101 using FT = GT::FT;
102 typedef FT (*Function)(Point_3);
103 using Surface_3 = CGAL::Implicit_surface_3<GT, Function>;
104 using Surface_mesh = CGAL::Surface_mesh<Point_3>;
105
106
107 auto cgal_implicit_function = [&](const Point_3 &p) {
108 return FT(
109 dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
110 };
111
112 Surface_3 surface(cgal_implicit_function,
113 Sphere_3(Point_3(interior_point[0],
114 interior_point[1],
115 interior_point[2]),
116 outer_ball_radius * outer_ball_radius));
117
118 Tr tr;
119 C2t3 c2t3(tr);
120 Surface_mesh mesh;
121
122 CGAL::Surface_mesh_default_criteria_3<Tr> criteria(data.angular_bound,
123 data.radius_bound,
124 data.distance_bound);
125 CGAL::make_surface_mesh(c2t3,
126 surface,
127 criteria,
128 CGAL::Non_manifold_tag());
129 CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, mesh);
130 CGALWrappers::cgal_surface_mesh_to_dealii_triangulation(mesh, tria);
131 }
132 else
133 {
134 Assert(false, ExcImpossibleInDim(dim));
135 }
136
137# else
138
139 (void)tria;
140 (void)dealii_implicit_function;
141 (void)data;
142 (void)interior_point;
143 (void)outer_ball_radius;
144 AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));
145
146# endif
147 }
148
149
150
151 void
153 Triangulation<3> & vol_tria,
155 {
156# ifdef DEAL_II_WITH_CGAL
157 Assert(
158 surface_tria.n_cells() > 0,
160 "The input triangulation cannot be empty when calling this function."));
161 Assert(
162 vol_tria.n_cells() == 0,
164 "The output triangulation must be empty when calling this function."));
165 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
166 using Point_3 = K::Point_3;
167
168 using Mesh_domain =
169 CGAL::Polyhedral_mesh_domain_with_features_3<K,
170 CGAL::Surface_mesh<Point_3>>;
171 using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,
172 CGAL::Default,
174 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
175 using C3t3 =
176 CGAL::Mesh_complex_3_in_triangulation_3<Tr,
177 Mesh_domain::Corner_index,
178 Mesh_domain::Curve_index>;
179
180 CGAL::Surface_mesh<Point_3> mesh;
181 // This function "fills" the missing arrow of the following diagram.
182 // Tria<2,3> Tria<3>
183 // | ^
184 // | |
185 // | |
186 // | |
187 // V |
188 // CGAL::Surface_mesh -----------> CGAL::C3t3
190 CGAL::Polygon_mesh_processing::triangulate_faces(mesh);
191 CGAL::Polygon_mesh_processing::stitch_borders(mesh);
192 Mesh_domain domain(mesh);
193 domain.detect_features();
194 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
195 CGAL::parameters::facet_angle = data.facet_angle,
196 CGAL::parameters::facet_distance =
197 data.facet_distance,
198 CGAL::parameters::cell_radius_edge_ratio =
200 CGAL::parameters::cell_size = data.cell_size);
201 const auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
202 CGALWrappers::cgal_triangulation_to_dealii_triangulation(cgal_triangulation,
203 vol_tria);
204
205# else
206
207 (void)surface_tria;
208 (void)vol_tria;
209 (void)data;
210 AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));
211
212# endif
213 }
214} // namespace GridGenerator
215
216// explicit instantiations
217# include "grid_generator_cgal.inst"
218
219#endif // DOXYGEN
220
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:112
unsigned int n_active_cells() const
unsigned int n_cells() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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::Surface_mesh< K_exact::Point_3 > Surface_mesh
CGAL::Sequential_tag ConcurrencyTag
Definition utilities.h:83
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)
const ::Triangulation< dim, spacedim > & tria