22 #ifdef DEAL_II_WITH_CGAL
34 # include <CGAL/Boolean_set_operations_2.h>
35 # include <CGAL/Cartesian.h>
36 # include <CGAL/Circular_kernel_intersections.h>
37 # include <CGAL/Constrained_Delaunay_triangulation_2.h>
38 # include <CGAL/Delaunay_mesh_face_base_2.h>
39 # include <CGAL/Delaunay_mesh_size_criteria_2.h>
40 # include <CGAL/Delaunay_mesher_2.h>
41 # include <CGAL/Delaunay_triangulation_2.h>
42 # include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
43 # include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
44 # include <CGAL/Kernel_traits.h>
45 # include <CGAL/Polygon_2.h>
46 # include <CGAL/Polygon_with_holes_2.h>
47 # include <CGAL/Projection_traits_xy_3.h>
48 # include <CGAL/Segment_3.h>
49 # include <CGAL/Simple_cartesian.h>
50 # include <CGAL/Tetrahedron_3.h>
51 # include <CGAL/Triangle_2.h>
52 # include <CGAL/Triangle_3.h>
53 # include <CGAL/Triangulation_2.h>
54 # include <CGAL/Triangulation_3.h>
55 # include <CGAL/Triangulation_face_base_with_id_2.h>
56 # include <CGAL/Triangulation_face_base_with_info_2.h>
61 # include <type_traits>
67 using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt;
68 using K_exact = CGAL::Exact_predicates_exact_constructions_kernel;
69 using K_inexact = CGAL::Exact_predicates_inexact_constructions_kernel;
102 using Vb = CGAL::Triangulation_vertex_base_2<K>;
103 using Fbb = CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>;
104 using CFb = CGAL::Constrained_triangulation_face_base_2<K, Fbb>;
105 using Fb = CGAL::Delaunay_mesh_face_base_2<K, CFb>;
106 using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
107 using Itag = CGAL::Exact_predicates_tag;
108 using CDT = CGAL::Constrained_Delaunay_triangulation_2<K, Tds, Itag>;
109 using Criteria = CGAL::Delaunay_mesh_size_criteria_2<CDT>;
121 template <
typename TargetVariant>
122 struct Repackage : boost::static_visitor<TargetVariant>
124 template <
typename T>
126 operator()(
const T &t)
const
128 return TargetVariant(t);
138 template <
typename... Types>
139 std_cxx17::optional<std_cxx17::variant<Types...>>
140 convert_boost_to_std(
const boost::optional<boost::variant<Types...>> &x)
149 using std_variant = std::variant<Types...>;
150 return boost::apply_visitor(Repackage<std_variant>(), *x);
167 std::list<CDT::Edge> &border)
169 if (start->info().nesting_level != -1)
173 std::list<Face_handle> queue;
174 queue.push_back(start);
175 while (!queue.empty())
179 if (fh->info().nesting_level == -1)
181 fh->info().nesting_level = index;
182 for (
int i = 0; i < 3; i++)
186 if (n->info().nesting_level == -1)
188 if (ct.is_constrained(
e))
205 f->info().nesting_level = -1;
207 std::list<CDT::Edge> border;
209 while (!border.empty())
211 CDT::Edge
e = border.front();
214 if (n->info().nesting_level == -1)
216 mark_domains(cdt, n,
e.first->info().nesting_level + 1, border);
230 std_cxx17::optional<std_cxx17::variant<
CGALPoint2,
233 std::vector<CGALPoint2>>>
235 const std::array<
Point<2>, 3> &second_simplex)
237 std::array<CGALPoint2, 3> pts0, pts1;
239 first_simplex.begin(),
243 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
247 second_simplex.begin(),
248 second_simplex.end(),
251 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
256 return convert_boost_to_std(CGAL::intersection(triangle1, triangle2));
261 std_cxx17::optional<std_cxx17::variant<CGALPoint2, CGALSegment2>>
263 const std::array<
Point<2>, 2> &second_simplex)
265 std::array<CGALPoint2, 3> pts0;
266 std::array<CGALPoint2, 2> pts1;
268 first_simplex.begin(),
272 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
276 second_simplex.begin(),
277 second_simplex.end(),
280 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
285 return convert_boost_to_std(CGAL::intersection(segm, triangle));
291 std::vector<Polygon_with_holes_2>
293 const std::array<
Point<2>, 4> &second_simplex)
295 std::array<CGALPoint2, 4> pts0, pts1;
297 first_simplex.begin(),
301 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
304 second_simplex.begin(),
305 second_simplex.end(),
308 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
313 std::vector<Polygon_with_holes_2> poly_list;
314 CGAL::intersection(
first,
second, std::back_inserter(poly_list));
320 std_cxx17::optional<std_cxx17::variant<CGALPoint3, CGALSegment3>>
322 const std::array<
Point<3>, 4> &second_simplex)
324 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
325 std::array<CGALPoint3, 4> pts0;
326 std::array<CGALPoint3, 2> pts1;
328 first_simplex.begin(),
332 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
337 second_simplex.begin(),
338 second_simplex.end(),
341 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
344 CGALTetra tetra{pts0[0], pts0[1], pts0[2], pts0[3]};
346 return convert_boost_to_std(CGAL::intersection(segm, tetra));
351 "This function requires a version of CGAL greater or equal than 5.5."));
353 (void)second_simplex;
360 std_cxx17::optional<std_cxx17::variant<
CGALPoint3,
363 std::vector<CGALPoint3>>>
365 const std::array<
Point<3>, 4> &second_simplex)
367 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
368 std::array<CGALPoint3, 4> pts0;
369 std::array<CGALPoint3, 3> pts1;
371 first_simplex.begin(),
375 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
379 second_simplex.begin(),
380 second_simplex.end(),
383 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
385 CGALTetra tetra{pts0[0], pts0[1], pts0[2], pts0[3]};
387 return convert_boost_to_std(CGAL::intersection(triangle, tetra));
393 "This function requires a version of CGAL greater or equal than 5.5."));
395 (void)second_simplex;
405 std::vector<std::array<Point<2>, 3>>
407 const std::array<
Point<2>, 4> &vertices0,
408 const std::array<
Point<2>, 4> &vertices1,
411 const auto intersection_test =
414 if (!intersection_test.empty())
416 const auto & poly = intersection_test[0].outer_boundary();
417 const unsigned int size_poly = poly.size();
423 {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)),
424 CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)),
425 CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(2))}}};
427 else if (size_poly >= 4)
430 std::vector<std::array<Point<2>, 3>> collection;
433 cdt.insert_constraint(poly.vertices_begin(),
440 for (
const Face_handle f : cdt.finite_face_handles())
442 if (f->info().in_domain() &&
443 CGAL::to_double(cdt.triangle(f).area()) > tol)
445 collection.push_back(
446 {{CGALWrappers::cgal_point_to_dealii_point<2>(
447 cdt.triangle(f).vertex(0)),
448 CGALWrappers::cgal_point_to_dealii_point<2>(
449 cdt.triangle(f).vertex(1)),
450 CGALWrappers::cgal_point_to_dealii_point<2>(
451 cdt.triangle(f).vertex(2))}});
472 std::vector<std::array<Point<2>, 2>>
474 const std::array<
Point<2>, 4> &vertices0,
475 const std::array<
Point<2>, 2> &vertices1,
478 std::array<CGALPoint2, 4> pts;
480 vertices0.begin(), vertices0.end(), pts.begin(), [&](
const Point<2> &p) {
481 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
487 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(vertices1[0]),
488 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(vertices1[1]));
490 cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(),
true);
491 std::vector<std::array<Point<2>, 2>>
vertices;
495 if (f->info().in_domain() &&
496 CGAL::to_double(cdt.triangle(f).area()) > tol &&
497 CGAL::do_intersect(segm, cdt.triangle(f)))
499 const auto intersection = CGAL::intersection(segm, cdt.triangle(f));
501 boost::get<CGALSegment2>(&*intersection))
504 {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]),
505 CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}});
514 std::vector<std::array<Point<3>, 2>>
516 const std::array<
Point<3>, 8> &vertices0,
517 const std::array<
Point<3>, 2> &vertices1,
520 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
521 std::array<CGALPoint3_exact, 8> pts;
523 vertices0.begin(), vertices0.end(), pts.begin(), [&](
const Point<3> &p) {
524 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
528 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(vertices1[0]),
529 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(vertices1[1]));
533 std::vector<std::array<Point<3>, 2>>
vertices;
535 tria.insert(pts.begin(), pts.end());
536 for (
const auto &c :
tria.finite_cell_handles())
538 const auto &tet =
tria.tetrahedron(c);
539 if (CGAL::do_intersect(segm, tet))
541 const auto intersection = CGAL::intersection(segm, tet);
543 boost::get<CGALSegment3_exact>(&*intersection))
545 if (s->squared_length() > tol * tol)
548 {{CGALWrappers::cgal_point_to_dealii_point<3>(
550 CGALWrappers::cgal_point_to_dealii_point<3>(
562 "This function requires a version of CGAL greater or equal than 5.5."));
571 std::vector<std::array<Point<3>, 3>>
573 const std::array<
Point<3>, 8> &vertices0,
574 const std::array<
Point<3>, 4> &vertices1,
577 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
578 std::array<CGALPoint3_exact, 8> pts_hex;
579 std::array<CGALPoint3_exact, 4> pts_quad;
585 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
593 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
597 std::vector<std::array<Point<3>, 3>>
vertices;
599 tria.insert(pts_hex.begin(), pts_hex.end());
603 tria_quad.insert(pts_quad.begin(), pts_quad.end());
605 for (
const auto &c :
tria.finite_cell_handles())
607 const auto &tet =
tria.tetrahedron(c);
609 for (
const auto &f : tria_quad.finite_facets())
611 if (CGAL::do_intersect(tet, tria_quad.triangle(f)))
613 const auto intersection =
614 CGAL::intersection(tria_quad.triangle(f), tet);
617 boost::get<CGALTriangle3_exact>(&*intersection))
619 if (CGAL::to_double(t->squared_area()) > tol * tol)
622 {{cgal_point_to_dealii_point<3>((*t)[0]),
623 cgal_point_to_dealii_point<3>((*t)[1]),
624 cgal_point_to_dealii_point<3>((*t)[2])}});
628 if (
const std::vector<CGALPoint3_exact> *vps =
629 boost::get<std::vector<CGALPoint3_exact>>(&*intersection))
632 tria_inter.insert(vps->begin(), vps->end());
634 for (
auto it = tria_inter.finite_facets_begin();
635 it != tria_inter.finite_facets_end();
638 const auto triangle = tria_inter.triangle(*it);
639 if (CGAL::to_double(triangle.squared_area()) >
642 std::array<Point<3>, 3> verts = {
643 {CGALWrappers::cgal_point_to_dealii_point<3>(
645 CGALWrappers::cgal_point_to_dealii_point<3>(
647 CGALWrappers::cgal_point_to_dealii_point<3>(
663 "This function requires a version of CGAL greater or equal than 5.5."));
674 std::vector<std::array<Point<3>, 4>>
676 const std::array<
Point<3>, 8> &vertices0,
677 const std::array<
Point<3>, 8> &vertices1,
681 std::array<CGALPoint3_inexact, 8> pts_hex0;
682 std::array<CGALPoint3_inexact, 8> pts_hex1;
688 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
696 return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
702 std::vector<std::array<Point<3>, 4>>
vertices;
705 tria0.insert(pts_hex0.begin(), pts_hex0.end());
706 tria1.insert(pts_hex1.begin(), pts_hex1.end());
708 for (
const auto &c0 : tria0.finite_cell_handles())
710 const auto & tet0 = tria1.tetrahedron(c0);
711 [[maybe_unused]]
const auto &tetg0 =
712 CGAL::make_tetrahedron(tet0.vertex(0),
717 for (
const auto &c1 : tria1.finite_cell_handles())
719 const auto & tet1 = tria1.tetrahedron(c1);
720 [[maybe_unused]]
const auto &tetg1 =
721 CGAL::make_tetrahedron(tet1.vertex(0),
726 namespace PMP = CGAL::Polygon_mesh_processing;
727 const bool test_intersection =
728 PMP::corefine_and_compute_intersection(surf0, surf1, sm);
733 tria.insert(sm.points().begin(), sm.points().end());
734 for (
const auto &c :
tria.finite_cell_handles())
736 const auto &tet =
tria.tetrahedron(c);
738 {{CGALWrappers::cgal_point_to_dealii_point<3>(
740 CGALWrappers::cgal_point_to_dealii_point<3>(
742 CGALWrappers::cgal_point_to_dealii_point<3>(
744 CGALWrappers::cgal_point_to_dealii_point<3>(
758 template <
int dim0,
int dim1,
int spacedim>
759 std::vector<std::array<Point<spacedim>, dim1 + 1>>
772 const auto vertices0 =
773 CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim0)>(
775 const auto vertices1 =
776 CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim1)>(
788 # include "intersections.inst"
803 std::vector<std::array<Point<spacedim>, dim1 + 1>>
815 template <
int dim0,
int dim1,
int spacedim>
816 std::vector<std::array<Point<spacedim>, dim1 + 1>>
Abstract base class for mapping classes.
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNeedsCGAL()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void mark_domains(CDT &ct, Face_handle start, int index, std::list< CDT::Edge > &border)
std_cxx17::optional< std_cxx17::variant< CGALPoint2, CGALSegment2, CGALTriangle2, std::vector< CGALPoint2 > > > compute_intersection(const std::array< Point< 2 >, 3 > &first_simplex, const std::array< Point< 2 >, 3 > &second_simplex)
std::vector< std::array< Point< 3 >, 2 > > compute_intersection_of_cells< 3, 1, 3, 8, 2 >(const std::array< Point< 3 >, 8 > &vertices0, const std::array< Point< 3 >, 2 > &vertices1, const double tol)
K::Tetrahedron_3 CGALTetra
K::Segment_2 CGALSegment2
CGAL::Exact_predicates_tag Itag
K_exact::Tetrahedron_3 CGALTetra_exact
K_exact::Triangle_3 CGALTriangle3_exact
CDT::Vertex_handle Vertex_handle
K::Triangle_3 CGALTriangle3
K::Segment_3 CGALSegment3
CDT::Face_handle Face_handle
CGAL::Triangulation_3< K_inexact > Triangulation3_inexact
K_inexact::Point_3 CGALPoint3_inexact
CGAL::Triangulation_3< K_exact > Triangulation3_exact
CGAL::Constrained_Delaunay_triangulation_2< K, Tds, Itag > CDT
CGAL::Surface_mesh< K_inexact::Point_3 > Surface_mesh
CGAL::Delaunay_mesh_size_criteria_2< CDT > Criteria
CGAL::Triangulation_vertex_base_2< K > Vb
std::vector< std::array< Point< 2 >, 2 > > compute_intersection_of_cells< 2, 1, 2, 4, 2 >(const std::array< Point< 2 >, 4 > &vertices0, const std::array< Point< 2 >, 2 > &vertices1, const double tol)
CGAL::Polygon_with_holes_2< K > Polygon_with_holes_2
std::vector< std::array< Point< 3 >, 3 > > compute_intersection_of_cells< 3, 2, 3, 8, 4 >(const std::array< Point< 3 >, 8 > &vertices0, const std::array< Point< 3 >, 4 > &vertices1, const double tol)
CGAL::Triangulation_data_structure_2< Vb, Fb > Tds
std::vector< std::array< Point< spacedim >, dim1+1 > > compute_intersection_of_cells(const typename Triangulation< dim0, spacedim >::cell_iterator &cell0, const typename Triangulation< dim1, spacedim >::cell_iterator &cell1, const Mapping< dim0, spacedim > &mapping0, const Mapping< dim1, spacedim > &mapping1, const double tol=1e-9)
CGAL::Polygon_2< K > CGALPolygon
CGAL::Triangulation_2< K > Triangulation2
CGAL::Exact_predicates_exact_constructions_kernel K_exact
std::vector< std::array< Point< 2 >, 3 > > compute_intersection_of_cells< 2, 2, 2, 4, 4 >(const std::array< Point< 2 >, 4 > &vertices0, const std::array< Point< 2 >, 4 > &vertices1, const double tol)
CGAL::Triangulation_3< K > Triangulation3
CGAL::Constrained_triangulation_face_base_2< K, Fbb > CFb
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
K::Triangle_2 CGALTriangle2
CGAL::Delaunay_mesh_face_base_2< K, CFb > Fb
K_exact::Segment_3 CGALSegment3_exact
CGAL::Exact_predicates_inexact_constructions_kernel K_inexact
std::vector< std::array< Point< 3 >, 4 > > compute_intersection_of_cells< 3, 3, 3, 8, 8 >(const std::array< Point< 3 >, 8 > &vertices0, const std::array< Point< 3 >, 8 > &vertices1, const double tol)
K_exact::Point_3 CGALPoint3_exact
CGAL::Triangulation_face_base_with_info_2< FaceInfo2, K > Fbb
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr T pow(const T base, const int iexp)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
const ::Triangulation< dim, spacedim > & tria