Merge pull request #14097 from luca-heltai/cgal-clean_utilities

Clean cgal utilities

include/deal.II/base/quadrature_lib.h

@@ -642,6 +642,23 @@ class QSimplex : public Quadrature<dim>
   Quadrature<spacedim>
   compute_affine_transformation(
     const std::array<Point<spacedim>, dim + 1> &vertices) const;
+
+  /**
+   *
+   * Given a partition of a cell into simplices, this function creates a
+   * quadrature rule on the cell by collecting Quadrature objects on each
+   * simplex. A simplex is identified by its vertices, which are stored into an
+   * array of Points. Hence, this function can provide quadrature rules on
+   * polygons (or polyhedra), as they can be split into simplices.
+   *
+   *
+   * @param simplices A std::vector where each entry is an array of `dim+1` points, which identifies the vertices of a simplex.
+   * @return A Quadrature object on the cell.
+   */
+  template <int spacedim = dim>
+  Quadrature<spacedim>
+  mapped_quadrature(
+    const std::vector<std::array<Point<spacedim>, dim + 1>> &simplices) const;
 };
 
 /**

include/deal.II/cgal/utilities.h

@@ -410,66 +410,22 @@ namespace CGALWrappers
 
     constexpr int spacedim =
       CGALTriangulationType::Point::Ambient_dimension::value;
-    QGaussSimplex<spacedim>                           quad(degree);
-    std::vector<dealii::Point<spacedim>>              pts;
-    std::vector<double>                               wts;
-    std::array<dealii::Point<spacedim>, spacedim + 1> vertices; // tets
+    std::vector<std::array<dealii::Point<spacedim>, spacedim + 1>>
+      vec_of_simplices; // tets
 
-    const auto is_c3t3 = boost::hana::is_valid(
-      [](auto &&obj) -> decltype(obj.cells_in_complex_begin()) {});
-
-    const auto is_tria3 = boost::hana::is_valid(
-      [](auto &&obj) -> decltype(obj.finite_cell_handles()) {});
-
-    if constexpr (decltype(is_c3t3(tria)){})
-      {
-        for (auto iter = tria.cells_in_complex_begin();
-             iter != tria.cells_in_complex_end();
-             ++iter)
-          {
-            for (unsigned int i = 0; i < (spacedim + 1); ++i)
-              {
-                vertices[i] = cgal_point_to_dealii_point<spacedim>(
-                  iter->vertex(i)->point());
-              }
-
-            auto local_quad = quad.compute_affine_transformation(vertices);
-            std::transform(local_quad.get_points().begin(),
-                           local_quad.get_points().end(),
-                           std::back_inserter(pts),
-                           [&pts](const auto &p) { return p; });
-            std::transform(local_quad.get_weights().begin(),
-                           local_quad.get_weights().end(),
-                           std::back_inserter(wts),
-                           [&wts](const double w) { return w; });
-          }
-      }
-    else if constexpr (decltype(is_tria3(tria)){})
+    std::array<dealii::Point<spacedim>, spacedim + 1> simplex;
+    for (const auto &f : tria.finite_cell_handles())
       {
-        for (const auto &f : tria.finite_cell_handles())
+        for (unsigned int i = 0; i < (spacedim + 1); ++i)
           {
-            for (unsigned int i = 0; i < (spacedim + 1); ++i)
-              {
-                vertices[i] =
-                  cgal_point_to_dealii_point<spacedim>(f->vertex(i)->point());
-              }
-
-            auto local_quad = quad.compute_affine_transformation(vertices);
-            std::transform(local_quad.get_points().begin(),
-                           local_quad.get_points().end(),
-                           std::back_inserter(pts),
-                           [&pts](const auto &p) { return p; });
-            std::transform(local_quad.get_weights().begin(),
-                           local_quad.get_weights().end(),
-                           std::back_inserter(wts),
-                           [&wts](const double w) { return w; });
+            simplex[i] =
+              cgal_point_to_dealii_point<spacedim>(f->vertex(i)->point());
           }
+
+        vec_of_simplices.push_back(simplex);
       }
-    else
-      {
-        Assert(false, ExcMessage("Not a valid CGAL Triangulation."));
-      }
-    return Quadrature<spacedim>(pts, wts);
+
+    return QGaussSimplex<spacedim>(degree).mapped_quadrature(vec_of_simplices);
   }
 
 
@@ -506,26 +462,21 @@ namespace CGALWrappers
       {
         using K         = CGAL::Exact_predicates_inexact_constructions_kernel;
         using CGALPoint = CGAL::Point_3<K>;
+        using CGALTriangulation = CGAL::Triangulation_3<K>;
         CGAL::Surface_mesh<CGALPoint> surface_1, surface_2, out_surface;
         dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surface_1);
         dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surface_2);
         // They have to be triangle meshes
         CGAL::Polygon_mesh_processing::triangulate_faces(surface_1);
         CGAL::Polygon_mesh_processing::triangulate_faces(surface_2);
-        Assert(CGAL::is_triangle_mesh(surface_1),
+        Assert(CGAL::is_triangle_mesh(surface_1) &&
+                 CGAL::is_triangle_mesh(surface_2),
                ExcMessage("The surface must be a triangle mesh."));
         compute_boolean_operation(surface_1, surface_2, bool_op, out_surface);
 
-        using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<
-          K,
-          CGAL::Surface_mesh<CGALPoint>>;
-        using Tr            = typename CGAL::Mesh_triangulation_3<Mesh_domain,
-                                                       CGAL::Default,
-                                                       ConcurrencyTag>::type;
-        using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
-        using C3t3          = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
-        C3t3 tria;
-        cgal_surface_mesh_to_cgal_triangulation(out_surface, tria);
+        CGALTriangulation tria;
+        tria.insert(out_surface.points().begin(), out_surface.points().end());
+
         return compute_quadrature(tria, degree);
       }
   }

source/base/quadrature_lib.cc

@@ -1260,6 +1260,36 @@ QSimplex<dim>::compute_affine_transformation(
 
 
 
+template <int dim>
+template <int spacedim>
+Quadrature<spacedim>
+QSimplex<dim>::mapped_quadrature(
+  const std::vector<std::array<Point<spacedim>, dim + 1>> &simplices) const
+{
+  Assert(!(dim == 1 && spacedim == 1),
+         ExcMessage("This function is not supposed to work in 1D-1D case."));
+  Assert(dim <= spacedim,
+         ExcMessage("Invalid combination of dim and spacedim ."));
+
+  std::vector<Point<spacedim>> qp;
+  std::vector<double>          ws;
+  for (const auto &simplex : simplices)
+    {
+      const auto rule = this->compute_affine_transformation(simplex);
+      std::transform(rule.get_points().begin(),
+                     rule.get_points().end(),
+                     std::back_inserter(qp),
+                     [&](const Point<spacedim> &p) { return p; });
+      std::transform(rule.get_weights().begin(),
+                     rule.get_weights().end(),
+                     std::back_inserter(ws),
+                     [&](const double w) { return w; });
+    }
+  return Quadrature<spacedim>(qp, ws);
+}
+
+
+
 QTrianglePolar::QTrianglePolar(const Quadrature<1> &radial_quadrature,
                                const Quadrature<1> &angular_quadrature)
   : QSimplex<2>(Quadrature<2>())
@@ -2227,3 +2257,26 @@ namespace dealii
   QSimplex<2>::compute_affine_transformation(
     const std::array<Point<3>, 2 + 1> &vertices) const;
 } // namespace dealii
+
+namespace dealii
+{
+  template Quadrature<2>
+  QSimplex<1>::mapped_quadrature(
+    const std::vector<std::array<Point<2>, 1 + 1>> &simplices) const;
+
+  template Quadrature<3>
+  QSimplex<1>::mapped_quadrature(
+    const std::vector<std::array<Point<3>, 1 + 1>> &simplices) const;
+
+  template Quadrature<2>
+  QSimplex<2>::mapped_quadrature(
+    const std::vector<std::array<Point<2>, 2 + 1>> &simplices) const;
+
+  template Quadrature<3>
+  QSimplex<2>::mapped_quadrature(
+    const std::vector<std::array<Point<3>, 2 + 1>> &simplices) const;
+
+  template Quadrature<3>
+  QSimplex<3>::mapped_quadrature(
+    const std::vector<std::array<Point<3>, 3 + 1>> &simplices) const;
+} // namespace dealii

tests/cgal/cgal_quadrature_01.cc

@@ -30,19 +30,11 @@
 
 // Compute volumes by splitting polyhedra into simplices.
 
-using K         = CGAL::Exact_predicates_inexact_constructions_kernel;
-using CGALPoint = CGAL::Point_3<K>;
+using K                 = CGAL::Exact_predicates_inexact_constructions_kernel;
+using CGALPoint         = CGAL::Point_3<K>;
+using CGALTriangulation = CGAL::Triangulation_3<K>;
 using namespace CGALWrappers;
-using Mesh_domain =
-  CGAL::Polyhedral_mesh_domain_with_features_3<K,
-                                               CGAL::Surface_mesh<CGALPoint>>;
-using Tr = typename CGAL::
-  Mesh_triangulation_3<Mesh_domain, CGAL::Default, ConcurrencyTag>::type;
 
-using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
-using C3t3          = CGAL::Mesh_complex_3_in_triangulation_3<Tr,
-                                                     Mesh_domain::Corner_index,
-                                                     Mesh_domain::Curve_index>;
 
 
 void
@@ -55,13 +47,13 @@ test()
                                         SOURCE_DIR
                                         "/input_grids/octahedron.off"};
   CGAL::Surface_mesh<CGALPoint>  sm;
-  C3t3                           tria;
+  CGALTriangulation              tria;
   constexpr int                  degree = 3;
   for (const auto &name : fnames)
     {
       std::ifstream input(name);
       input >> sm;
-      cgal_surface_mesh_to_cgal_triangulation(sm, tria);
+      tria.insert(sm.points().begin(), sm.points().end());
       auto b = compute_quadrature(tria, degree);
       deallog << "Volume of poly with Quadrature: " << std::setprecision(12)
               << std::accumulate(b.get_weights().begin(),