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