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