]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix linker error, due to compute_intersection_of_cells() 14362/head
authorJohannes Heinz <johannes.heinz@tuwien.ac.at>
Wed, 19 Oct 2022 07:47:18 +0000 (09:47 +0200)
committerJohannes Heinz <johannes.heinz@tuwien.ac.at>
Mon, 20 Mar 2023 05:10:33 +0000 (06:10 +0100)
Co-authored-by: Marco Feder <marco.feder@sissa.it>
include/deal.II/cgal/intersections.h
include/deal.II/cgal/utilities.h
source/cgal/intersections.cc
source/cgal/intersections.inst.in
tests/cgal/cgal_intersect_simplices_1d_2d.cc
tests/cgal/cgal_intersect_simplices_1d_3d.cc
tests/cgal/cgal_intersect_simplices_2d_2d.cc
tests/cgal/cgal_intersect_simplices_2d_3d.cc
tests/cgal/cgal_intersect_simplices_3d_3d.cc
tests/cgal/cgal_intersect_simplices_with_vertices.cc [new file with mode: 0644]
tests/cgal/cgal_intersect_simplices_with_vertices.output [new file with mode: 0644]

index 0435b04af932dece58e3cdcef045e76b6f8ff02e..1162a25ec9112ae18d553795608614bfe671b328 100644 (file)
@@ -18,6 +18,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/array_view.h>
 #include <deal.II/base/point.h>
 
 #include <deal.II/fe/mapping.h>
@@ -60,21 +61,20 @@ namespace CGALWrappers
 
 
   /**
-   * Same function as above, but working directly with vertices.
+   * Same function as above, but working directly with vertices. It's assumed
+   * that the first vertices are coming from a Triangulation<dim0,spacedim>,
+   * while the other vertices are from a Triangulation<dim1,spacedim>, with
+   * @p dim0 > @p dim1.
+   *
+   * @note The vertices have to be given in CGAL order.
    */
-  template <int dim0, int dim1, int spacedim, int n_vertices0, int n_vertices1>
+  template <int dim0, int dim1, int spacedim>
   std::vector<std::array<Point<spacedim>, dim1 + 1>>
   compute_intersection_of_cells(
-    const std::array<Point<spacedim>, n_vertices0> &vertices0,
-    const std::array<Point<spacedim>, n_vertices1> &vertices1,
-    const double                                    tol = 1e-9)
-  {
-    (void)vertices0;
-    (void)vertices1;
-    (void)tol;
-    Assert(false, ExcMessage("No explicit template instantiation available"));
-    return {};
-  }
+    const ArrayView<const Point<spacedim>> &vertices0,
+    const ArrayView<const Point<spacedim>> &vertices1,
+    const double                            tol = 1e-9);
+
 } // namespace CGALWrappers
 
 DEAL_II_NAMESPACE_CLOSE
index b6f4634012800694226711e4201e8d00f0a9417b..041b10c7b888e25273e2044c47f8b864483e740c 100644 (file)
@@ -601,28 +601,28 @@ namespace CGALWrappers
    * @param mapping Mapping object for the cell.
    * @return  Array of vertices in CGAL order.
    */
-  template <int n_vertices, int dim, int spacedim>
-  std::array<Point<spacedim>, n_vertices>
+  template <int dim, int spacedim>
+  std::vector<Point<spacedim>>
   get_vertices_in_cgal_order(
     const typename dealii::Triangulation<dim, spacedim>::cell_iterator &cell,
     const Mapping<dim, spacedim> &                                      mapping)
   {
     // Elements have to be rectangular or simplices
-    Assert(n_vertices == std::pow(2, dim) || n_vertices == dim + 1,
+    const unsigned int n_vertices = cell->n_vertices();
+    Assert((n_vertices == ReferenceCells::get_hypercube<dim>().n_vertices()) ||
+             (n_vertices == ReferenceCells::get_simplex<dim>().n_vertices()),
            ExcNotImplemented());
-    AssertDimension(mapping.get_vertices(cell).size(), n_vertices);
-
-    std::array<Point<spacedim>, n_vertices> vertices;
 
+    std::vector<Point<spacedim>> ordered_vertices(n_vertices);
     std::copy_n(mapping.get_vertices(cell).begin(),
-                vertices.size(),
-                vertices.begin());
+                n_vertices,
+                ordered_vertices.begin());
 
     if (ReferenceCell::n_vertices_to_type(dim, n_vertices) ==
         ReferenceCells::Quadrilateral)
-      std::swap(vertices[2], vertices[3]);
+      std::swap(ordered_vertices[2], ordered_vertices[3]);
 
-    return vertices;
+    return ordered_vertices;
   }
 
   /**
index a1634706f303e2de69ceed51e40c0f98242a02ad..245103a74815c764bbdf8fba943a0c553cb8f9b9 100644 (file)
@@ -229,126 +229,127 @@ namespace CGALWrappers
                                            CGALSegment2,
                                            CGALTriangle2,
                                            std::vector<CGALPoint2>>>
-    compute_intersection(const std::array<Point<2>, 3> &first_simplex,
-                         const std::array<Point<2>, 3> &second_simplex)
+    compute_intersection_triangle_triangle(
+      const ArrayView<const Point<2>> &triangle0,
+      const ArrayView<const Point<2>> &triangle1)
     {
+      AssertDimension(triangle0.size(), 3);
+      AssertDimension(triangle0.size(), triangle1.size());
+
       std::array<CGALPoint2, 3> pts0, pts1;
-      std::transform(
-        first_simplex.begin(),
-        first_simplex.end(),
-        pts0.begin(),
-        [&](const Point<2> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
-        });
 
-      std::transform(
-        second_simplex.begin(),
-        second_simplex.end(),
-        pts1.begin(),
-        [&](const Point<2> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
-        });
-
-      CGALTriangle2 triangle1{pts0[0], pts0[1], pts0[2]};
-      CGALTriangle2 triangle2{pts1[0], pts1[1], pts1[2]};
-      return convert_boost_to_std(CGAL::intersection(triangle1, triangle2));
-    }
+      std::transform(triangle0.begin(),
+                     triangle0.end(),
+                     pts0.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+      std::transform(triangle1.begin(),
+                     triangle1.end(),
+                     pts1.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
 
+      CGALTriangle2 cgal_triangle0{pts0[0], pts0[1], pts0[2]};
+      CGALTriangle2 cgal_triangle1{pts1[0], pts1[1], pts1[2]};
+      return convert_boost_to_std(
+        CGAL::intersection(cgal_triangle0, cgal_triangle1));
+    }
 
 
     std_cxx17::optional<std_cxx17::variant<CGALPoint2, CGALSegment2>>
-    compute_intersection(const std::array<Point<2>, 3> &first_simplex,
-                         const std::array<Point<2>, 2> &second_simplex)
+    compute_intersection_triangle_segment(
+      const ArrayView<const Point<2>> &triangle,
+      const ArrayView<const Point<2>> &segment)
     {
+      AssertDimension(triangle.size(), 3);
+      AssertDimension(segment.size(), 2);
+
       std::array<CGALPoint2, 3> pts0;
       std::array<CGALPoint2, 2> pts1;
-      std::transform(
-        first_simplex.begin(),
-        first_simplex.end(),
-        pts0.begin(),
-        [&](const Point<2> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
-        });
 
-      std::transform(
-        second_simplex.begin(),
-        second_simplex.end(),
-        pts1.begin(),
-        [&](const Point<2> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
-        });
-
-      CGALTriangle2 triangle{pts0[0], pts0[1], pts0[2]};
-      CGALSegment2  segm{pts1[0], pts1[1]};
-      return convert_boost_to_std(CGAL::intersection(segm, triangle));
+      std::transform(triangle.begin(),
+                     triangle.end(),
+                     pts0.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+      std::transform(segment.begin(),
+                     segment.end(),
+                     pts1.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+      CGALTriangle2 cgal_triangle{pts0[0], pts0[1], pts0[2]};
+      CGALSegment2  cgal_segment{pts1[0], pts1[1]};
+      return convert_boost_to_std(
+        CGAL::intersection(cgal_segment, cgal_triangle));
     }
 
 
 
     // rectangle-rectangle
     std::vector<Polygon_with_holes_2>
-    compute_intersection(const std::array<Point<2>, 4> &first_simplex,
-                         const std::array<Point<2>, 4> &second_simplex)
+    compute_intersection_rect_rect(const ArrayView<const Point<2>> &rectangle0,
+                                   const ArrayView<const Point<2>> &rectangle1)
     {
+      AssertDimension(rectangle0.size(), 4);
+      AssertDimension(rectangle0.size(), rectangle1.size());
+
       std::array<CGALPoint2, 4> pts0, pts1;
-      std::transform(
-        first_simplex.begin(),
-        first_simplex.end(),
-        pts0.begin(),
-        [&](const Point<2> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
-        });
-      std::transform(
-        second_simplex.begin(),
-        second_simplex.end(),
-        pts1.begin(),
-        [&](const Point<2> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
-        });
-      const CGALPolygon first{pts0.begin(), pts0.end()};
-      const CGALPolygon second{pts1.begin(), pts1.end()};
+
+      std::transform(rectangle0.begin(),
+                     rectangle0.end(),
+                     pts0.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+      std::transform(rectangle1.begin(),
+                     rectangle1.end(),
+                     pts1.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+      const CGALPolygon first_poly{pts0.begin(), pts0.end()};
+      const CGALPolygon second_poly{pts1.begin(), pts1.end()};
 
       std::vector<Polygon_with_holes_2> poly_list;
-      CGAL::intersection(first, second, std::back_inserter(poly_list));
+      CGAL::intersection(first_poly,
+                         second_poly,
+                         std::back_inserter(poly_list));
       return poly_list;
     }
 
 
 
     std_cxx17::optional<std_cxx17::variant<CGALPoint3, CGALSegment3>>
-    compute_intersection(const std::array<Point<3>, 2> &first_simplex,
-                         const std::array<Point<3>, 4> &second_simplex)
+    compute_intersection_tetra_segment(
+      const ArrayView<const Point<3>> &tetrahedron,
+      const ArrayView<const Point<3>> &segment)
     {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
+
+      AssertDimension(tetrahedron.size(), 4);
+      AssertDimension(segment.size(), 2);
+
       std::array<CGALPoint3, 4> pts0;
       std::array<CGALPoint3, 2> pts1;
-      std::transform(
-        first_simplex.begin(),
-        first_simplex.end(),
-        pts0.begin(),
-        [&](const Point<3> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
-        });
 
+      std::transform(tetrahedron.begin(),
+                     tetrahedron.end(),
+                     pts0.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
 
-      std::transform(
-        second_simplex.begin(),
-        second_simplex.end(),
-        pts1.begin(),
-        [&](const Point<3> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
-        });
-
-      CGALTetra    tetra{pts0[0], pts0[1], pts0[2], pts0[3]};
-      CGALSegment3 segm{pts1[0], pts1[1]};
-      return convert_boost_to_std(CGAL::intersection(segm, tetra));
+      std::transform(segment.begin(),
+                     segment.end(),
+                     pts1.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
+
+      CGALTetra    cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
+      CGALSegment3 cgal_segment{pts1[0], pts1[1]};
+      return convert_boost_to_std(
+        CGAL::intersection(cgal_segment, cgal_tetrahedron));
 #  else
       Assert(
         false,
         ExcMessage(
           "This function requires a version of CGAL greater or equal than 5.5."));
-      (void)first_simplex;
-      (void)second_simplex;
+      (void)tetrahedron;
+      (void)segment;
       return {};
 #  endif
     }
@@ -359,400 +360,473 @@ namespace CGALWrappers
                                            CGALSegment3,
                                            CGALTriangle3,
                                            std::vector<CGALPoint3>>>
-    compute_intersection(const std::array<Point<3>, 3> &first_simplex,
-                         const std::array<Point<3>, 4> &second_simplex)
+    compute_intersection_tetra_triangle(
+      const ArrayView<const Point<3>> &tetrahedron,
+      const ArrayView<const Point<3>> &triangle)
     {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
+
+      AssertDimension(tetrahedron.size(), 4);
+      AssertDimension(triangle.size(), 3);
+
       std::array<CGALPoint3, 4> pts0;
       std::array<CGALPoint3, 3> pts1;
-      std::transform(
-        first_simplex.begin(),
-        first_simplex.end(),
-        pts0.begin(),
-        [&](const Point<3> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
-        });
 
-      std::transform(
-        second_simplex.begin(),
-        second_simplex.end(),
-        pts1.begin(),
-        [&](const Point<3> &p) {
-          return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
-        });
-      CGALTetra     tetra{pts0[0], pts0[1], pts0[2], pts0[3]};
-      CGALTriangle3 triangle{pts1[0], pts1[1], pts1[2]};
-      return convert_boost_to_std(CGAL::intersection(triangle, tetra));
+      std::transform(tetrahedron.begin(),
+                     tetrahedron.end(),
+                     pts0.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
+
+      std::transform(triangle.begin(),
+                     triangle.end(),
+                     pts1.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
+
+      CGALTetra     cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]};
+      CGALTriangle3 cgal_triangle{pts1[0], pts1[1], pts1[2]};
+      return convert_boost_to_std(
+        CGAL::intersection(cgal_triangle, cgal_tetrahedron));
 #  else
 
       Assert(
         false,
         ExcMessage(
           "This function requires a version of CGAL greater or equal than 5.5."));
-      (void)first_simplex;
-      (void)second_simplex;
+      (void)tetrahedron;
+      (void)triangle;
       return {};
 #  endif
     }
-  } // namespace internal
 
+    // quad-quad
+    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)
+    {
+      AssertDimension(quad0.size(), 4);
+      AssertDimension(quad0.size(), quad1.size());
 
+      const auto intersection_test =
+        internal::compute_intersection_rect_rect(quad0, quad1);
 
-  // Specialization for quads
-  template <>
-  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)
-  {
-    const auto intersection_test =
-      internal::compute_intersection(vertices0, vertices1);
+      if (!intersection_test.empty())
+        {
+          const auto &       poly      = intersection_test[0].outer_boundary();
+          const unsigned int size_poly = poly.size();
+          if (size_poly == 3)
+            {
+              // intersection is a triangle itself, so directly return its
+              // vertices.
+              return {
+                {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)),
+                  CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)),
+                  CGALWrappers::cgal_point_to_dealii_point<2>(
+                    poly.vertex(2))}}};
+            }
+          else if (size_poly >= 4)
+            {
+              // intersection is a polygon, need to triangulate it.
+              std::vector<std::array<Point<2>, 3>> collection;
 
-    if (!intersection_test.empty())
-      {
-        const auto &       poly      = intersection_test[0].outer_boundary();
-        const unsigned int size_poly = poly.size();
-        if (size_poly == 3)
-          {
-            // intersection is a triangle itself, so directly return its
-            // vertices.
-            return {
-              {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)),
-                CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)),
-                CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(2))}}};
-          }
-        else if (size_poly >= 4)
-          {
-            // intersection is a polygon, need to triangulate it.
-            std::vector<std::array<Point<2>, 3>> collection;
-
-            CDT cdt;
-            cdt.insert_constraint(poly.vertices_begin(),
-                                  poly.vertices_end(),
-                                  true);
-
-            internal::mark_domains(cdt);
-            std::array<Point<2>, 3> vertices;
-
-            for (const Face_handle f : cdt.finite_face_handles())
-              {
-                if (f->info().in_domain() &&
-                    CGAL::to_double(cdt.triangle(f).area()) > tol)
-                  {
-                    collection.push_back(
-                      {{CGALWrappers::cgal_point_to_dealii_point<2>(
-                          cdt.triangle(f).vertex(0)),
-                        CGALWrappers::cgal_point_to_dealii_point<2>(
-                          cdt.triangle(f).vertex(1)),
-                        CGALWrappers::cgal_point_to_dealii_point<2>(
-                          cdt.triangle(f).vertex(2))}});
-                  }
-              }
-            return collection;
-          }
-        else
-          {
-            Assert(false, ExcMessage("The polygon is degenerate."));
-            return {};
-          }
-      }
-    else
-      {
-        return {};
-      }
-  }
+              CDT cdt;
+              cdt.insert_constraint(poly.vertices_begin(),
+                                    poly.vertices_end(),
+                                    true);
 
+              internal::mark_domains(cdt);
+              std::array<Point<2>, 3> vertices;
 
+              for (Face_handle f : cdt.finite_face_handles())
+                {
+                  if (f->info().in_domain() &&
+                      CGAL::to_double(cdt.triangle(f).area()) > tol)
+                    {
+                      collection.push_back(
+                        {{CGALWrappers::cgal_point_to_dealii_point<2>(
+                            cdt.triangle(f).vertex(0)),
+                          CGALWrappers::cgal_point_to_dealii_point<2>(
+                            cdt.triangle(f).vertex(1)),
+                          CGALWrappers::cgal_point_to_dealii_point<2>(
+                            cdt.triangle(f).vertex(2))}});
+                    }
+                }
+              return collection;
+            }
+          else
+            {
+              Assert(false, ExcMessage("The polygon is degenerate."));
+              return {};
+            }
+        }
+      else
+        {
+          return {};
+        }
+    }
 
-  // Specialization for quad \cap line
-  template <>
-  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)
-  {
-    std::array<CGALPoint2, 4> pts;
-    std::transform(
-      vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<2> &p) {
-        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
-      });
-
-    CGALPolygon poly(pts.begin(), pts.end());
-
-    CGALSegment2 segm(
-      CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(vertices1[0]),
-      CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(vertices1[1]));
-    CDT cdt;
-    cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true);
-    std::vector<std::array<Point<2>, 2>> vertices;
-    internal::mark_domains(cdt);
-    for (Face_handle f : cdt.finite_face_handles())
-      {
-        if (f->info().in_domain() &&
-            CGAL::to_double(cdt.triangle(f).area()) > tol &&
-            CGAL::do_intersect(segm, cdt.triangle(f)))
-          {
-            const auto intersection = CGAL::intersection(segm, cdt.triangle(f));
-            if (const CGALSegment2 *s =
-                  boost::get<CGALSegment2>(&*intersection))
-              {
-                vertices.push_back(
-                  {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]),
-                    CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}});
-              }
-          }
-      }
-    return vertices;
-  }
+    // Specialization for quad \cap line
+    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)
+    {
+      AssertDimension(quad.size(), 4);
+      AssertDimension(line.size(), 2);
+
+      std::array<CGALPoint2, 4> pts;
+
+      std::transform(quad.begin(),
+                     quad.end(),
+                     pts.begin(),
+                     &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+      CGALPolygon poly(pts.begin(), pts.end());
+
+      CGALSegment2 segm(
+        CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[0]),
+        CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[1]));
+      CDT cdt;
+      cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true);
+      std::vector<std::array<Point<2>, 2>> vertices;
+      internal::mark_domains(cdt);
+      for (Face_handle f : cdt.finite_face_handles())
+        {
+          if (f->info().in_domain() &&
+              CGAL::to_double(cdt.triangle(f).area()) > tol &&
+              CGAL::do_intersect(segm, cdt.triangle(f)))
+            {
+              const auto intersection =
+                CGAL::intersection(segm, cdt.triangle(f));
+              if (const CGALSegment2 *s =
+                    boost::get<CGALSegment2>(&*intersection))
+                {
+                  vertices.push_back(
+                    {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]),
+                      CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}});
+                }
+            }
+        }
 
-  // specialization for hex \cap line
-  template <>
-  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)
-  {
+      return vertices;
+    }
+
+    // specialization for hex \cap line
+    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)
+    {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
-    std::array<CGALPoint3_exact, 8> pts;
-    std::transform(
-      vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<3> &p) {
-        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
-      });
-
-    CGALSegment3_exact segm(
-      CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(vertices1[0]),
-      CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(vertices1[1]));
-
-    // Subdivide the hex into tetrahedrons, and intersect each one of them with
-    // the line
-    std::vector<std::array<Point<3>, 2>> vertices;
-    Triangulation3_exact                 tria;
-    tria.insert(pts.begin(), pts.end());
-    for (const auto &c : tria.finite_cell_handles())
-      {
-        const auto &tet = tria.tetrahedron(c);
-        if (CGAL::do_intersect(segm, tet))
-          {
-            const auto intersection = CGAL::intersection(segm, tet);
-            if (const CGALSegment3_exact *s =
-                  boost::get<CGALSegment3_exact>(&*intersection))
-              {
-                if (s->squared_length() > tol * tol)
-                  {
-                    vertices.push_back(
-                      {{CGALWrappers::cgal_point_to_dealii_point<3>(
-                          s->vertex(0)),
-                        CGALWrappers::cgal_point_to_dealii_point<3>(
-                          s->vertex(1))}});
-                  }
-              }
-          }
-      }
 
-    return vertices;
+      AssertDimension(hexa.size(), 8);
+      AssertDimension(line.size(), 2);
+
+      std::array<CGALPoint3_exact, 8> pts;
+
+      std::transform(
+        hexa.begin(),
+        hexa.end(),
+        pts.begin(),
+        &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
+
+      CGALSegment3_exact cgal_segment(
+        CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[0]),
+        CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[1]));
+
+      // Subdivide the hex into tetrahedrons, and intersect each one of them
+      // with the line
+      std::vector<std::array<Point<3>, 2>> vertices;
+      Triangulation3_exact                 cgal_triangulation;
+      cgal_triangulation.insert(pts.begin(), pts.end());
+      for (const auto &c : cgal_triangulation.finite_cell_handles())
+        {
+          const auto &cgal_tetrahedron = cgal_triangulation.tetrahedron(c);
+          if (CGAL::do_intersect(cgal_segment, cgal_tetrahedron))
+            {
+              const auto intersection =
+                CGAL::intersection(cgal_segment, cgal_tetrahedron);
+              if (const CGALSegment3_exact *s =
+                    boost::get<CGALSegment3_exact>(&*intersection))
+                {
+                  if (s->squared_length() > tol * tol)
+                    {
+                      vertices.push_back(
+                        {{CGALWrappers::cgal_point_to_dealii_point<3>(
+                            s->vertex(0)),
+                          CGALWrappers::cgal_point_to_dealii_point<3>(
+                            s->vertex(1))}});
+                    }
+                }
+            }
+        }
+      return vertices;
 #  else
-    Assert(
-      false,
-      ExcMessage(
-        "This function requires a version of CGAL greater or equal than 5.5."));
-    (void)vertices0;
-    (void)vertices1;
-    (void)tol;
-    return {};
+      Assert(
+        false,
+        ExcMessage(
+          "This function requires a version of CGAL greater or equal than 5.5."));
+      (void)hexa;
+      (void)line;
+      (void)tol;
+      return {};
 #  endif
-  }
+    }
 
-  template <>
-  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)
-  {
+    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)
+    {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
-    std::array<CGALPoint3_exact, 8> pts_hex;
-    std::array<CGALPoint3_exact, 4> pts_quad;
-    std::transform(
-      vertices0.begin(),
-      vertices0.end(),
-      pts_hex.begin(),
-      [&](const Point<3> &p) {
-        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
-      });
-
-    std::transform(
-      vertices1.begin(),
-      vertices1.end(),
-      pts_quad.begin(),
-      [&](const Point<3> &p) {
-        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
-      });
-
-    // Subdivide hex into tetrahedrons
-    std::vector<std::array<Point<3>, 3>> vertices;
-    Triangulation3_exact                 tria;
-    tria.insert(pts_hex.begin(), pts_hex.end());
-
-    // Subdivide quad into triangles
-    Triangulation3_exact tria_quad;
-    tria_quad.insert(pts_quad.begin(), pts_quad.end());
-
-    for (const auto &c : tria.finite_cell_handles())
-      {
-        const auto &tet = tria.tetrahedron(c);
 
-        for (const auto &f : tria_quad.finite_facets())
-          {
-            if (CGAL::do_intersect(tet, tria_quad.triangle(f)))
-              {
-                const auto intersection =
-                  CGAL::intersection(tria_quad.triangle(f), tet);
-
-                if (const CGALTriangle3_exact *t =
-                      boost::get<CGALTriangle3_exact>(&*intersection))
-                  {
-                    if (CGAL::to_double(t->squared_area()) > tol * tol)
-                      {
-                        vertices.push_back(
-                          {{cgal_point_to_dealii_point<3>((*t)[0]),
-                            cgal_point_to_dealii_point<3>((*t)[1]),
-                            cgal_point_to_dealii_point<3>((*t)[2])}});
-                      }
-                  }
-
-                if (const std::vector<CGALPoint3_exact> *vps =
-                      boost::get<std::vector<CGALPoint3_exact>>(&*intersection))
-                  {
-                    Triangulation3_exact tria_inter;
-                    tria_inter.insert(vps->begin(), vps->end());
-
-                    for (auto it = tria_inter.finite_facets_begin();
-                         it != tria_inter.finite_facets_end();
-                         ++it)
-                      {
-                        const auto triangle = tria_inter.triangle(*it);
-                        if (CGAL::to_double(triangle.squared_area()) >
-                            tol * tol)
-                          {
-                            std::array<Point<3>, 3> verts = {
-                              {CGALWrappers::cgal_point_to_dealii_point<3>(
-                                 triangle[0]),
-                               CGALWrappers::cgal_point_to_dealii_point<3>(
-                                 triangle[1]),
-                               CGALWrappers::cgal_point_to_dealii_point<3>(
-                                 triangle[2])}};
-
-                            vertices.push_back(verts);
-                          }
-                      }
-                  }
-              }
-          }
-      }
+      AssertDimension(hexa.size(), 8);
+      AssertDimension(quad.size(), 4);
+
+      std::array<CGALPoint3_exact, 8> pts_hex;
+      std::array<CGALPoint3_exact, 4> pts_quad;
+
+      std::transform(
+        hexa.begin(),
+        hexa.end(),
+        pts_hex.begin(),
+        &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
+
+      std::transform(
+        quad.begin(),
+        quad.end(),
+        pts_quad.begin(),
+        &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
+
+      // Subdivide hex into tetrahedrons
+      std::vector<std::array<Point<3>, 3>> vertices;
+      Triangulation3_exact                 triangulation_hexa;
+      triangulation_hexa.insert(pts_hex.begin(), pts_hex.end());
+
+      // Subdivide quad into triangles
+      Triangulation3_exact triangulation_quad;
+      triangulation_quad.insert(pts_quad.begin(), pts_quad.end());
+
+      for (const auto &c : triangulation_hexa.finite_cell_handles())
+        {
+          const auto &tet = triangulation_hexa.tetrahedron(c);
+
+          for (const auto &f : triangulation_quad.finite_facets())
+            {
+              if (CGAL::do_intersect(tet, triangulation_quad.triangle(f)))
+                {
+                  const auto intersection =
+                    CGAL::intersection(triangulation_quad.triangle(f), tet);
+
+                  if (const CGALTriangle3_exact *t =
+                        boost::get<CGALTriangle3_exact>(&*intersection))
+                    {
+                      if (CGAL::to_double(t->squared_area()) > tol * tol)
+                        {
+                          vertices.push_back(
+                            {{cgal_point_to_dealii_point<3>((*t)[0]),
+                              cgal_point_to_dealii_point<3>((*t)[1]),
+                              cgal_point_to_dealii_point<3>((*t)[2])}});
+                        }
+                    }
+
+                  if (const std::vector<CGALPoint3_exact> *vps =
+                        boost::get<std::vector<CGALPoint3_exact>>(
+                          &*intersection))
+                    {
+                      Triangulation3_exact tria_inter;
+                      tria_inter.insert(vps->begin(), vps->end());
+
+                      for (auto it = tria_inter.finite_facets_begin();
+                           it != tria_inter.finite_facets_end();
+                           ++it)
+                        {
+                          const auto triangle = tria_inter.triangle(*it);
+                          if (CGAL::to_double(triangle.squared_area()) >
+                              tol * tol)
+                            {
+                              std::array<Point<3>, 3> verts = {
+                                {CGALWrappers::cgal_point_to_dealii_point<3>(
+                                   triangle[0]),
+                                 CGALWrappers::cgal_point_to_dealii_point<3>(
+                                   triangle[1]),
+                                 CGALWrappers::cgal_point_to_dealii_point<3>(
+                                   triangle[2])}};
+
+                              vertices.push_back(verts);
+                            }
+                        }
+                    }
+                }
+            }
+        }
 
-    return vertices;
+      return vertices;
 #  else
-    Assert(
-      false,
-      ExcMessage(
-        "This function requires a version of CGAL greater or equal than 5.5."));
-    (void)vertices0;
-    (void)vertices1;
-    (void)tol;
-    return {};
+      Assert(
+        false,
+        ExcMessage(
+          "This function requires a version of CGAL greater or equal than 5.5."));
+      (void)hexa;
+      (void)quad;
+      (void)tol;
+      return {};
 #  endif
-  }
+    }
+
+    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)
+    {
+      AssertDimension(hexa0.size(), 8);
+      AssertDimension(hexa0.size(), hexa1.size());
 
+      std::array<CGALPoint3_inexact, 8> pts_hex0;
+      std::array<CGALPoint3_inexact, 8> pts_hex1;
 
+      std::transform(
+        hexa0.begin(),
+        hexa0.end(),
+        pts_hex0.begin(),
+        &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact, 3>);
+
+      std::transform(
+        hexa1.begin(),
+        hexa1.end(),
+        pts_hex1.begin(),
+        &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact, 3>);
+
+      Surface_mesh surf0, surf1, sm;
+      // Subdivide hex into tetrahedrons
+      std::vector<std::array<Point<3>, 4>> vertices;
+      Triangulation3_inexact               tria0, tria1;
+
+      tria0.insert(pts_hex0.begin(), pts_hex0.end());
+      tria1.insert(pts_hex1.begin(), pts_hex1.end());
+
+      for (const auto &c0 : tria0.finite_cell_handles())
+        {
+          const auto &tet0  = tria1.tetrahedron(c0);
+          const auto &tetg0 = CGAL::make_tetrahedron(tet0.vertex(0),
+                                                     tet0.vertex(1),
+                                                     tet0.vertex(2),
+                                                     tet0.vertex(3),
+                                                     surf0);
+          (void)tetg0; // instead of C++ 17s [[maybe unused]]
+          for (const auto &c1 : tria1.finite_cell_handles())
+            {
+              const auto &tet1  = tria1.tetrahedron(c1);
+              const auto &tetg1 = CGAL::make_tetrahedron(tet1.vertex(0),
+                                                         tet1.vertex(1),
+                                                         tet1.vertex(2),
+                                                         tet1.vertex(3),
+                                                         surf1);
+              (void)tetg1; // instead of C++ 17s [[maybe unused]]
+              namespace PMP = CGAL::Polygon_mesh_processing;
+              const bool test_intersection =
+                PMP::corefine_and_compute_intersection(surf0, surf1, sm);
+              if (PMP::volume(sm) > tol && test_intersection)
+                {
+                  // Collect tetrahedrons
+                  Triangulation3_inexact triangulation_hexa;
+                  triangulation_hexa.insert(sm.points().begin(),
+                                            sm.points().end());
+                  for (const auto &c : triangulation_hexa.finite_cell_handles())
+                    {
+                      const auto &tet = triangulation_hexa.tetrahedron(c);
+                      vertices.push_back(
+                        {{CGALWrappers::cgal_point_to_dealii_point<3>(
+                            tet.vertex(0)),
+                          CGALWrappers::cgal_point_to_dealii_point<3>(
+                            tet.vertex(1)),
+                          CGALWrappers::cgal_point_to_dealii_point<3>(
+                            tet.vertex(2)),
+                          CGALWrappers::cgal_point_to_dealii_point<3>(
+                            tet.vertex(3))}});
+                    }
+                }
+              surf1.clear();
+              sm.clear();
+            }
+          surf0.clear();
+        }
+      return vertices;
+    }
+
+  } // namespace internal
 
-  template <>
-  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)
+  template <int dim0, int dim1, int spacedim>
+  std::vector<std::array<Point<spacedim>, dim1 + 1>>
+  compute_intersection_of_cells(
+    const ArrayView<const Point<spacedim>> &vertices0,
+    const ArrayView<const Point<spacedim>> &vertices1,
+    const double                            tol)
   {
-    std::array<CGALPoint3_inexact, 8> pts_hex0;
-    std::array<CGALPoint3_inexact, 8> pts_hex1;
-    std::transform(
-      vertices0.begin(),
-      vertices0.end(),
-      pts_hex0.begin(),
-      [&](const Point<3> &p) {
-        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
-      });
-
-    std::transform(
-      vertices1.begin(),
-      vertices1.end(),
-      pts_hex1.begin(),
-      [&](const Point<3> &p) {
-        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
-      });
-
-
-    Surface_mesh surf0, surf1, sm;
-    // Subdivide hex into tetrahedrons
-    std::vector<std::array<Point<3>, 4>> vertices;
-    Triangulation3_inexact               tria0, tria1;
-
-    tria0.insert(pts_hex0.begin(), pts_hex0.end());
-    tria1.insert(pts_hex1.begin(), pts_hex1.end());
-
-    for (const auto &c0 : tria0.finite_cell_handles())
+    const unsigned int n_vertices0 = vertices0.size();
+    const unsigned int n_vertices1 = vertices1.size();
+
+    Assert(
+      n_vertices0 > 0 || n_vertices1 > 0,
+      ExcMessage(
+        "The intersection cannot be computed as at least one of the two cells has no vertices."));
+
+    if constexpr (dim0 == 2 && dim1 == 2 && spacedim == 2)
       {
-        const auto &                 tet0 = tria1.tetrahedron(c0);
-        [[maybe_unused]] const auto &tetg0 =
-          CGAL::make_tetrahedron(tet0.vertex(0),
-                                 tet0.vertex(1),
-                                 tet0.vertex(2),
-                                 tet0.vertex(3),
-                                 surf0);
-        for (const auto &c1 : tria1.finite_cell_handles())
+        if (n_vertices0 == 4 && n_vertices1 == 4)
           {
-            const auto &                 tet1 = tria1.tetrahedron(c1);
-            [[maybe_unused]] const auto &tetg1 =
-              CGAL::make_tetrahedron(tet1.vertex(0),
-                                     tet1.vertex(1),
-                                     tet1.vertex(2),
-                                     tet1.vertex(3),
-                                     surf1);
-            namespace PMP = CGAL::Polygon_mesh_processing;
-            const bool test_intersection =
-              PMP::corefine_and_compute_intersection(surf0, surf1, sm);
-            if (PMP::volume(sm) > tol && test_intersection)
-              {
-                // Collect tetrahedrons
-                Triangulation3_inexact tria;
-                tria.insert(sm.points().begin(), sm.points().end());
-                for (const auto &c : tria.finite_cell_handles())
-                  {
-                    const auto &tet = tria.tetrahedron(c);
-                    vertices.push_back(
-                      {{CGALWrappers::cgal_point_to_dealii_point<3>(
-                          tet.vertex(0)),
-                        CGALWrappers::cgal_point_to_dealii_point<3>(
-                          tet.vertex(1)),
-                        CGALWrappers::cgal_point_to_dealii_point<3>(
-                          tet.vertex(2)),
-                        CGALWrappers::cgal_point_to_dealii_point<3>(
-                          tet.vertex(3))}});
-                  }
-              }
-            surf1.clear();
-            sm.clear();
+            return internal::compute_intersection_quad_quad(vertices0,
+                                                            vertices1,
+                                                            tol);
           }
-        surf0.clear();
       }
-    return vertices;
+    else if constexpr (dim0 == 2 && dim1 == 1 && spacedim == 2)
+      {
+        if (n_vertices0 == 4 && n_vertices1 == 2)
+          {
+            return internal::compute_intersection_quad_line(vertices0,
+                                                            vertices1,
+                                                            tol);
+          }
+      }
+    else if constexpr (dim0 == 3 && dim1 == 1 && spacedim == 3)
+      {
+        if (n_vertices0 == 8 && n_vertices1 == 2)
+          {
+            return internal::compute_intersection_hexa_line(vertices0,
+                                                            vertices1,
+                                                            tol);
+          }
+      }
+    else if constexpr (dim0 == 3 && dim1 == 2 && spacedim == 3)
+      {
+        if (n_vertices0 == 8 && n_vertices1 == 4)
+          {
+            return internal::compute_intersection_hexa_quad(vertices0,
+                                                            vertices1,
+                                                            tol);
+          }
+      }
+    else if constexpr (dim0 == 3 && dim1 == 3 && spacedim == 3)
+      {
+        if (n_vertices0 == 8 && n_vertices1 == 8)
+          {
+            return internal::compute_intersection_hexa_hexa(vertices0,
+                                                            vertices1,
+                                                            tol);
+          }
+      }
+    else
+      {
+        Assert(false, ExcNotImplemented());
+        return {};
+      }
+    (void)tol;
+    return {};
   }
 
 
-
   template <int dim0, int dim1, int spacedim>
   std::vector<std::array<Point<spacedim>, dim1 + 1>>
   compute_intersection_of_cells(
@@ -762,25 +836,21 @@ namespace CGALWrappers
     const Mapping<dim1, spacedim> &                              mapping1,
     const double                                                 tol)
   {
-    Assert(mapping0.get_vertices(cell0).size() == std::pow(2, dim0),
+    Assert(mapping0.get_vertices(cell0).size() ==
+             ReferenceCells::get_hypercube<dim0>().n_vertices(),
            ExcNotImplemented());
-    Assert(mapping1.get_vertices(cell1).size() == std::pow(2, dim1),
+    Assert(mapping1.get_vertices(cell1).size() ==
+             ReferenceCells::get_hypercube<dim1>().n_vertices(),
            ExcNotImplemented());
 
-    const auto vertices0 =
-      CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim0)>(
-        cell0, mapping0);
-    const auto vertices1 =
-      CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim1)>(
-        cell1, mapping1);
-
-    return compute_intersection_of_cells<dim0,
-                                         dim1,
-                                         spacedim,
-                                         Utilities::pow(2, dim0),
-                                         Utilities::pow(2, dim1)>(vertices0,
-                                                                  vertices1,
-                                                                  tol);
+    const auto &vertices0 =
+      CGALWrappers::get_vertices_in_cgal_order(cell0, mapping0);
+    const auto &vertices1 =
+      CGALWrappers::get_vertices_in_cgal_order(cell1, mapping1);
+
+    return compute_intersection_of_cells<dim0, dim1, spacedim>(vertices0,
+                                                               vertices1,
+                                                               tol);
   }
 
 #  include "intersections.inst"
index d550252e0f72b0402025d22f4c543b2f6b8e1d94..3024147791ca246a8dc9597dbd399facf2d65eec 100644 (file)
@@ -27,5 +27,11 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS)
       const Mapping<dim1, spacedim> &                              mapping1,
       const double                                                 tol);
 
+    template std::vector<std::array<Point<spacedim>, dim1 + 1>>
+    compute_intersection_of_cells<dim0, dim1, spacedim>(
+      const ArrayView<const Point<spacedim>> &vertices0,
+      const ArrayView<const Point<spacedim>> &vertices1,
+      const double                            tol);
+
 #endif
   }
index b0d00ba33cb8cf32331b19bda60c03059ad9ac08..fdb33fa109a2f4d4c2152d5513ae86e8ea4d6830 100644 (file)
@@ -42,11 +42,8 @@ test_inside_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1)
   const auto   cell1         = tria1.begin_active();
   // cell1->vertex(1)+=Point<2>(1.0,1.);
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<2, 1, 2>(cell0,
-                                                         cell1,
-                                                         MappingQ1<2>(),
-                                                         MappingQ1<1, 2>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<2>(), MappingQ1<1, 2>());
 
   const auto   quad = qgauss.mapped_quadrature(vec_of_arrays);
   const double sum =
@@ -67,11 +64,8 @@ test_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1)
   const auto   cell0         = tria0.begin_active();
   const auto   cell1         = tria1.begin_active();
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<2, 1, 2>(cell0,
-                                                         cell1,
-                                                         MappingQ1<2>(),
-                                                         MappingQ1<1, 2>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<2>(), MappingQ1<1, 2>());
 
   const auto   quad = qgauss.mapped_quadrature(vec_of_arrays);
   const double sum =
@@ -91,11 +85,8 @@ test_failing_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1)
   const auto   cell0         = tria0.begin_active();
   const auto   cell1         = tria1.begin_active();
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<2, 1, 2>(cell0,
-                                                         cell1,
-                                                         MappingQ1<2>(),
-                                                         MappingQ1<1, 2>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<2>(), MappingQ1<1, 2>());
 
   const auto   quad = qgauss.mapped_quadrature(vec_of_arrays);
   const double sum =
index d7bd153ac43bb9cf94ba1dd87d3ec7e0a60d6332..7e64cc84266fc021a06cd7e1f555128d9fbff945 100644 (file)
@@ -41,10 +41,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<1, 3> &tria1)
     (tria1.begin_active()->vertex(1) - tria1.begin_active()->vertex(0)).norm();
 
   const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<3, 1, 3>(tria0.begin_active(),
-                                                         tria1.begin_active(),
-                                                         MappingQ1<3>(),
-                                                         MappingQ1<1, 3>());
+    CGALWrappers::compute_intersection_of_cells(tria0.begin_active(),
+                                                tria1.begin_active(),
+                                                MappingQ1<3>(),
+                                                MappingQ1<1, 3>());
 
   const auto   quad = QGaussSimplex<1>(1).mapped_quadrature(vec_of_arrays);
   const double sum =
@@ -65,11 +65,8 @@ test_intersection(Triangulation<3> &tria0, Triangulation<1, 3> &tria1)
   cell1->vertex(1)              = Point<3>(1.5, 1.5, 1.5);
   const double expected_measure = std::sqrt(3.);
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<3, 1, 3>(tria0.begin_active(),
-                                                         cell1,
-                                                         MappingQ1<3>(),
-                                                         MappingQ1<1, 3>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    tria0.begin_active(), cell1, MappingQ1<3>(), MappingQ1<1, 3>());
 
   const auto   quad = QGaussSimplex<1>(1).mapped_quadrature(vec_of_arrays);
   const double sum =
@@ -90,11 +87,8 @@ test_failing_intersection(Triangulation<3> &tria0, Triangulation<1, 3> &tria1)
   const auto   cell1            = tria1.begin_active();
   const double expected_measure = 0.;
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<3, 1, 3>(cell0,
-                                                         cell1,
-                                                         MappingQ1<3>(),
-                                                         MappingQ1<1, 3>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<3>(), MappingQ1<1, 3>());
 
   const auto   quad = QGaussSimplex<1>(1).mapped_quadrature(vec_of_arrays);
   const double sum =
index c4e99b60b6845849aa65be58ec66b88cf7038065..1d77622d5163cbf70bb4f7dea2cc36cec4150f39 100644 (file)
@@ -41,11 +41,8 @@ test_inside_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1)
   const auto   cell0         = tria0.begin_active();
   const auto   cell1         = tria1.begin_active();
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<2, 2, 2>(cell0,
-                                                         cell1,
-                                                         MappingQ1<2>(),
-                                                         MappingQ1<2>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<2>(), MappingQ1<2>());
 
 
   const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
@@ -67,11 +64,8 @@ test_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1)
   const auto cell0 = tria0.begin_active();
   const auto cell1 = tria1.begin_active();
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<2, 2, 2>(cell0,
-                                                         cell1,
-                                                         MappingQ1<2>(),
-                                                         MappingQ1<2>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<2>(), MappingQ1<2>());
 
   const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
   const double sum =
@@ -92,11 +86,8 @@ test_failing_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1)
   const auto cell0 = tria0.begin_active();
   const auto cell1 = tria1.begin_active();
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<2, 2, 2>(cell0,
-                                                         cell1,
-                                                         MappingQ1<2>(),
-                                                         MappingQ1<2>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<2>(), MappingQ1<2>());
 
   const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
   const double sum =
index be43f35fd7ff47ce7e8c2f32abe6e691beae168a..225637f71c8ea953d37bf042cd2db1899893bd72 100644 (file)
@@ -58,10 +58,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
   {
     const double expected_measure = 0.36;
     const auto   vec_of_arrays =
-      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                           tria1.begin_active(),
-                                                           MappingQ1<3>(),
-                                                           MappingQ1<2, 3>());
+      CGALWrappers::compute_intersection_of_cells(tria0.begin_active(),
+                                                  tria1.begin_active(),
+                                                  MappingQ1<3>(),
+                                                  MappingQ1<2, 3>());
 
     const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
     const double sum =
@@ -73,10 +73,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
     GridTools::transform(swap_coordinates(1, 2), tria1);
     const double expected_measure = 0.36;
     const auto   vec_of_arrays =
-      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                           tria1.begin_active(),
-                                                           MappingQ1<3>(),
-                                                           MappingQ1<2, 3>());
+      CGALWrappers::compute_intersection_of_cells(tria0.begin_active(),
+                                                  tria1.begin_active(),
+                                                  MappingQ1<3>(),
+                                                  MappingQ1<2, 3>());
 
     const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
     const double sum =
@@ -88,10 +88,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
     GridTools::transform(swap_coordinates(0, 1), tria1);
     const double expected_measure = 0.36;
     const auto   vec_of_arrays =
-      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                           tria1.begin_active(),
-                                                           MappingQ1<3>(),
-                                                           MappingQ1<2, 3>());
+      CGALWrappers::compute_intersection_of_cells(tria0.begin_active(),
+                                                  tria1.begin_active(),
+                                                  MappingQ1<3>(),
+                                                  MappingQ1<2, 3>());
 
     const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
     const double sum =
@@ -112,10 +112,10 @@ test_intersection(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
   {
     const double expected_measure = 0.25;
     const auto   vec_of_arrays =
-      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                           tria1.begin_active(),
-                                                           MappingQ1<3>(),
-                                                           MappingQ1<2, 3>());
+      CGALWrappers::compute_intersection_of_cells(tria0.begin_active(),
+                                                  tria1.begin_active(),
+                                                  MappingQ1<3>(),
+                                                  MappingQ1<2, 3>());
 
     const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
     const double sum =
@@ -127,10 +127,10 @@ test_intersection(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
     GridTools::transform(swap_coordinates(1, 2), tria1);
     const double expected_measure = 0.25;
     const auto   vec_of_arrays =
-      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                           tria1.begin_active(),
-                                                           MappingQ1<3>(),
-                                                           MappingQ1<2, 3>());
+      CGALWrappers::compute_intersection_of_cells(tria0.begin_active(),
+                                                  tria1.begin_active(),
+                                                  MappingQ1<3>(),
+                                                  MappingQ1<2, 3>());
 
     const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
     const double sum =
@@ -142,10 +142,10 @@ test_intersection(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
     GridTools::transform(swap_coordinates(0, 1), tria1);
     const double expected_measure = 0.25;
     const auto   vec_of_arrays =
-      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                           tria1.begin_active(),
-                                                           MappingQ1<3>(),
-                                                           MappingQ1<2, 3>());
+      CGALWrappers::compute_intersection_of_cells(tria0.begin_active(),
+                                                  tria1.begin_active(),
+                                                  MappingQ1<3>(),
+                                                  MappingQ1<2, 3>());
 
     const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
     const double sum =
index d19f3f3ef9c00f701b5e9504170ddfb397b229aa..e181eba7f9313c406f6014afa6f71780c0a27096 100644 (file)
@@ -36,11 +36,8 @@ test_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1)
   const auto              cell0            = tria0.begin_active();
   const auto              cell1            = tria1.begin_active();
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<3, 3, 3>(cell0,
-                                                         cell1,
-                                                         MappingQ1<3>(),
-                                                         MappingQ1<3>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<3>(), MappingQ1<3>());
 
   const auto   quad = QGaussSimplex<3>(1).mapped_quadrature(vec_of_arrays);
   const double sum =
@@ -60,11 +57,8 @@ test_failing_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1)
   const auto              cell0            = tria0.begin_active();
   const auto              cell1            = tria1.begin_active();
 
-  const auto vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<3, 3, 3>(cell0,
-                                                         cell1,
-                                                         MappingQ1<3>(),
-                                                         MappingQ1<3>());
+  const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells(
+    cell0, cell1, MappingQ1<3>(), MappingQ1<3>());
 
   const auto   quad = QGaussSimplex<3>(1).mapped_quadrature(vec_of_arrays);
   const double sum =
diff --git a/tests/cgal/cgal_intersect_simplices_with_vertices.cc b/tests/cgal/cgal_intersect_simplices_with_vertices.cc
new file mode 100644 (file)
index 0000000..a494e91
--- /dev/null
@@ -0,0 +1,136 @@
+// ---------------------------------------------------------------------
+
+// Copyright (C) 2022 by the deal.II authors
+
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// This test is the same as cgal_intersection_simplices_2d_2d() but directly
+// working with vertices While the functionality of the intersections is tested
+// with the iterator interface this test ensures there are no linker errors if
+// the vertex version is called directly.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/cgal/intersections.h>
+#include <deal.II/cgal/utilities.h>
+
+#include "../tests.h"
+
+using namespace CGALWrappers;
+
+void
+test_inside_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1)
+{
+  GridGenerator::hyper_cube(tria0, -1., 1.);
+  GridGenerator::hyper_cube(tria1, -0.45, 0.45);
+  GridTools::rotate(numbers::PI_4, tria1);
+  const double expected_area = 0.81;
+  const auto   cell0         = tria0.begin_active();
+  const auto   cell1         = tria1.begin_active();
+
+  const auto &vertices0 =
+    CGALWrappers::get_vertices_in_cgal_order(cell0, MappingQ1<2>());
+  const auto &vertices1 =
+    CGALWrappers::get_vertices_in_cgal_order(cell1, MappingQ1<2>());
+
+  const auto &vec_of_arrays =
+    CGALWrappers::compute_intersection_of_cells<2, 2, 2>(vertices0, vertices1);
+
+  const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
+  const double sum =
+    std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.);
+  assert(std::abs(expected_area - sum) < 1e-15);
+  deallog << "OK" << std::endl;
+}
+
+
+
+void
+test_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1)
+{
+  GridGenerator::hyper_cube(tria0, -1., 1.);
+  GridGenerator::hyper_cube(tria1, .5, 1.45);
+  const double expected_area = 0.25;
+
+  const auto cell0 = tria0.begin_active();
+  const auto cell1 = tria1.begin_active();
+
+  const auto &vertices0 =
+    CGALWrappers::get_vertices_in_cgal_order(cell0, MappingQ1<2>());
+  const auto &vertices1 =
+    CGALWrappers::get_vertices_in_cgal_order(cell1, MappingQ1<2>());
+
+  const auto &vec_of_arrays =
+    CGALWrappers::compute_intersection_of_cells<2, 2, 2>(vertices0, vertices1);
+
+
+  const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
+  const double sum =
+    std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.);
+  assert(std::abs(expected_area - sum) < 1e-15);
+  deallog << "OK" << std::endl;
+}
+
+
+
+void
+test_failing_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1)
+{
+  GridGenerator::hyper_cube(tria0, -1., 1.);
+  GridGenerator::hyper_cube(tria1, 1.5, 2.5);
+  const double expected_area = 0.;
+
+  const auto cell0 = tria0.begin_active();
+  const auto cell1 = tria1.begin_active();
+
+  const auto &vertices0 =
+    CGALWrappers::get_vertices_in_cgal_order(cell0, MappingQ1<2>());
+  const auto &vertices1 =
+    CGALWrappers::get_vertices_in_cgal_order(cell1, MappingQ1<2>());
+
+  const auto &vec_of_arrays =
+    CGALWrappers::compute_intersection_of_cells<2, 2, 2>(vertices0, vertices1);
+
+  const auto   quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays);
+  const double sum =
+    std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.);
+  assert(std::abs(expected_area - sum) < 1e-15);
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+  Triangulation<2> tria0;
+  Triangulation<2> tria1;
+
+  test_inside_intersection(tria0, tria1);
+
+  tria0.clear();
+  tria1.clear();
+
+  test_intersection(tria0, tria1);
+
+  tria0.clear();
+  tria1.clear();
+
+  test_failing_intersection(tria0, tria1);
+}
diff --git a/tests/cgal/cgal_intersect_simplices_with_vertices.output b/tests/cgal/cgal_intersect_simplices_with_vertices.output
new file mode 100644 (file)
index 0000000..fb71de2
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.