add additional interface and fix subminor version

diff --git a/cmake/modules/FindCGAL.cmake b/cmake/modules/FindCGAL.cmake
index 820f27a57f2d4da8342f75afe2752ab94fd076ec..34692bee59b627d099078679c39723e85ab94cf8 100644 (file)
--- a/cmake/modules/FindCGAL.cmake
+++ b/cmake/modules/FindCGAL.cmake
@@ -71,7 +71,7 @@ IF(DEAL_II_HAVE_CXX17 AND NOT FEATURE_BOOST_BUNDLED_CONFIGURED)
         CGAL_VERSION_SUBMINOR "${CGAL_VERSION}"
       )
 
-      IF("${CGAL_VERSION_SUBMINOR}" STREQUAL "")
+      IF("${CGAL_VERSION_SUBMINOR}" MATCHES "^(|${CGAL_VERSION})$")
         SET(CGAL_VERSION_SUBMINOR "0")
       ENDIF()
     ENDIF()

diff --git a/include/deal.II/cgal/intersections.h b/include/deal.II/cgal/intersections.h
index 5a7f3b4b17e3fe282297c6762ab3f7f0fb0b278b..8254d10e35c74629830b033fed2c60fc3dd21c53 100644 (file)
--- a/include/deal.II/cgal/intersections.h
+++ b/include/deal.II/cgal/intersections.h
@@ -28,6 +28,24 @@
 DEAL_II_NAMESPACE_OPEN
 namespace CGALWrappers
 {
+  template <int dim0,
+            int dim1,
+            int spacedim,
+            int n_components0,
+            int n_components1>
+  std::vector<std::array<Point<spacedim>, dim1 + 1>>
+  compute_intersection_of_cells(
+    const std::array<Point<spacedim>, n_components0> &vertices0,
+    const std::array<Point<spacedim>, n_components1> &vertices1,
+    const double                                      tol)
+  {
+    (void)vertices0;
+    (void)vertices1;
+    (void)tol;
+    Assert(false, ExcMessage("No explicit template instantiation available"));
+    return {};
+  }
+
   /**
    * Given two deal.II affine cells, compute the intersection and return a
    * vector of simplices, each one identified by an array of deal.II Points. All

diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h
index 916a8e44378b93bad90ee3afcea95460d5bd2086..4e8bc23a80ab1b382a90e032eba6e76bed39f436 100644 (file)
--- a/include/deal.II/cgal/utilities.h
+++ b/include/deal.II/cgal/utilities.h
@@ -592,6 +592,39 @@ namespace CGALWrappers
         Assert(false, ExcInternalError());
       }
   }
+
+
+  /**
+   * Get vertices of cell in CGAL ordering
+   *
+   * @param [in] cell A cell_iterator to a deal.II cell.
+   * @param [in] mapping Mapping object for the cell.
+   * @return [out] Array of vertices in CGAL order.
+   */
+  template <int n_vertices, int dim, int spacedim>
+  std::array<Point<spacedim>, n_vertices>
+  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,
+           ExcNotImplemented());
+    AssertDimension(mapping.get_vertices(cell).size(), n_vertices);
+
+    std::array<Point<spacedim>, n_vertices> vertices;
+
+    std::copy_n(mapping.get_vertices(cell).begin(),
+                vertices.size(),
+                vertices.begin());
+
+    if (ReferenceCell::n_vertices_to_type(dim, n_vertices) ==
+        ReferenceCells::Quadrilateral)
+      std::swap(vertices[2], vertices[3]);
+
+    return vertices;
+  }
+
 } // namespace CGALWrappers
 #  endif
 

diff --git a/source/cgal/CMakeLists.txt b/source/cgal/CMakeLists.txt
index 2c951a6b36344cd4d7c84b2d1003e0aa04449083..3d164833397b42ad5f4af15dd240957e4b80cf73 100644 (file)
--- a/source/cgal/CMakeLists.txt
+++ b/source/cgal/CMakeLists.txt
@@ -26,6 +26,7 @@ SET(_src
 
 SET(_inst
 surface_mesh.inst.in
+intersections.inst.in
 )
 
 FILE(GLOB _header

diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc
index a9b23fe29ba94334b8682e93b6d00af5c77ceb0f..7a859849b92fb583545cb907eb6370a1b4d7f53f 100644 (file)
--- a/source/cgal/intersections.cc
+++ b/source/cgal/intersections.cc
@@ -19,7 +19,6 @@
 
 #include <algorithm>
 
-
 #ifdef DEAL_II_WITH_CGAL
 
 #  include <deal.II/base/quadrature_lib.h>
@@ -56,8 +55,6 @@
 #  include <fstream>
 #  include <type_traits>
 
-
-
 DEAL_II_NAMESPACE_OPEN
 
 namespace CGALWrappers
@@ -144,8 +141,6 @@ namespace CGALWrappers
         }
     }
 
-
-
     void
     mark_domains(CDT &cdt)
     {
@@ -205,8 +200,6 @@ namespace CGALWrappers
       return CGAL::intersection(triangle1, triangle2);
     }
 
-
-
     boost::optional<boost::variant<CGALPoint2, CGALSegment2>>
     compute_intersection(const std::array<Point<2>, 3> &first_simplex,
                          const std::array<Point<2>, 2> &second_simplex)
@@ -234,8 +227,6 @@ namespace CGALWrappers
       return CGAL::intersection(segm, triangle);
     }
 
-
-
     // rectangle-rectangle
     std::vector<Polygon_with_holes_2>
     compute_intersection(const std::array<Point<2>, 4> &first_simplex,
@@ -264,14 +255,11 @@ namespace CGALWrappers
       return poly_list;
     }
 
-
-
     boost::optional<boost::variant<CGALPoint3, CGALSegment3>>
     compute_intersection(const std::array<Point<3>, 2> &first_simplex,
                          const std::array<Point<3>, 4> &second_simplex)
     {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 1, 5)
-
       std::array<CGALPoint3, 4> pts0;
       std::array<CGALPoint3, 2> pts1;
       std::transform(
@@ -294,9 +282,7 @@ namespace CGALWrappers
       CGALTetra    tetra{pts0[0], pts0[1], pts0[2], pts0[3]};
       CGALSegment3 segm{pts1[0], pts1[1]};
       return CGAL::intersection(segm, tetra);
-
 #  else
-
       Assert(
         false,
         ExcMessage(
@@ -307,8 +293,6 @@ namespace CGALWrappers
 #  endif
     }
 
-
-
     // tetra, triangle
     boost::optional<boost::variant<CGALPoint3,
                                    CGALSegment3,
@@ -351,33 +335,14 @@ namespace CGALWrappers
     }
   } // namespace internal
 
-
-
   // Specialization for quads
   template <>
   std::vector<std::array<Point<2>, 3>>
-  compute_intersection_of_cells<2, 2, 2>(
-    const typename Triangulation<2, 2>::cell_iterator &cell0,
-    const typename Triangulation<2, 2>::cell_iterator &cell1,
-    const Mapping<2, 2> &                              mapping0,
-    const Mapping<2, 2> &                              mapping1,
-    const double                                       tol)
+  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)
   {
-    Assert(ReferenceCell::n_vertices_to_type(
-             2, mapping0.get_vertices(cell0).size()) ==
-             ReferenceCells::Quadrilateral,
-           ExcNotImplemented());
-    Assert(ReferenceCell::n_vertices_to_type(
-             2, mapping1.get_vertices(cell1).size()) ==
-             ReferenceCells::Quadrilateral,
-           ExcNotImplemented());
-
-    std::array<Point<2>, 4> vertices0, vertices1;
-    std::copy_n(mapping0.get_vertices(cell0).begin(), 4, vertices0.begin());
-    std::copy_n(mapping1.get_vertices(cell1).begin(), 4, vertices1.begin());
-
-    std::swap(vertices0[2], vertices0[3]);
-    std::swap(vertices1[2], vertices1[3]);
     const auto intersection_test =
       internal::compute_intersection(vertices0, vertices1);
 
@@ -435,33 +400,14 @@ namespace CGALWrappers
       }
   }
 
-
-
   // Specialization for quad \cap line
   template <>
   std::vector<std::array<Point<2>, 2>>
-  compute_intersection_of_cells<2, 1, 2>(
-    const typename Triangulation<2, 2>::cell_iterator &cell0,
-    const typename Triangulation<1, 2>::cell_iterator &cell1,
-    const Mapping<2, 2> &                              mapping0,
-    const Mapping<1, 2> &                              mapping1,
-    const double                                       tol)
+  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)
   {
-    Assert(ReferenceCell::n_vertices_to_type(
-             2, mapping0.get_vertices(cell0).size()) ==
-             ReferenceCells::Quadrilateral,
-           ExcNotImplemented());
-    Assert(ReferenceCell::n_vertices_to_type(
-             1, mapping1.get_vertices(cell1).size()) == ReferenceCells::Line,
-           ExcNotImplemented());
-
-    std::array<Point<2>, 4> vertices0;
-    std::array<Point<2>, 2> vertices1;
-    std::copy_n(mapping0.get_vertices(cell0).begin(), 4, vertices0.begin());
-    std::copy_n(mapping1.get_vertices(cell1).begin(), 2, vertices1.begin());
-
-    std::swap(vertices0[2], vertices0[3]);
-
     std::array<CGALPoint2, 4> pts;
     std::transform(
       vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<2> &p) {
@@ -496,34 +442,15 @@ namespace CGALWrappers
     return vertices;
   }
 
-
-
   // specialization for hex \cap line
   template <>
   std::vector<std::array<Point<3>, 2>>
-  compute_intersection_of_cells<3, 1, 3>(
-    const typename Triangulation<3, 3>::cell_iterator &cell0,
-    const typename Triangulation<1, 3>::cell_iterator &cell1,
-    const Mapping<3, 3> &                              mapping0,
-    const Mapping<1, 3> &                              mapping1,
-    const double                                       tol)
+  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)
   {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 1, 5)
-
-    Assert(ReferenceCell::n_vertices_to_type(
-             3, mapping0.get_vertices(cell0).size()) ==
-             ReferenceCells::Hexahedron,
-           ExcNotImplemented());
-    Assert(ReferenceCell::n_vertices_to_type(
-             1, mapping1.get_vertices(cell1).size()) == ReferenceCells::Line,
-           ExcNotImplemented());
-    (void)tol;
-
-    std::array<Point<3>, 8> vertices0; // 8 vertices of the hex
-    std::array<Point<3>, 2> vertices1; // 2 endpoints of the segment
-    std::copy_n(mapping0.get_vertices(cell0).begin(), 8, vertices0.begin());
-    std::copy_n(mapping1.get_vertices(cell1).begin(), 2, vertices1.begin());
-
     std::array<CGALPoint3_exact, 8> pts;
     std::transform(
       vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<3> &p) {
@@ -566,43 +493,21 @@ namespace CGALWrappers
       false,
       ExcMessage(
         "This function requires a version of CGAL greater or equal than 5.1.5."));
-    (void)cell0;
-    (void)cell1;
-    (void)mapping0;
-    (void)mapping1;
+    (void)vertices0;
+    (void)vertices1;
     (void)tol;
     return {};
-
 #  endif
   }
 
-
-
   template <>
   std::vector<std::array<Point<3>, 3>>
-  compute_intersection_of_cells<3, 2, 3>(
-    const typename Triangulation<3, 3>::cell_iterator &cell0,
-    const typename Triangulation<2, 3>::cell_iterator &cell1,
-    const Mapping<3, 3> &                              mapping0,
-    const Mapping<2, 3> &                              mapping1,
-    const double                                       tol)
+  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)
   {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 1, 5)
-
-    Assert(ReferenceCell::n_vertices_to_type(
-             3, mapping0.get_vertices(cell0).size()) ==
-             ReferenceCells::Hexahedron,
-           ExcNotImplemented());
-    Assert(ReferenceCell::n_vertices_to_type(
-             2, mapping1.get_vertices(cell1).size()) ==
-             ReferenceCells::Quadrilateral,
-           ExcNotImplemented());
-
-    std::array<Point<3>, 8> vertices0; // 8 vertices of the hex
-    std::array<Point<3>, 4> vertices1; // 4 vertices of the quad
-    std::copy_n(mapping0.get_vertices(cell0).begin(), 8, vertices0.begin());
-    std::copy_n(mapping1.get_vertices(cell1).begin(), 4, vertices1.begin());
-
     std::array<CGALPoint3_exact, 8> pts_hex;
     std::array<CGALPoint3_exact, 4> pts_quad;
     std::transform(
@@ -621,7 +526,6 @@ namespace CGALWrappers
         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;
@@ -685,71 +589,96 @@ namespace CGALWrappers
       }
 
     return vertices;
-
 #  else
-
     Assert(
       false,
       ExcMessage(
         "This function requires a version of CGAL greater or equal than 5.1.5."));
-    (void)cell0;
-    (void)cell1;
-    (void)mapping0;
-    (void)mapping1;
+    (void)vertices0;
+    (void)vertices1;
     (void)tol;
     return {};
 #  endif
   }
 
-
-
   template <>
   std::vector<std::array<Point<3>, 4>>
-  compute_intersection_of_cells<3, 3, 3>(
-    const typename Triangulation<3, 3>::cell_iterator &cell0,
-    const typename Triangulation<3, 3>::cell_iterator &cell1,
-    const Mapping<3, 3> &                              mapping0,
-    const Mapping<3, 3> &                              mapping1,
-    const double                                       tol)
+  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)
   {
-    Assert(ReferenceCell::n_vertices_to_type(
-             3, mapping0.get_vertices(cell0).size()) ==
-             ReferenceCells::Hexahedron,
+    // Surface_mesh surf0, surf1, sm;
+    // CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surf0);
+    // CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surf1);
+    // CGAL::Polygon_mesh_processing::triangulate_faces(surf0);
+    // CGAL::Polygon_mesh_processing::triangulate_faces(surf1);
+    // CGALWrappers::compute_boolean_operation(
+    //   surf0, surf1, CGALWrappers::BooleanOperation::compute_intersection,
+    //   sm);
+    // std::vector<std::array<Point<3>, 4>> vertices;
+    // if (CGAL::Polygon_mesh_processing::volume(sm) > tol)
+    //   {
+    //     // 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))}});
+    //       }
+    //     return vertices;
+    //   }
+    // else
+    //   {
+    //     return vertices;
+    //   }
+
+    // TODO: implememnt 3d/3d cut
+    AssertThrow(false, ExcMessage("Not yet implemented"));
+    (void)vertices0;
+    (void)vertices1;
+    (void)tol;
+    return {};
+  }
+
+  template <int dim0, int dim1, int spacedim>
+  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)
+  {
+    Assert(mapping0.get_vertices(cell0).size() == std::pow(2, dim0),
            ExcNotImplemented());
-    Assert(ReferenceCell::n_vertices_to_type(
-             3, mapping1.get_vertices(cell1).size()) ==
-             ReferenceCells::Hexahedron,
+    Assert(mapping1.get_vertices(cell1).size() == std::pow(2, dim1),
            ExcNotImplemented());
-    Surface_mesh surf0, surf1, sm;
-    CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surf0);
-    CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surf1);
-    CGAL::Polygon_mesh_processing::triangulate_faces(surf0);
-    CGAL::Polygon_mesh_processing::triangulate_faces(surf1);
-    CGALWrappers::compute_boolean_operation(
-      surf0, surf1, CGALWrappers::BooleanOperation::compute_intersection, sm);
-    std::vector<std::array<Point<3>, 4>> vertices;
-    if (CGAL::Polygon_mesh_processing::volume(sm) > tol)
-      {
-        // 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))}});
-          }
-        return vertices;
-      }
-    else
-      {
-        return vertices;
-      }
+
+    const auto vertices0 =
+      CGALWrappers::get_vertices_in_cgal_order<int(std::pow(2, dim0))>(
+        cell0, mapping0);
+
+    const auto vertices1 =
+      CGALWrappers::get_vertices_in_cgal_order<int(std::pow(2, dim1))>(
+        cell1, mapping1);
+
+    return compute_intersection_of_cells<dim0,
+                                         dim1,
+                                         spacedim,
+                                         int(std::pow(2, dim0)),
+                                         int(std::pow(2, dim1))>(vertices0,
+                                                                 vertices1,
+                                                                 tol);
   }
 
+#  include "intersections.inst"
+
 } // namespace CGALWrappers
 
 DEAL_II_NAMESPACE_CLOSE
@@ -758,6 +687,23 @@ DEAL_II_NAMESPACE_CLOSE
 
 DEAL_II_NAMESPACE_OPEN
 
+template <int dim0,
+          int dim1,
+          int spacedim,
+          int n_components0,
+          int n_components1>
+std::vector<std::array<Point<spacedim>, dim1 + 1>>
+compute_intersection_of_cells(
+  const std::array<Point<spacedim>, n_components0> &vertices0,
+  const std::array<Point<spacedim>, n_components1> &vertices1,
+  const double                                      tol)
+{
+  (void)vertices0;
+  (void)vertices1;
+  (void)tol;
+  AssertThrow(false, ExcNeedsCGAL());
+}
+
 template <int dim0, int dim1, int spacedim>
 std::vector<std::array<Point<spacedim>, dim1 + 1>>
 compute_intersection_of_cells(

diff --git a/source/cgal/intersections.inst.in b/source/cgal/intersections.inst.in
new file mode 100644 (file)
index 0000000..d550252
--- /dev/null
+++ b/source/cgal/intersections.inst.in
@@ -0,0 +1,31 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS)
+  {
+#if dim0 <= spacedim && dim0 >= dim1
+
+    template std::vector<std::array<Point<spacedim>, dim1 + 1>>
+    compute_intersection_of_cells<dim0, dim1, spacedim>(
+      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);
+
+#endif
+  }