21#ifdef DEAL_II_WITH_CGAL
31# include <CGAL/Boolean_set_operations_2.h>
34# include <CGAL/Cartesian.h>
35# include <CGAL/Circular_kernel_intersections.h>
36# include <CGAL/Constrained_Delaunay_triangulation_2.h>
37# include <CGAL/Delaunay_mesh_face_base_2.h>
38# include <CGAL/Delaunay_mesh_size_criteria_2.h>
39# include <CGAL/Delaunay_mesher_2.h>
40# include <CGAL/Delaunay_triangulation_2.h>
41# include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
42# include <CGAL/Kernel_traits.h>
43# include <CGAL/Polygon_2.h>
44# include <CGAL/Polygon_with_holes_2.h>
45# include <CGAL/Projection_traits_xy_3.h>
46# include <CGAL/Segment_3.h>
47# include <CGAL/Simple_cartesian.h>
48# include <CGAL/Surface_mesh/Surface_mesh.h>
49# include <CGAL/Tetrahedron_3.h>
50# include <CGAL/Triangle_2.h>
51# include <CGAL/Triangle_3.h>
52# include <CGAL/Triangulation_2.h>
53# include <CGAL/Triangulation_3.h>
54# include <CGAL/Triangulation_face_base_with_id_2.h>
55# include <CGAL/Triangulation_face_base_with_info_2.h>
65 using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt;
66 using K_exact = CGAL::Exact_predicates_exact_constructions_kernel;
96 using Vb = CGAL::Triangulation_vertex_base_2<K>;
97 using Fbb = CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>;
98 using CFb = CGAL::Constrained_triangulation_face_base_2<K, Fbb>;
99 using Fb = CGAL::Delaunay_mesh_face_base_2<K, CFb>;
100 using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
101 using Itag = CGAL::Exact_predicates_tag;
102 using CDT = CGAL::Constrained_Delaunay_triangulation_2<K, Tds, Itag>;
103 using Criteria = CGAL::Delaunay_mesh_size_criteria_2<CDT>;
107 template <
class T,
class... Types>
111 return std::get_if<T>(v);
114 template <
class T,
class... Types>
118 return boost::get<T>(v);
129 template <
typename TargetVariant>
130 struct Repackage : boost::static_visitor<TargetVariant>
132 template <
typename T>
134 operator()(
const T &t)
const
136 return TargetVariant(t);
146 template <
typename... Types>
147 std::optional<std::variant<Types...>>
148 convert_boost_to_std(
const boost::optional<boost::variant<Types...>> &x)
157 using std_variant = std::variant<Types...>;
158 return boost::apply_visitor(Repackage<std_variant>(), *x);
168 template <
typename... Types>
169 const std::optional<std::variant<Types...>> &
170 convert_boost_to_std(
const std::optional<std::variant<Types...>> &opt)
182 std::list<CDT::Edge> &border)
184 if (start->info().nesting_level != -1)
188 std::list<Face_handle> queue;
189 queue.push_back(start);
190 while (!queue.empty())
194 if (fh->info().nesting_level == -1)
196 fh->info().nesting_level = index;
197 for (
int i = 0; i < 3; i++)
201 if (n->info().nesting_level == -1)
203 if (ct.is_constrained(e))
218 for (CDT::Face_handle f : cdt.all_face_handles())
220 f->info().nesting_level = -1;
222 std::list<CDT::Edge> border;
224 while (!border.empty())
226 CDT::Edge e = border.front();
229 if (n->info().nesting_level == -1)
231 mark_domains(cdt, n, e.first->info().nesting_level + 1, border);
248 std::vector<CGALPoint2>>>
256 std::array<CGALPoint2, 3> pts0, pts1;
258 std::transform(triangle0.begin(),
261 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
263 std::transform(triangle1.begin(),
266 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
270 return convert_boost_to_std(
271 CGAL::intersection(cgal_triangle0, cgal_triangle1));
275 std::optional<std::variant<CGALPoint2, CGALSegment2>>
283 std::array<CGALPoint2, 3> pts0;
284 std::array<CGALPoint2, 2> pts1;
286 std::transform(triangle.begin(),
289 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
291 std::transform(segment.begin(),
294 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
298 return convert_boost_to_std(
299 CGAL::intersection(cgal_segment, cgal_triangle));
305 std::vector<Polygon_with_holes_2>
312 std::array<CGALPoint2, 4> pts0, pts1;
314 std::transform(rectangle0.begin(),
317 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
319 std::transform(rectangle1.begin(),
322 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
324 const CGALPolygon first_poly{pts0.begin(), pts0.end()};
325 const CGALPolygon second_poly{pts1.begin(), pts1.end()};
327 std::vector<Polygon_with_holes_2> poly_list;
328 CGAL::intersection(first_poly,
330 std::back_inserter(poly_list));
336 std::optional<std::variant<CGALPoint3, CGALSegment3>>
341# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
346 std::array<CGALPoint3, 4> pts0;
347 std::array<CGALPoint3, 2> pts1;
349 std::transform(tetrahedron.begin(),
352 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
354 std::transform(segment.begin(),
357 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
359 CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
361 return convert_boost_to_std(
362 CGAL::intersection(cgal_segment, cgal_tetrahedron));
367 "This function requires a version of CGAL greater or equal than 5.5."));
379 std::vector<CGALPoint3>>>
384# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
389 std::array<CGALPoint3, 4> pts0;
390 std::array<CGALPoint3, 3> pts1;
392 std::transform(tetrahedron.begin(),
395 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
397 std::transform(triangle.begin(),
400 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
402 CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
404 return convert_boost_to_std(
405 CGAL::intersection(cgal_triangle, cgal_tetrahedron));
411 "This function requires a version of CGAL greater or equal than 5.5."));
419 std::vector<std::array<Point<2>, 3>>
427 const auto intersection_test =
430 if (!intersection_test.empty())
432 const auto &poly = intersection_test[0].outer_boundary();
433 const unsigned int size_poly = poly.size();
439 {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)),
440 CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)),
441 CGALWrappers::cgal_point_to_dealii_point<2>(
444 else if (size_poly >= 4)
447 std::vector<std::array<Point<2>, 3>> collection;
450 cdt.insert_constraint(poly.vertices_begin(),
458 if (f->info().in_domain() &&
459 CGAL::to_double(cdt.triangle(f).area()) > tol)
461 collection.push_back(
462 {{CGALWrappers::cgal_point_to_dealii_point<2>(
463 cdt.triangle(f).vertex(0)),
464 CGALWrappers::cgal_point_to_dealii_point<2>(
465 cdt.triangle(f).vertex(1)),
466 CGALWrappers::cgal_point_to_dealii_point<2>(
467 cdt.triangle(f).vertex(2))}});
485 std::vector<std::array<Point<2>, 2>>
493 std::array<CGALPoint2, 4> pts;
495 std::transform(quad.begin(),
498 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
503 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[0]),
504 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[1]));
506 cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(),
true);
507 std::vector<std::array<Point<2>, 2>> vertices;
511 if (f->info().in_domain() &&
512 CGAL::to_double(cdt.triangle(f).area()) > tol &&
513 CGAL::do_intersect(segm, cdt.triangle(f)))
515 const auto intersection =
516 CGAL::intersection(segm, cdt.triangle(f));
517 if (
const CGALSegment2 *s = get_if_<CGALSegment2>(&*intersection))
520 {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]),
521 CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}});
530 std::vector<std::array<Point<3>, 2>>
535# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
540 std::array<CGALPoint3_exact, 8> pts;
546 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
549 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[0]),
550 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[1]));
554 std::vector<std::array<Point<3>, 2>> vertices;
556 cgal_triangulation.insert(pts.begin(), pts.end());
557 for (
const auto &c : cgal_triangulation.finite_cell_handles())
559 const auto &cgal_tetrahedron = cgal_triangulation.tetrahedron(c);
560 if (CGAL::do_intersect(cgal_segment, cgal_tetrahedron))
562 const auto intersection =
563 CGAL::intersection(cgal_segment, cgal_tetrahedron);
565 get_if_<CGALSegment3_exact>(&*intersection))
567 if (s->squared_length() > tol * tol)
570 {{CGALWrappers::cgal_point_to_dealii_point<3>(
572 CGALWrappers::cgal_point_to_dealii_point<3>(
583 "This function requires a version of CGAL greater or equal than 5.5."));
591 std::vector<std::array<Point<3>, 3>>
596# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
601 std::array<CGALPoint3_exact, 8> pts_hex;
602 std::array<CGALPoint3_exact, 4> pts_quad;
608 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
614 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
617 std::vector<std::array<Point<3>, 3>> vertices;
619 triangulation_hexa.insert(pts_hex.begin(), pts_hex.end());
623 triangulation_quad.insert(pts_quad.begin(), pts_quad.end());
625 for (
const auto &c : triangulation_hexa.finite_cell_handles())
627 const auto &tet = triangulation_hexa.tetrahedron(c);
629 for (
const auto &f : triangulation_quad.finite_facets())
631 if (CGAL::do_intersect(tet, triangulation_quad.triangle(f)))
633 const auto intersection =
634 CGAL::intersection(triangulation_quad.triangle(f), tet);
637 get_if_<CGALTriangle3_exact>(&*intersection))
639 if (CGAL::to_double(t->squared_area()) > tol * tol)
642 {{cgal_point_to_dealii_point<3>((*t)[0]),
643 cgal_point_to_dealii_point<3>((*t)[1]),
644 cgal_point_to_dealii_point<3>((*t)[2])}});
648 if (
const std::vector<CGALPoint3_exact> *vps =
649 get_if_<std::vector<CGALPoint3_exact>>(&*intersection))
652 tria_inter.insert(vps->begin(), vps->end());
654 for (
auto it = tria_inter.finite_facets_begin();
655 it != tria_inter.finite_facets_end();
658 const auto triangle = tria_inter.triangle(*it);
659 if (CGAL::to_double(triangle.squared_area()) >
662 std::array<Point<3>, 3> verts = {
663 {CGALWrappers::cgal_point_to_dealii_point<3>(
665 CGALWrappers::cgal_point_to_dealii_point<3>(
667 CGALWrappers::cgal_point_to_dealii_point<3>(
670 vertices.push_back(verts);
683 "This function requires a version of CGAL greater or equal than 5.5."));
691 std::vector<std::array<Point<3>, 4>>
699 std::array<CGALPoint3_exact, 8> pts_hex0;
700 std::array<CGALPoint3_exact, 8> pts_hex1;
706 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
712 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
716 std::vector<std::array<Point<3>, 4>> vertices;
719 tria0.insert(pts_hex0.begin(), pts_hex0.end());
720 tria1.insert(pts_hex1.begin(), pts_hex1.end());
722 for (
const auto &c0 : tria0.finite_cell_handles())
724 const auto &tet0 = tria1.tetrahedron(c0);
725 const auto &tetg0 = CGAL::make_tetrahedron(tet0.vertex(0),
731 for (
const auto &c1 : tria1.finite_cell_handles())
733 const auto &tet1 = tria1.tetrahedron(c1);
734 const auto &tetg1 = CGAL::make_tetrahedron(tet1.vertex(0),
740 namespace PMP = CGAL::Polygon_mesh_processing;
741 const bool test_intersection =
742 PMP::corefine_and_compute_intersection(surf0, surf1, sm);
743 if (PMP::volume(sm) > tol && test_intersection)
747 triangulation_hexa.insert(sm.points().begin(),
749 for (
const auto &c : triangulation_hexa.finite_cell_handles())
751 const auto &tet = triangulation_hexa.tetrahedron(c);
753 {{CGALWrappers::cgal_point_to_dealii_point<3>(
755 CGALWrappers::cgal_point_to_dealii_point<3>(
757 CGALWrappers::cgal_point_to_dealii_point<3>(
759 CGALWrappers::cgal_point_to_dealii_point<3>(
774 template <
int structdim0,
int structdim1,
int spacedim>
775 std::vector<std::array<Point<spacedim>, structdim1 + 1>>
781 const unsigned int n_vertices0 = vertices0.size();
782 const unsigned int n_vertices1 = vertices1.size();
785 n_vertices0 > 0 || n_vertices1 > 0,
787 "The intersection cannot be computed as at least one of the two cells has no vertices."));
789 if constexpr (structdim0 == 2 && structdim1 == 2 && spacedim == 2)
791 if (n_vertices0 == 4 && n_vertices1 == 4)
798 else if constexpr (structdim0 == 2 && structdim1 == 1 && spacedim == 2)
800 if (n_vertices0 == 4 && n_vertices1 == 2)
807 else if constexpr (structdim0 == 3 && structdim1 == 1 && spacedim == 3)
809 if (n_vertices0 == 8 && n_vertices1 == 2)
816 else if constexpr (structdim0 == 3 && structdim1 == 2 && spacedim == 3)
818 if (n_vertices0 == 8 && n_vertices1 == 4)
825 else if constexpr (structdim0 == 3 && structdim1 == 3 && spacedim == 3)
827 if (n_vertices0 == 8 && n_vertices1 == 8)
844 template <
int structdim0,
int structdim1,
int spacedim>
845 std::vector<std::array<Point<spacedim>, structdim1 + 1>>
854 ReferenceCells::get_hypercube<structdim0>().n_vertices(),
857 ReferenceCells::get_hypercube<structdim1>().n_vertices(),
860 const auto &vertices0 =
861 CGALWrappers::get_vertices_in_cgal_order(cell0, mapping0);
862 const auto &vertices1 =
863 CGALWrappers::get_vertices_in_cgal_order(cell1, mapping1);
865 return compute_intersection_of_cells<structdim0, structdim1, spacedim>(
866 vertices0, vertices1, tol);
869# include "cgal/intersections.inst"
879template <
int structdim0,
884std::vector<std::array<Point<spacedim>, structdim1 + 1>>
896template <
int structdim0,
int structdim1,
int spacedim>
897std::vector<std::array<Point<spacedim>, structdim1 + 1>>
Abstract base class for mapping classes.
virtual boost::container::small_vector< Point< spacedim >, ReferenceCells::max_n_vertices< dim >() > 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
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcNeedsCGAL()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< std::array< Point< 2 >, 2 > > compute_intersection_quad_line(const ArrayView< const Point< 2 > > &quad, const ArrayView< const Point< 2 > > &line, const double tol)
std::vector< std::array< Point< 3 >, 3 > > compute_intersection_hexa_quad(const ArrayView< const Point< 3 > > &hexa, const ArrayView< const Point< 3 > > &quad, const double tol)
void mark_domains(CDT &ct, Face_handle start, int index, std::list< CDT::Edge > &border)
std::optional< std::variant< CGALPoint3, CGALSegment3 > > compute_intersection_tetra_segment(const ArrayView< const Point< 3 > > &tetrahedron, const ArrayView< const Point< 3 > > &segment)
std::vector< Polygon_with_holes_2 > compute_intersection_rect_rect(const ArrayView< const Point< 2 > > &rectangle0, const ArrayView< const Point< 2 > > &rectangle1)
std::vector< std::array< Point< 3 >, 4 > > compute_intersection_hexa_hexa(const ArrayView< const Point< 3 > > &hexa0, const ArrayView< const Point< 3 > > &hexa1, const double tol)
std::vector< std::array< Point< 3 >, 2 > > compute_intersection_hexa_line(const ArrayView< const Point< 3 > > &hexa, const ArrayView< const Point< 3 > > &line, const double tol)
std::optional< std::variant< CGALPoint3, CGALSegment3, CGALTriangle3, std::vector< CGALPoint3 > > > compute_intersection_tetra_triangle(const ArrayView< const Point< 3 > > &tetrahedron, const ArrayView< const Point< 3 > > &triangle)
std::optional< std::variant< CGALPoint2, CGALSegment2 > > compute_intersection_triangle_segment(const ArrayView< const Point< 2 > > &triangle, const ArrayView< const Point< 2 > > &segment)
std::optional< std::variant< CGALPoint2, CGALSegment2, CGALTriangle2, std::vector< CGALPoint2 > > > compute_intersection_triangle_triangle(const ArrayView< const Point< 2 > > &triangle0, const ArrayView< const Point< 2 > > &triangle1)
std::vector< std::array< Point< 2 >, 3 > > compute_intersection_quad_quad(const ArrayView< const Point< 2 > > &quad0, const ArrayView< const Point< 2 > > &quad1, const double tol)
K::Tetrahedron_3 CGALTetra
CGAL::Surface_mesh< K_exact::Point_3 > Surface_mesh
K::Segment_2 CGALSegment2
const T * get_if_(const std::variant< Types... > *v)
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_exact > Triangulation3_exact
CGAL::Constrained_Delaunay_triangulation_2< K, Tds, Itag > CDT
std::vector< std::array< Point< spacedim >, structdim1+1 > > compute_intersection_of_cells(const typename Triangulation< structdim0, spacedim >::cell_iterator &cell0, const typename Triangulation< structdim1, spacedim >::cell_iterator &cell1, const Mapping< structdim0, spacedim > &mapping0, const Mapping< structdim1, spacedim > &mapping1, const double tol=1e-9)
CGAL::Delaunay_mesh_size_criteria_2< CDT > Criteria
CGAL::Triangulation_vertex_base_2< K > Vb
CGAL::Polygon_with_holes_2< K > Polygon_with_holes_2
CGAL::Triangulation_data_structure_2< Vb, Fb > Tds
CGAL::Polygon_2< K > CGALPolygon
CGAL::Triangulation_2< K > Triangulation2
CGAL::Exact_predicates_exact_constructions_kernel K_exact
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
K_exact::Point_3 CGALPoint3_exact
CGAL::Triangulation_face_base_with_info_2< FaceInfo2, K > Fbb