Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
intersections.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2022 - 2023 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/config.h>
17
19
20#include <algorithm>
21
22#ifdef DEAL_II_WITH_CGAL
23
28
29# include <deal.II/fe/mapping.h>
30
31# include <deal.II/grid/tria.h>
32
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/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/Tetrahedron_3.h>
49# include <CGAL/Triangle_2.h>
50# include <CGAL/Triangle_3.h>
51# include <CGAL/Triangulation_2.h>
52# include <CGAL/Triangulation_3.h>
53# include <CGAL/Triangulation_face_base_with_id_2.h>
54# include <CGAL/Triangulation_face_base_with_info_2.h>
56
57# include <fstream>
58# include <type_traits>
59
61
62namespace CGALWrappers
63{
64 using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt;
65 using K_exact = CGAL::Exact_predicates_exact_constructions_kernel;
66 using CGALPolygon = CGAL::Polygon_2<K>;
67 using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2<K>;
68 using CGALTriangle2 = K::Triangle_2;
69 using CGALTriangle3 = K::Triangle_3;
70 using CGALTriangle3_exact = K_exact::Triangle_3;
71 using CGALPoint2 = K::Point_2;
72 using CGALPoint3 = K::Point_3;
73 using CGALPoint3_exact = K_exact::Point_3;
74 using CGALSegment2 = K::Segment_2;
75 using Surface_mesh = CGAL::Surface_mesh<K_exact::Point_3>;
76 using CGALSegment3 = K::Segment_3;
77 using CGALSegment3_exact = K_exact::Segment_3;
78 using CGALTetra = K::Tetrahedron_3;
79 using CGALTetra_exact = K_exact::Tetrahedron_3;
80 using Triangulation2 = CGAL::Triangulation_2<K>;
81 using Triangulation3 = CGAL::Triangulation_3<K>;
82 using Triangulation3_exact = CGAL::Triangulation_3<K_exact>;
83
84 struct FaceInfo2
85 {
87 {}
89 bool
91 {
92 return nesting_level % 2 == 1;
93 }
94 };
95
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>;
104 using Vertex_handle = CDT::Vertex_handle;
105 using Face_handle = CDT::Face_handle;
106
107 namespace internal
108 {
109 namespace
110 {
115 template <typename TargetVariant>
116 struct Repackage : boost::static_visitor<TargetVariant>
117 {
118 template <typename T>
119 TargetVariant
120 operator()(const T &t) const
121 {
122 return TargetVariant(t);
123 }
124 };
125
132 template <typename... Types>
133 std_cxx17::optional<std_cxx17::variant<Types...>>
134 convert_boost_to_std(const boost::optional<boost::variant<Types...>> &x)
135 {
136 if (x)
137 {
138 // The boost::optional object contains an object of type
139 // boost::variant. We need to unpack which type the
140 // variant contains, and re-package that into a
141 // std::variant. This is easily done using a visitor
142 // object.
143 using std_variant = std::variant<Types...>;
144 return boost::apply_visitor(Repackage<std_variant>(), *x);
145 }
146 else
147 {
148 // The boost::optional object was empty. Return an empty
149 // std::optional object.
150 return {};
151 }
152 }
153 } // namespace
154
155
156
157 void
159 Face_handle start,
160 int index,
161 std::list<CDT::Edge> &border)
162 {
163 if (start->info().nesting_level != -1)
164 {
165 return;
166 }
167 std::list<Face_handle> queue;
168 queue.push_back(start);
169 while (!queue.empty())
170 {
171 Face_handle fh = queue.front();
172 queue.pop_front();
173 if (fh->info().nesting_level == -1)
174 {
175 fh->info().nesting_level = index;
176 for (int i = 0; i < 3; i++)
177 {
178 CDT::Edge e(fh, i);
179 Face_handle n = fh->neighbor(i);
180 if (n->info().nesting_level == -1)
181 {
182 if (ct.is_constrained(e))
183 border.push_back(e);
184 else
185 queue.push_back(n);
186 }
187 }
188 }
189 }
190 }
191
192
193
194 void
196 {
197 for (CDT::Face_handle f : cdt.all_face_handles())
198 {
199 f->info().nesting_level = -1;
200 }
201 std::list<CDT::Edge> border;
202 mark_domains(cdt, cdt.infinite_face(), 0, border);
203 while (!border.empty())
204 {
205 CDT::Edge e = border.front();
206 border.pop_front();
207 Face_handle n = e.first->neighbor(e.second);
208 if (n->info().nesting_level == -1)
209 {
210 mark_domains(cdt, n, e.first->info().nesting_level + 1, border);
211 }
212 }
213 }
214
215 // Collection of utilities that compute intersection between simplices
216 // identified by array of points. The return type is the one of
217 // CGAL::intersection(), i.e. a std_cxx17::optional<std_cxx17::variant<>>.
218 // Intersection between 2d and 3d objects and 1d/3d objects are available
219 // only with CGAL versions greater or equal than 5.5, hence the
220 // corresponding functions are guarded by #ifdef directives. All the
221 // signatures follow the convection that the first entity has an intrinsic
222 // dimension higher than the second one.
223
224 std_cxx17::optional<std_cxx17::variant<CGALPoint2,
227 std::vector<CGALPoint2>>>
229 const ArrayView<const Point<2>> &triangle0,
230 const ArrayView<const Point<2>> &triangle1)
231 {
232 AssertDimension(triangle0.size(), 3);
233 AssertDimension(triangle0.size(), triangle1.size());
234
235 std::array<CGALPoint2, 3> pts0, pts1;
236
237 std::transform(triangle0.begin(),
238 triangle0.end(),
239 pts0.begin(),
240 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
241
242 std::transform(triangle1.begin(),
243 triangle1.end(),
244 pts1.begin(),
245 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
246
247 CGALTriangle2 cgal_triangle0{pts0[0], pts0[1], pts0[2]};
248 CGALTriangle2 cgal_triangle1{pts1[0], pts1[1], pts1[2]};
249 return convert_boost_to_std(
250 CGAL::intersection(cgal_triangle0, cgal_triangle1));
251 }
252
253
254 std_cxx17::optional<std_cxx17::variant<CGALPoint2, CGALSegment2>>
256 const ArrayView<const Point<2>> &triangle,
257 const ArrayView<const Point<2>> &segment)
258 {
259 AssertDimension(triangle.size(), 3);
260 AssertDimension(segment.size(), 2);
261
262 std::array<CGALPoint2, 3> pts0;
263 std::array<CGALPoint2, 2> pts1;
264
265 std::transform(triangle.begin(),
266 triangle.end(),
267 pts0.begin(),
268 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
269
270 std::transform(segment.begin(),
271 segment.end(),
272 pts1.begin(),
273 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
274
275 CGALTriangle2 cgal_triangle{pts0[0], pts0[1], pts0[2]};
276 CGALSegment2 cgal_segment{pts1[0], pts1[1]};
277 return convert_boost_to_std(
278 CGAL::intersection(cgal_segment, cgal_triangle));
279 }
280
281
282
283 // rectangle-rectangle
284 std::vector<Polygon_with_holes_2>
286 const ArrayView<const Point<2>> &rectangle1)
287 {
288 AssertDimension(rectangle0.size(), 4);
289 AssertDimension(rectangle0.size(), rectangle1.size());
290
291 std::array<CGALPoint2, 4> pts0, pts1;
292
293 std::transform(rectangle0.begin(),
294 rectangle0.end(),
295 pts0.begin(),
296 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
297
298 std::transform(rectangle1.begin(),
299 rectangle1.end(),
300 pts1.begin(),
301 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
302
303 const CGALPolygon first_poly{pts0.begin(), pts0.end()};
304 const CGALPolygon second_poly{pts1.begin(), pts1.end()};
305
306 std::vector<Polygon_with_holes_2> poly_list;
307 CGAL::intersection(first_poly,
308 second_poly,
309 std::back_inserter(poly_list));
310 return poly_list;
311 }
312
313
314
315 std_cxx17::optional<std_cxx17::variant<CGALPoint3, CGALSegment3>>
317 const ArrayView<const Point<3>> &tetrahedron,
318 const ArrayView<const Point<3>> &segment)
319 {
320# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
321
322 AssertDimension(tetrahedron.size(), 4);
323 AssertDimension(segment.size(), 2);
324
325 std::array<CGALPoint3, 4> pts0;
326 std::array<CGALPoint3, 2> pts1;
327
328 std::transform(tetrahedron.begin(),
329 tetrahedron.end(),
330 pts0.begin(),
331 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
332
333 std::transform(segment.begin(),
334 segment.end(),
335 pts1.begin(),
336 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
337
338 CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
339 CGALSegment3 cgal_segment{pts1[0], pts1[1]};
340 return convert_boost_to_std(
341 CGAL::intersection(cgal_segment, cgal_tetrahedron));
342# else
343 Assert(
344 false,
346 "This function requires a version of CGAL greater or equal than 5.5."));
347 (void)tetrahedron;
348 (void)segment;
349 return {};
350# endif
351 }
352
353
354 // tetra, triangle
355 std_cxx17::optional<std_cxx17::variant<CGALPoint3,
358 std::vector<CGALPoint3>>>
360 const ArrayView<const Point<3>> &tetrahedron,
361 const ArrayView<const Point<3>> &triangle)
362 {
363# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
364
365 AssertDimension(tetrahedron.size(), 4);
366 AssertDimension(triangle.size(), 3);
367
368 std::array<CGALPoint3, 4> pts0;
369 std::array<CGALPoint3, 3> pts1;
370
371 std::transform(tetrahedron.begin(),
372 tetrahedron.end(),
373 pts0.begin(),
374 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
375
376 std::transform(triangle.begin(),
377 triangle.end(),
378 pts1.begin(),
379 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
380
381 CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
382 CGALTriangle3 cgal_triangle{pts1[0], pts1[1], pts1[2]};
383 return convert_boost_to_std(
384 CGAL::intersection(cgal_triangle, cgal_tetrahedron));
385# else
386
387 Assert(
388 false,
390 "This function requires a version of CGAL greater or equal than 5.5."));
391 (void)tetrahedron;
392 (void)triangle;
393 return {};
394# endif
395 }
396
397 // quad-quad
398 std::vector<std::array<Point<2>, 3>>
400 const ArrayView<const Point<2>> &quad1,
401 const double tol)
402 {
403 AssertDimension(quad0.size(), 4);
404 AssertDimension(quad0.size(), quad1.size());
405
406 const auto intersection_test =
408
409 if (!intersection_test.empty())
410 {
411 const auto & poly = intersection_test[0].outer_boundary();
412 const unsigned int size_poly = poly.size();
413 if (size_poly == 3)
414 {
415 // intersection is a triangle itself, so directly return its
416 // vertices.
417 return {
418 {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)),
419 CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)),
420 CGALWrappers::cgal_point_to_dealii_point<2>(
421 poly.vertex(2))}}};
422 }
423 else if (size_poly >= 4)
424 {
425 // intersection is a polygon, need to triangulate it.
426 std::vector<std::array<Point<2>, 3>> collection;
427
428 CDT cdt;
429 cdt.insert_constraint(poly.vertices_begin(),
430 poly.vertices_end(),
431 true);
432
434 std::array<Point<2>, 3> vertices;
435
436 for (Face_handle f : cdt.finite_face_handles())
437 {
438 if (f->info().in_domain() &&
439 CGAL::to_double(cdt.triangle(f).area()) > tol)
440 {
441 collection.push_back(
442 {{CGALWrappers::cgal_point_to_dealii_point<2>(
443 cdt.triangle(f).vertex(0)),
444 CGALWrappers::cgal_point_to_dealii_point<2>(
445 cdt.triangle(f).vertex(1)),
446 CGALWrappers::cgal_point_to_dealii_point<2>(
447 cdt.triangle(f).vertex(2))}});
448 }
449 }
450 return collection;
451 }
452 else
453 {
454 Assert(false, ExcMessage("The polygon is degenerate."));
455 return {};
456 }
457 }
458 else
459 {
460 return {};
461 }
462 }
463
464 // Specialization for quad \cap line
465 std::vector<std::array<Point<2>, 2>>
467 const ArrayView<const Point<2>> &line,
468 const double tol)
469 {
470 AssertDimension(quad.size(), 4);
471 AssertDimension(line.size(), 2);
472
473 std::array<CGALPoint2, 4> pts;
474
475 std::transform(quad.begin(),
476 quad.end(),
477 pts.begin(),
478 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
479
480 CGALPolygon poly(pts.begin(), pts.end());
481
482 CGALSegment2 segm(
483 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[0]),
484 CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[1]));
485 CDT cdt;
486 cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true);
487 std::vector<std::array<Point<2>, 2>> vertices;
489 for (Face_handle f : cdt.finite_face_handles())
490 {
491 if (f->info().in_domain() &&
492 CGAL::to_double(cdt.triangle(f).area()) > tol &&
493 CGAL::do_intersect(segm, cdt.triangle(f)))
494 {
495 const auto intersection =
496 CGAL::intersection(segm, cdt.triangle(f));
497 if (const CGALSegment2 *s =
498 boost::get<CGALSegment2>(&*intersection))
499 {
500 vertices.push_back(
501 {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]),
502 CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}});
503 }
504 }
505 }
506
507 return vertices;
508 }
509
510 // specialization for hex \cap line
511 std::vector<std::array<Point<3>, 2>>
513 const ArrayView<const Point<3>> &line,
514 const double tol)
515 {
516# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
517
518 AssertDimension(hexa.size(), 8);
519 AssertDimension(line.size(), 2);
520
521 std::array<CGALPoint3_exact, 8> pts;
522
523 std::transform(
524 hexa.begin(),
525 hexa.end(),
526 pts.begin(),
527 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
528
529 CGALSegment3_exact cgal_segment(
530 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[0]),
531 CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[1]));
532
533 // Subdivide the hex into tetrahedrons, and intersect each one of them
534 // with the line
535 std::vector<std::array<Point<3>, 2>> vertices;
536 Triangulation3_exact cgal_triangulation;
537 cgal_triangulation.insert(pts.begin(), pts.end());
538 for (const auto &c : cgal_triangulation.finite_cell_handles())
539 {
540 const auto &cgal_tetrahedron = cgal_triangulation.tetrahedron(c);
541 if (CGAL::do_intersect(cgal_segment, cgal_tetrahedron))
542 {
543 const auto intersection =
544 CGAL::intersection(cgal_segment, cgal_tetrahedron);
545 if (const CGALSegment3_exact *s =
546 boost::get<CGALSegment3_exact>(&*intersection))
547 {
548 if (s->squared_length() > tol * tol)
549 {
550 vertices.push_back(
551 {{CGALWrappers::cgal_point_to_dealii_point<3>(
552 s->vertex(0)),
553 CGALWrappers::cgal_point_to_dealii_point<3>(
554 s->vertex(1))}});
555 }
556 }
557 }
558 }
559 return vertices;
560# else
561 Assert(
562 false,
564 "This function requires a version of CGAL greater or equal than 5.5."));
565 (void)hexa;
566 (void)line;
567 (void)tol;
568 return {};
569# endif
570 }
571
572 std::vector<std::array<Point<3>, 3>>
574 const ArrayView<const Point<3>> &quad,
575 const double tol)
576 {
577# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
578
579 AssertDimension(hexa.size(), 8);
580 AssertDimension(quad.size(), 4);
581
582 std::array<CGALPoint3_exact, 8> pts_hex;
583 std::array<CGALPoint3_exact, 4> pts_quad;
584
585 std::transform(
586 hexa.begin(),
587 hexa.end(),
588 pts_hex.begin(),
589 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
590
591 std::transform(
592 quad.begin(),
593 quad.end(),
594 pts_quad.begin(),
595 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
596
597 // Subdivide hex into tetrahedrons
598 std::vector<std::array<Point<3>, 3>> vertices;
599 Triangulation3_exact triangulation_hexa;
600 triangulation_hexa.insert(pts_hex.begin(), pts_hex.end());
601
602 // Subdivide quad into triangles
603 Triangulation3_exact triangulation_quad;
604 triangulation_quad.insert(pts_quad.begin(), pts_quad.end());
605
606 for (const auto &c : triangulation_hexa.finite_cell_handles())
607 {
608 const auto &tet = triangulation_hexa.tetrahedron(c);
609
610 for (const auto &f : triangulation_quad.finite_facets())
611 {
612 if (CGAL::do_intersect(tet, triangulation_quad.triangle(f)))
613 {
614 const auto intersection =
615 CGAL::intersection(triangulation_quad.triangle(f), tet);
616
617 if (const CGALTriangle3_exact *t =
618 boost::get<CGALTriangle3_exact>(&*intersection))
619 {
620 if (CGAL::to_double(t->squared_area()) > tol * tol)
621 {
622 vertices.push_back(
623 {{cgal_point_to_dealii_point<3>((*t)[0]),
624 cgal_point_to_dealii_point<3>((*t)[1]),
625 cgal_point_to_dealii_point<3>((*t)[2])}});
626 }
627 }
628
629 if (const std::vector<CGALPoint3_exact> *vps =
630 boost::get<std::vector<CGALPoint3_exact>>(
631 &*intersection))
632 {
633 Triangulation3_exact tria_inter;
634 tria_inter.insert(vps->begin(), vps->end());
635
636 for (auto it = tria_inter.finite_facets_begin();
637 it != tria_inter.finite_facets_end();
638 ++it)
639 {
640 const auto triangle = tria_inter.triangle(*it);
641 if (CGAL::to_double(triangle.squared_area()) >
642 tol * tol)
643 {
644 std::array<Point<3>, 3> verts = {
645 {CGALWrappers::cgal_point_to_dealii_point<3>(
646 triangle[0]),
647 CGALWrappers::cgal_point_to_dealii_point<3>(
648 triangle[1]),
649 CGALWrappers::cgal_point_to_dealii_point<3>(
650 triangle[2])}};
651
652 vertices.push_back(verts);
653 }
654 }
655 }
656 }
657 }
658 }
659
660 return vertices;
661# else
662 Assert(
663 false,
665 "This function requires a version of CGAL greater or equal than 5.5."));
666 (void)hexa;
667 (void)quad;
668 (void)tol;
669 return {};
670# endif
671 }
672
673 std::vector<std::array<Point<3>, 4>>
675 const ArrayView<const Point<3>> &hexa1,
676 const double tol)
677 {
678 AssertDimension(hexa0.size(), 8);
679 AssertDimension(hexa0.size(), hexa1.size());
680
681 std::array<CGALPoint3_exact, 8> pts_hex0;
682 std::array<CGALPoint3_exact, 8> pts_hex1;
683
684 std::transform(
685 hexa0.begin(),
686 hexa0.end(),
687 pts_hex0.begin(),
688 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
689
690 std::transform(
691 hexa1.begin(),
692 hexa1.end(),
693 pts_hex1.begin(),
694 &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
695
696 Surface_mesh surf0, surf1, sm;
697 // Subdivide hex into tetrahedrons
698 std::vector<std::array<Point<3>, 4>> vertices;
699 Triangulation3_exact tria0, tria1;
700
701 tria0.insert(pts_hex0.begin(), pts_hex0.end());
702 tria1.insert(pts_hex1.begin(), pts_hex1.end());
703
704 for (const auto &c0 : tria0.finite_cell_handles())
705 {
706 const auto &tet0 = tria1.tetrahedron(c0);
707 const auto &tetg0 = CGAL::make_tetrahedron(tet0.vertex(0),
708 tet0.vertex(1),
709 tet0.vertex(2),
710 tet0.vertex(3),
711 surf0);
712 (void)tetg0; // instead of C++ 17s [[maybe unused]]
713 for (const auto &c1 : tria1.finite_cell_handles())
714 {
715 const auto &tet1 = tria1.tetrahedron(c1);
716 const auto &tetg1 = CGAL::make_tetrahedron(tet1.vertex(0),
717 tet1.vertex(1),
718 tet1.vertex(2),
719 tet1.vertex(3),
720 surf1);
721 (void)tetg1; // instead of C++ 17s [[maybe unused]]
722 namespace PMP = CGAL::Polygon_mesh_processing;
723 const bool test_intersection =
724 PMP::corefine_and_compute_intersection(surf0, surf1, sm);
725 if (PMP::volume(sm) > tol && test_intersection)
726 {
727 // Collect tetrahedrons
728 Triangulation3_exact triangulation_hexa;
729 triangulation_hexa.insert(sm.points().begin(),
730 sm.points().end());
731 for (const auto &c : triangulation_hexa.finite_cell_handles())
732 {
733 const auto &tet = triangulation_hexa.tetrahedron(c);
734 vertices.push_back(
735 {{CGALWrappers::cgal_point_to_dealii_point<3>(
736 tet.vertex(0)),
737 CGALWrappers::cgal_point_to_dealii_point<3>(
738 tet.vertex(1)),
739 CGALWrappers::cgal_point_to_dealii_point<3>(
740 tet.vertex(2)),
741 CGALWrappers::cgal_point_to_dealii_point<3>(
742 tet.vertex(3))}});
743 }
744 }
745 surf1.clear();
746 sm.clear();
747 }
748 surf0.clear();
749 }
750 return vertices;
751 }
752
753 } // namespace internal
754
755
756 template <int dim0, int dim1, int spacedim>
757 std::vector<std::array<Point<spacedim>, dim1 + 1>>
759 const ArrayView<const Point<spacedim>> &vertices0,
760 const ArrayView<const Point<spacedim>> &vertices1,
761 const double tol)
762 {
763 const unsigned int n_vertices0 = vertices0.size();
764 const unsigned int n_vertices1 = vertices1.size();
765
766 Assert(
767 n_vertices0 > 0 || n_vertices1 > 0,
769 "The intersection cannot be computed as at least one of the two cells has no vertices."));
770
771 if constexpr (dim0 == 2 && dim1 == 2 && spacedim == 2)
772 {
773 if (n_vertices0 == 4 && n_vertices1 == 4)
774 {
776 vertices1,
777 tol);
778 }
779 }
780 else if constexpr (dim0 == 2 && dim1 == 1 && spacedim == 2)
781 {
782 if (n_vertices0 == 4 && n_vertices1 == 2)
783 {
785 vertices1,
786 tol);
787 }
788 }
789 else if constexpr (dim0 == 3 && dim1 == 1 && spacedim == 3)
790 {
791 if (n_vertices0 == 8 && n_vertices1 == 2)
792 {
794 vertices1,
795 tol);
796 }
797 }
798 else if constexpr (dim0 == 3 && dim1 == 2 && spacedim == 3)
799 {
800 if (n_vertices0 == 8 && n_vertices1 == 4)
801 {
803 vertices1,
804 tol);
805 }
806 }
807 else if constexpr (dim0 == 3 && dim1 == 3 && spacedim == 3)
808 {
809 if (n_vertices0 == 8 && n_vertices1 == 8)
810 {
812 vertices1,
813 tol);
814 }
815 }
816 else
817 {
818 Assert(false, ExcNotImplemented());
819 return {};
820 }
821 (void)tol;
822 return {};
823 }
824
825
826 template <int dim0, int dim1, int spacedim>
827 std::vector<std::array<Point<spacedim>, dim1 + 1>>
831 const Mapping<dim0, spacedim> & mapping0,
832 const Mapping<dim1, spacedim> & mapping1,
833 const double tol)
834 {
835 Assert(mapping0.get_vertices(cell0).size() ==
836 ReferenceCells::get_hypercube<dim0>().n_vertices(),
838 Assert(mapping1.get_vertices(cell1).size() ==
839 ReferenceCells::get_hypercube<dim1>().n_vertices(),
841
842 const auto &vertices0 =
843 CGALWrappers::get_vertices_in_cgal_order(cell0, mapping0);
844 const auto &vertices1 =
845 CGALWrappers::get_vertices_in_cgal_order(cell1, mapping1);
846
847 return compute_intersection_of_cells<dim0, dim1, spacedim>(vertices0,
848 vertices1,
849 tol);
850 }
851
852# include "intersections.inst"
853
854} // namespace CGALWrappers
855
857
858#else
859
861
862template <int dim0,
863 int dim1,
864 int spacedim,
865 int n_components0,
866 int n_components1>
867std::vector<std::array<Point<spacedim>, dim1 + 1>>
869 const std::array<Point<spacedim>, n_components0> &vertices0,
870 const std::array<Point<spacedim>, n_components1> &vertices1,
871 const double tol)
872{
873 (void)vertices0;
874 (void)vertices1;
875 (void)tol;
876 AssertThrow(false, ExcNeedsCGAL());
877}
878
879template <int dim0, int dim1, int spacedim>
880std::vector<std::array<Point<spacedim>, dim1 + 1>>
884 const Mapping<dim0, spacedim> & mapping0,
885 const Mapping<dim1, spacedim> & mapping1,
886 const double tol)
887{
888 (void)cell0;
889 (void)cell1;
890 (void)mapping0;
891 (void)mapping1;
892 (void)tol;
893 AssertThrow(false, ExcNeedsCGAL());
894}
895
897
898#endif
Abstract base class for mapping classes.
Definition mapping.h:317
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition point.h:112
constexpr void clear()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcNeedsCGAL()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
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)
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_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 > > compute_intersection_triangle_segment(const ArrayView< const Point< 2 > > &triangle, const ArrayView< const Point< 2 > > &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_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 >, 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
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
K::Point_3 CGALPoint3
CGAL::Triangulation_3< K_exact > Triangulation3_exact
CGAL::Constrained_Delaunay_triangulation_2< K, Tds, Itag > CDT
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
K::Point_2 CGALPoint2
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