20#ifdef DEAL_II_WITH_CGAL
22# include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
23# include <CGAL/Labeled_mesh_domain_3.h>
24# include <CGAL/Mesh_triangulation_3.h>
43 const double &outer_ball_radius)
45# ifdef DEAL_II_WITH_CGAL
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."));
55 ExcMessage(
"The triangulation must be empty."));
57 if constexpr (dim == 3)
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;
64 using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
66 CGAL::Mesh_triangulation_3<Mesh_domain,
69 using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr, int, int>;
70 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
73 auto cgal_implicit_function = [&](
const Point_3 &p) {
75 dealii_implicit_function.
value(
Point<3>(p.x(), p.y(), p.z())));
78 Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
79 cgal_implicit_function,
81 Point_3(interior_point[0], interior_point[1], interior_point[2]),
82 outer_ball_radius * outer_ball_radius));
84 Mesh_criteria criteria(CGAL::parameters::facet_size =
data.facet_size,
85 CGAL::parameters::facet_angle =
data.facet_angle,
86 CGAL::parameters::facet_distance =
88 CGAL::parameters::cell_radius_edge_ratio =
89 data.cell_radius_edge_ratio,
90 CGAL::parameters::cell_size =
data.cell_size);
92 auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
93 CGALWrappers::cgal_triangulation_to_dealii_triangulation(
94 cgal_triangulation, tria);
96 else if constexpr (dim == 2)
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;
106 using Surface_3 = CGAL::Implicit_surface_3<GT, Function>;
110 auto cgal_implicit_function = [&](
const Point_3 &p) {
112 dealii_implicit_function.
value(
Point<3>(p.x(), p.y(), p.z())));
115 Surface_3 surface(cgal_implicit_function,
116 Sphere_3(Point_3(interior_point[0],
119 outer_ball_radius * outer_ball_radius));
125 CGAL::Surface_mesh_default_criteria_3<Tr> criteria(
data.angular_bound,
127 data.distance_bound);
128 CGAL::make_surface_mesh(c2t3,
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);
143 (void)dealii_implicit_function;
145 (void)interior_point;
146 (void)outer_ball_radius;
159# ifdef DEAL_II_WITH_CGAL
163 "The input triangulation cannot be empty when calling this function."));
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;
172 CGAL::Polyhedral_mesh_domain_with_features_3<
K,
173 CGAL::Surface_mesh<Point_3>>;
174 using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,
177 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
179 CGAL::Mesh_complex_3_in_triangulation_3<Tr,
180 Mesh_domain::Corner_index,
181 Mesh_domain::Curve_index>;
183 CGAL::Surface_mesh<Point_3> mesh;
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 =
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,
220# include "grid/grid_generator_cgal.inst"
const unsigned int n_components
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
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
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
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)