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