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