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