22 #ifdef DEAL_II_WITH_CGAL
33 # 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/Exact_predicates_inexact_constructions_kernel.h>
43 # include <CGAL/Kernel_traits.h>
44 # include <CGAL/Polygon_2.h>
45 # include <CGAL/Polygon_with_holes_2.h>
46 # include <CGAL/Projection_traits_xy_3.h>
47 # include <CGAL/Segment_3.h>
48 # include <CGAL/Simple_cartesian.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>
59 # include <type_traits>
65 using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt;
66 using K_exact = CGAL::Exact_predicates_exact_constructions_kernel;
67 using K_inexact = CGAL::Exact_predicates_inexact_constructions_kernel;
100 using Vb = CGAL::Triangulation_vertex_base_2<K>;
101 using Fbb = CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>;
102 using CFb = CGAL::Constrained_triangulation_face_base_2<K, Fbb>;
103 using Fb = CGAL::Delaunay_mesh_face_base_2<K, CFb>;
104 using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
105 using Itag = CGAL::Exact_predicates_tag;
106 using CDT = CGAL::Constrained_Delaunay_triangulation_2<K, Tds, Itag>;
107 using Criteria = CGAL::Delaunay_mesh_size_criteria_2<CDT>;
119 template <
typename TargetVariant>
120 struct Repackage : boost::static_visitor<TargetVariant>
122 template <
typename T>
124 operator()(
const T &t)
const
126 return TargetVariant(t);
136 template <
typename... Types>
137 std_cxx17::optional<std_cxx17::variant<Types...>>
138 convert_boost_to_std(
const boost::optional<boost::variant<Types...>> &x)
147 using std_variant = std::variant<Types...>;
148 return boost::apply_visitor(Repackage<std_variant>(), *x);
165 std::list<CDT::Edge> &border)
167 if (start->info().nesting_level != -1)
171 std::list<Face_handle> queue;
172 queue.push_back(start);
173 while (!queue.empty())
177 if (fh->info().nesting_level == -1)
179 fh->info().nesting_level = index;
180 for (
int i = 0; i < 3; i++)
184 if (n->info().nesting_level == -1)
186 if (ct.is_constrained(
e))
203 f->info().nesting_level = -1;
205 std::list<CDT::Edge> border;
207 while (!border.empty())
209 CDT::Edge
e = border.front();
212 if (n->info().nesting_level == -1)
214 mark_domains(cdt, n,
e.first->info().nesting_level + 1, border);
228 std_cxx17::optional<std_cxx17::variant<
CGALPoint2,
231 std::vector<CGALPoint2>>>
239 std::array<CGALPoint2, 3> pts0, pts1;
244 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
249 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
253 return convert_boost_to_std(
254 CGAL::intersection(cgal_triangle0, cgal_triangle1));
258 std_cxx17::optional<std_cxx17::variant<CGALPoint2, CGALSegment2>>
266 std::array<CGALPoint2, 3> pts0;
267 std::array<CGALPoint2, 2> pts1;
272 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
277 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
281 return convert_boost_to_std(
282 CGAL::intersection(cgal_segment, cgal_triangle));
288 std::vector<Polygon_with_holes_2>
295 std::array<CGALPoint2, 4> pts0, pts1;
300 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
305 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
307 const CGALPolygon first_poly{pts0.begin(), pts0.end()};
308 const CGALPolygon second_poly{pts1.begin(), pts1.end()};
310 std::vector<Polygon_with_holes_2> poly_list;
311 CGAL::intersection(first_poly,
313 std::back_inserter(poly_list));
319 std_cxx17::optional<std_cxx17::variant<CGALPoint3, CGALSegment3>>
324 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
329 std::array<CGALPoint3, 4> pts0;
330 std::array<CGALPoint3, 2> pts1;
335 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
340 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
342 CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
344 return convert_boost_to_std(
345 CGAL::intersection(cgal_segment, cgal_tetrahedron));
350 "This function requires a version of CGAL greater or equal than 5.5."));
359 std_cxx17::optional<std_cxx17::variant<
CGALPoint3,
362 std::vector<CGALPoint3>>>
367 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
372 std::array<CGALPoint3, 4> pts0;
373 std::array<CGALPoint3, 3> pts1;
378 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
383 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
385 CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
387 return convert_boost_to_std(
388 CGAL::intersection(cgal_triangle, cgal_tetrahedron));
394 "This function requires a version of CGAL greater or equal than 5.5."));
402 std::vector<std::array<Point<2>, 3>>
410 const auto intersection_test =
413 if (!intersection_test.empty())
415 const auto & poly = intersection_test[0].outer_boundary();
416 const unsigned int size_poly = poly.size();
422 {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)),
423 CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)),
424 CGALWrappers::cgal_point_to_dealii_point<2>(
427 else if (size_poly >= 4)
430 std::vector<std::array<Point<2>, 3>> collection;
433 cdt.insert_constraint(poly.vertices_begin(),
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))}});
469 std::vector<std::array<Point<2>, 2>>
477 std::array<CGALPoint2, 4> pts;
482 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
487 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[0]),
488 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[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 =
500 CGAL::intersection(segm, cdt.triangle(f));
502 boost::get<CGALSegment2>(&*intersection))
505 {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]),
506 CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}});
515 std::vector<std::array<Point<3>, 2>>
520 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
525 std::array<CGALPoint3_exact, 8> pts;
531 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
534 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[0]),
535 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[1]));
539 std::vector<std::array<Point<3>, 2>>
vertices;
541 cgal_triangulation.insert(pts.begin(), pts.end());
542 for (
const auto &c : cgal_triangulation.finite_cell_handles())
544 const auto &cgal_tetrahedron = cgal_triangulation.tetrahedron(c);
545 if (CGAL::do_intersect(cgal_segment, cgal_tetrahedron))
547 const auto intersection =
548 CGAL::intersection(cgal_segment, cgal_tetrahedron);
550 boost::get<CGALSegment3_exact>(&*intersection))
552 if (s->squared_length() > tol * tol)
555 {{CGALWrappers::cgal_point_to_dealii_point<3>(
557 CGALWrappers::cgal_point_to_dealii_point<3>(
568 "This function requires a version of CGAL greater or equal than 5.5."));
576 std::vector<std::array<Point<3>, 3>>
581 # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
586 std::array<CGALPoint3_exact, 8> pts_hex;
587 std::array<CGALPoint3_exact, 4> pts_quad;
593 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
599 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
602 std::vector<std::array<Point<3>, 3>>
vertices;
604 triangulation_hexa.insert(pts_hex.begin(), pts_hex.end());
608 triangulation_quad.insert(pts_quad.begin(), pts_quad.end());
610 for (
const auto &c : triangulation_hexa.finite_cell_handles())
612 const auto &tet = triangulation_hexa.tetrahedron(c);
614 for (
const auto &f : triangulation_quad.finite_facets())
616 if (CGAL::do_intersect(tet, triangulation_quad.triangle(f)))
618 const auto intersection =
619 CGAL::intersection(triangulation_quad.triangle(f), tet);
622 boost::get<CGALTriangle3_exact>(&*intersection))
624 if (CGAL::to_double(t->squared_area()) > tol * tol)
627 {{cgal_point_to_dealii_point<3>((*t)[0]),
628 cgal_point_to_dealii_point<3>((*t)[1]),
629 cgal_point_to_dealii_point<3>((*t)[2])}});
633 if (
const std::vector<CGALPoint3_exact> *vps =
634 boost::get<std::vector<CGALPoint3_exact>>(
638 tria_inter.insert(vps->begin(), vps->end());
640 for (
auto it = tria_inter.finite_facets_begin();
641 it != tria_inter.finite_facets_end();
644 const auto triangle = tria_inter.triangle(*it);
645 if (CGAL::to_double(triangle.squared_area()) >
648 std::array<Point<3>, 3> verts = {
649 {CGALWrappers::cgal_point_to_dealii_point<3>(
651 CGALWrappers::cgal_point_to_dealii_point<3>(
653 CGALWrappers::cgal_point_to_dealii_point<3>(
669 "This function requires a version of CGAL greater or equal than 5.5."));
677 std::vector<std::array<Point<3>, 4>>
685 std::array<CGALPoint3_inexact, 8> pts_hex0;
686 std::array<CGALPoint3_inexact, 8> pts_hex1;
692 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact, 3>);
698 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact, 3>);
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 const auto &tetg0 = CGAL::make_tetrahedron(tet0.vertex(0),
717 for (
const auto &c1 : tria1.finite_cell_handles())
719 const auto &tet1 = tria1.tetrahedron(c1);
720 const auto &tetg1 = 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 triangulation_hexa.insert(sm.points().begin(),
735 for (
const auto &c : triangulation_hexa.finite_cell_handles())
737 const auto &tet = triangulation_hexa.tetrahedron(c);
739 {{CGALWrappers::cgal_point_to_dealii_point<3>(
741 CGALWrappers::cgal_point_to_dealii_point<3>(
743 CGALWrappers::cgal_point_to_dealii_point<3>(
745 CGALWrappers::cgal_point_to_dealii_point<3>(
760 template <
int dim0,
int dim1,
int spacedim>
761 std::vector<std::array<Point<spacedim>, dim1 + 1>>
767 const unsigned int n_vertices0 = vertices0.size();
768 const unsigned int n_vertices1 = vertices1.size();
771 n_vertices0 > 0 || n_vertices1 > 0,
773 "The intersection cannot be computed as at least one of the two cells has no vertices."));
775 if constexpr (dim0 == 2 && dim1 == 2 && spacedim == 2)
777 if (n_vertices0 == 4 && n_vertices1 == 4)
784 else if constexpr (dim0 == 2 && dim1 == 1 && spacedim == 2)
786 if (n_vertices0 == 4 && n_vertices1 == 2)
793 else if constexpr (dim0 == 3 && dim1 == 1 && spacedim == 3)
795 if (n_vertices0 == 8 && n_vertices1 == 2)
802 else if constexpr (dim0 == 3 && dim1 == 2 && spacedim == 3)
804 if (n_vertices0 == 8 && n_vertices1 == 4)
811 else if constexpr (dim0 == 3 && dim1 == 3 && spacedim == 3)
813 if (n_vertices0 == 8 && n_vertices1 == 8)
830 template <
int dim0,
int dim1,
int spacedim>
831 std::vector<std::array<Point<spacedim>, dim1 + 1>>
840 ReferenceCells::get_hypercube<dim0>().n_vertices(),
843 ReferenceCells::get_hypercube<dim1>().n_vertices(),
846 const auto &vertices0 =
847 CGALWrappers::get_vertices_in_cgal_order(cell0, mapping0);
848 const auto &vertices1 =
849 CGALWrappers::get_vertices_in_cgal_order(cell1, mapping1);
851 return compute_intersection_of_cells<dim0, dim1, spacedim>(vertices0,
856 # include "intersections.inst"
871 std::vector<std::array<Point<spacedim>, dim1 + 1>>
883 template <
int dim0,
int dim1,
int spacedim>
884 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_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcNeedsCGAL()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std_cxx17::optional< std_cxx17::variant< CGALPoint3, CGALSegment3, CGALTriangle3, std::vector< CGALPoint3 > > > compute_intersection_tetra_triangle(const ArrayView< const Point< 3 >> &tetrahedron, const ArrayView< const Point< 3 >> &triangle)
std_cxx17::optional< std_cxx17::variant< CGALPoint2, CGALSegment2, CGALTriangle2, std::vector< CGALPoint2 > > > compute_intersection_triangle_triangle(const ArrayView< const Point< 2 >> &triangle0, const ArrayView< const Point< 2 >> &triangle1)
void mark_domains(CDT &ct, Face_handle start, int index, std::list< CDT::Edge > &border)
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< 2 >, 3 > > compute_intersection_quad_quad(const ArrayView< const Point< 2 >> &quad0, const ArrayView< const Point< 2 >> &quad1, const double tol)
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 >, 2 > > compute_intersection_hexa_line(const ArrayView< const Point< 3 >> &hexa, const ArrayView< const Point< 3 >> &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)
std_cxx17::optional< std_cxx17::variant< CGALPoint3, CGALSegment3 > > compute_intersection_tetra_segment(const ArrayView< const Point< 3 >> &tetrahedron, const ArrayView< const Point< 3 >> &segment)
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_cxx17::optional< std_cxx17::variant< CGALPoint2, CGALSegment2 > > compute_intersection_triangle_segment(const ArrayView< const Point< 2 >> &triangle, const ArrayView< const Point< 2 >> &segment)
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
CGAL::Polygon_with_holes_2< K > Polygon_with_holes_2
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
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
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)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Function &function, const unsigned int grainsize)