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