#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>
# 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>
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 CGALTriangulationType>
+ dealii::Quadrature<CGALTriangulationType::Point::Ambient_dimension::value>
+ compute_quadrature(const CGALTriangulationType &tria,
+ const unsigned int degree);
+
+ /**
+ * Compute a Quadrature formula over the polygonal/polyhedral region described
+ * by a BooleanOperation between two deal.II cell_iterator objects. 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
+ * affine 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
Assert(res,
ExcMessage("The boolean operation was not successfully computed."));
}
+
+
+
+ template <typename CGALTriangulationType>
+ dealii::Quadrature<CGALTriangulationType::Point::Ambient_dimension::value>
+ compute_quadrature(const CGALTriangulationType &tria,
+ const unsigned int degree)
+ {
+ Assert(tria.is_valid(), ExcMessage("The triangulation is not valid."));
+ Assert(CGALTriangulationType::Point::Ambient_dimension::value == 3,
+ ExcNotImplemented());
+ Assert(degree > 0,
+ ExcMessage("The degree of the Quadrature formula is not positive."));
+
+ 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
+
+ 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Intersect reference cells and compute the volume of the intersection by
+// integration.
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.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 of intersections of polyhedra by splitting each intersection
+// into simplices.
+
+using namespace CGALWrappers;
+
+void
+test()
+{
+ using namespace ReferenceCells;
+ std::vector<std::pair<ReferenceCell, ReferenceCell>> ref_pairs = {
+ {Hexahedron, Pyramid},
+ {Hexahedron, Tetrahedron},
+ {Tetrahedron, Pyramid},
+ {Pyramid, Pyramid}};
+
+ constexpr int degree = 3;
+ constexpr auto bool_op = BooleanOperation::compute_intersection;
+ for (const auto &pair : ref_pairs)
+ {
+ const auto ref_cell0 = pair.first;
+ const auto ref_cell1 = pair.second;
+ Triangulation<3> tria0;
+ Triangulation<3> tria1;
+ GridGenerator::reference_cell(tria0, ref_cell0);
+ GridGenerator::reference_cell(tria1, ref_cell1);
+ const auto mapping0 = ref_cell0.template get_default_mapping<3>(1);
+ const auto mapping1 = ref_cell1.template get_default_mapping<3>(1);
+ const auto cell0 = tria0.begin_active();
+ const auto cell1 = tria1.begin_active();
+
+ // Shift each vertex by the same amount
+ const unsigned int n_vertices = cell1->n_vertices();
+ for (unsigned int i = 0; i < n_vertices; ++i)
+ {
+ cell1->vertex(i) += Point<3>(0.1, 0.1, 0.1);
+ }
+
+ auto test_quad = compute_quadrature_on_boolean_operation<3, 3, 3>(
+ cell0, cell1, degree, bool_op, *mapping0, *mapping1);
+ deallog << "Volume of poly with Quadrature: " << std::setprecision(12)
+ << std::accumulate(test_quad.get_weights().begin(),
+ test_quad.get_weights().end(),
+ 0.)
+ << std::endl;
+ }
+}
+
+int
+main()
+{
+ initlog();
+ test();
+}
projects / dealii.git / commitdiff