Reference documentation for deal.II version GIT d77e5ebb0a 2023-01-27 22:35:02+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\}}\)
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 
16 #include <deal.II/base/ndarray.h>
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 
32 namespace 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,
44  ExcMessage(
45  "The implicit function must have exactly one component."));
46  Assert(dealii_implicit_function.value(interior_point) < 0,
47  ExcMessage(
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."));
51  Assert(tria.n_active_cells() == 0,
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,
159  ExcMessage(
160  "The input triangulation cannot be empty when calling this function."));
161  Assert(
162  vol_tria.n_cells() == 0,
163  ExcMessage(
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
unsigned int n_active_cells() const
unsigned int n_cells() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
#define Assert(cond, exc)
Definition: exceptions.h:1509
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1619
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)
void dealii_tria_to_cgal_surface_mesh(const ::Triangulation< dim, spacedim > &triangulation, CGAL::Surface_mesh< CGALPointType > &mesh)
CGAL::Sequential_tag ConcurrencyTag
Definition: utilities.h:79
CGAL::Surface_mesh< K_inexact::Point_3 > Surface_mesh
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 >{})
const ::Triangulation< dim, spacedim > & tria