Compute Quadrature formula over a general poly

include/deal.II/cgal/utilities.h

@@ -19,6 +19,12 @@
 #include <deal.II/base/config.h>
 
 #ifdef DEAL_II_WITH_CGAL
+#  include <deal.II/base/quadrature_lib.h>
+
+#  include <deal.II/grid/tria.h>
+
+#  include <boost/hana.hpp>
+
 #  include <CGAL/Cartesian.h>
 #  include <CGAL/Complex_2_in_triangulation_3.h>
 #  include <CGAL/Exact_predicates_exact_constructions_kernel.h>
@@ -28,13 +34,16 @@
 #  include <CGAL/Mesh_criteria_3.h>
 #  include <CGAL/Mesh_triangulation_3.h>
 #  include <CGAL/Polygon_mesh_processing/corefinement.h>
+#  include <CGAL/Polygon_mesh_processing/measure.h>
 #  include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
 #  include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
 #  include <CGAL/Simple_cartesian.h>
 #  include <CGAL/Surface_mesh.h>
 #  include <CGAL/Triangulation_3.h>
+#  include <CGAL/convex_hull_3.h>
 #  include <CGAL/make_mesh_3.h>
 #  include <CGAL/make_surface_mesh.h>
+#  include <deal.II/cgal/surface_mesh.h>
 
 #  include <type_traits>
 
@@ -198,6 +207,73 @@ namespace CGALWrappers
     const CGAL::Surface_mesh<CGALPointType> &surface_mesh_2,
     const BooleanOperation &                 boolean_operation,
     CGAL::Surface_mesh<CGALPointType> &      output_surface_mesh);
+
+  /**
+   * Given a CGAL Triangulation describing a polyhedral region, create
+   * a Quadrature rule to integrate over the polygon by looping trough all the
+   * vertices and exploiting QGaussSimplex.
+   *
+   * @param[in] tria The CGAL triangulation object describing the polyhedral
+   * region.
+   * @param[in] degree Desired degree of the Quadrature rule on each element of
+   * the polyhedral.
+   * @return A global Quadrature rule that allows to integrate over the polyhedron.
+   */
+  template <typename Tr>
+  dealii::Quadrature<Tr::Point::Ambient_dimension::value>
+  compute_quadrature(const Tr &tria, const unsigned int degree);
+
+  /**
+   * Compute a Quadrature formula over the polygonal/polyhedral region described
+   * by a BooleanOperation between two deal.II cell_iterator. The degree of the
+   * Quadrature
+   * rule is specified by the argument @p degree. This function splits the polygon/polyhedron
+   * into simplices and collects QGaussSimplex quadrature rules.
+   *
+   * @param [in] cell0 A cell_iterator to the first deal.II cell.
+   * @param [in] cell1 A cell_iterator to the second deal.II cell.
+   * @param [in] degree The degree of accuracy you wish to get for the global quadrature formula.
+   * @param [in] bool_op The BooleanOperation to be performed.
+   * @param [in] mapping0 Mapping object for the first cell.
+   * @param [in] mapping1 Mapping object for the first cell.
+   * @return [out] Quadrature<spacedim> The global quadrature rule on the polygon/polyhedron.
+   */
+  template <int dim0, int dim1, int spacedim>
+  dealii::Quadrature<spacedim>
+  compute_quadrature_on_boolean_operation(
+    const typename dealii::Triangulation<dim0, spacedim>::cell_iterator &cell0,
+    const typename dealii::Triangulation<dim1, spacedim>::cell_iterator &cell1,
+    const unsigned int                                                   degree,
+    const BooleanOperation &       bool_op,
+    const Mapping<dim0, spacedim> &mapping0 =
+      (dealii::ReferenceCells::get_hypercube<dim0>()
+         .template get_default_linear_mapping<dim0, spacedim>()),
+    const Mapping<dim1, spacedim> &mapping1 =
+      (dealii::ReferenceCells::get_hypercube<dim1>()
+         .template get_default_linear_mapping<dim1, spacedim>()));
+  /**
+   * A specialization of the function above when the BooleanOperation is an
+   * intersection. The rationale behind this specialization is that deal.II
+   * cells are convex sets, and as the intersection of convex sets is itself
+   * convex, this function internally exploits this to use a cheaper way to mesh
+   * the inside.
+   *
+   *
+   * @param [in] cell0 A cell_iterator to the first deal.II cell.
+   * @param [in] cell1 A cell_iterator to the second deal.II cell.
+   * @param [in] mapping0 Mapping object for the first cell.
+   * @param [in] mapping1 Mapping object for the first cell.
+   * @param [in] degree The degree of accuracy you wish to get for the global quadrature formula.
+   * @return [out] Quadrature<spacedim> The global quadrature rule on the polygon/polyhedron.
+   */
+  template <int dim0, int dim1, int spacedim>
+  dealii::Quadrature<spacedim>
+  compute_quadrature_on_intersection(
+    const typename dealii::Triangulation<dim0, spacedim>::cell_iterator &cell0,
+    const typename dealii::Triangulation<dim1, spacedim>::cell_iterator &cell1,
+    const unsigned int                                                   degree,
+    const Mapping<dim0, spacedim> &mapping0,
+    const Mapping<dim1, spacedim> &mapping1);
 } // namespace CGALWrappers
 
 #  ifndef DOXYGEN
@@ -330,6 +406,175 @@ namespace CGALWrappers
     Assert(res,
            ExcMessage("The boolean operation was not successfully computed."));
   }
+
+
+
+  template <typename Tr>
+  dealii::Quadrature<Tr::Point::Ambient_dimension::value>
+  compute_quadrature(const Tr &tria, const unsigned int degree)
+  {
+    Assert(tria.is_valid(), ExcMessage("The triangulation is not valid."));
+    Assert(Tr::Point::Ambient_dimension::value == 3, ExcNotImplemented());
+    Assert(degree > 0,
+           ExcMessage("The degree of the Quadrature formula is not positive."));
+
+    constexpr int           spacedim = Tr::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
+
+    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)){})
+      {
+        for (const auto &f : tria.finite_cell_handles())
+          {
+            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; });
+          }
+      }
+    else
+      {
+        Assert(false, ExcMessage("Not a valid CGAL Triangulation."));
+      }
+    return Quadrature<spacedim>(pts, wts);
+  }
+
+
+
+  template <int dim0, int dim1, int spacedim>
+  dealii::Quadrature<spacedim>
+  compute_quadrature_on_boolean_operation(
+    const typename dealii::Triangulation<dim0, spacedim>::cell_iterator &cell0,
+    const typename dealii::Triangulation<dim1, spacedim>::cell_iterator &cell1,
+    const unsigned int                                                   degree,
+    const BooleanOperation &       bool_op,
+    const Mapping<dim0, spacedim> &mapping0,
+    const Mapping<dim1, spacedim> &mapping1)
+  {
+    Assert(dim0 == 3 && dim1 == 3 && spacedim == 3,
+           ExcNotImplemented("2D geometries are not yet supported."));
+    if (dim1 > dim0)
+      {
+        return compute_quadrature_on_boolean_operation(
+          cell1,
+          cell0,
+          degree,
+          bool_op,
+          mapping1,
+          mapping0); // This function works for dim1<=dim0, so swap them if this
+                     // is not the case.
+      }
+    if (bool_op == BooleanOperation::compute_intersection)
+      {
+        return compute_quadrature_on_intersection(
+          cell0, cell1, degree, mapping0, mapping1);
+      }
+    else
+      {
+        using K         = CGAL::Exact_predicates_inexact_constructions_kernel;
+        using CGALPoint = CGAL::Point_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),
+               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_coarse_triangulation(out_surface, tria);
+        return compute_quadrature(tria, degree);
+      }
+  }
+
+
+
+  template <int dim0, int dim1, int spacedim>
+  dealii::Quadrature<spacedim>
+  compute_quadrature_on_intersection(
+    const typename dealii::Triangulation<dim0, spacedim>::cell_iterator &cell0,
+    const typename dealii::Triangulation<dim1, spacedim>::cell_iterator &cell1,
+    const unsigned int                                                   degree,
+    const Mapping<dim0, spacedim> &mapping0,
+    const Mapping<dim1, spacedim> &mapping1)
+  {
+    Assert(dim0 == 3 && dim1 == 3 && spacedim == 3,
+           ExcNotImplemented("2D geometries are not yet supported."));
+    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);
+
+    compute_boolean_operation(surface_1,
+                              surface_2,
+                              BooleanOperation::compute_intersection,
+                              out_surface);
+    CGAL::Surface_mesh<CGALPoint> dummy;
+    CGALTriangulation             tr;
+    CGAL::convex_hull_3(out_surface.points().begin(),
+                        out_surface.points().end(),
+                        dummy);
+    tr.insert(dummy.points().begin(), dummy.points().end());
+    return compute_quadrature(tr, degree);
+  }
 } // namespace CGALWrappers
 #  endif
 

tests/cgal/cgal_quadrature_01.cc

@@ -0,0 +1,81 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Read a Surface_mesh from an .off file, then create a coarse mesh filled with
+// tets
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <CGAL/IO/File_medit.h>
+#include <CGAL/IO/io.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
+#include <deal.II/cgal/utilities.h>
+
+#include "../tests.h"
+
+// Compute volumes by splitting polyhedra into simplices.
+
+using K         = CGAL::Exact_predicates_inexact_constructions_kernel;
+using CGALPoint = CGAL::Point_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
+test()
+{
+  const std::vector<std::string> fnames{"input_grids/cube.off",
+                                        "input_grids/tetrahedron.off",
+                                        "input_grids/hedra.off",
+                                        "input_grids/octahedron.off"};
+  CGAL::Surface_mesh<CGALPoint>  sm;
+  C3t3                           tria;
+  constexpr int                  degree = 3;
+  for (const auto &name : fnames)
+    {
+      std::ifstream input(name);
+      input >> sm;
+      cgal_surface_mesh_to_cgal_coarse_triangulation(sm, tria);
+      auto b = compute_quadrature(tria, degree);
+      deallog << "Volume of poly with Quadrature: " << std::setprecision(12)
+              << std::accumulate(b.get_weights().begin(),
+                                 b.get_weights().end(),
+                                 0.)
+              << "\t Expected:" << std::setprecision(12)
+              << CGAL::to_double(CGAL::Polygon_mesh_processing::volume(sm))
+              << std::endl;
+      sm.clear(); // reset surface
+      tria.clear();
+    }
+}
+
+int
+main()
+{
+  initlog();
+  test();
+}

tests/cgal/cgal_quadrature_01.output

@@ -0,0 +1,5 @@
+
+DEAL::Volume of poly with Quadrature: 8.00000000000	 Expected:8.00000000000
+DEAL::Volume of poly with Quadrature: 0.166666666667	 Expected:0.166666666667
+DEAL::Volume of poly with Quadrature: 6.38820454400	 Expected:6.38820454400
+DEAL::Volume of poly with Quadrature: 10.6666666667	 Expected:10.6666666667