From: Marco Feder Date: Sat, 21 May 2022 19:37:09 +0000 (+0200) Subject: Compute Quadrature formula over a general poly X-Git-Tag: v9.4.0-rc1~103^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=198c9dd9c68c12f909952bd5653bd3e5240bb1d3;p=dealii.git Compute Quadrature formula over a general poly --- diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 4c37536268..d6025b78eb 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -19,6 +19,12 @@ #include #ifdef DEAL_II_WITH_CGAL +# include + +# include + +# include + # include # include # include @@ -28,13 +34,16 @@ # include # include # include +# include # include # include # include # include # include +# include # include # include +# include # include @@ -198,6 +207,75 @@ namespace CGALWrappers const CGAL::Surface_mesh &surface_mesh_2, const BooleanOperation & boolean_operation, CGAL::Surface_mesh & 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 + dealii::Quadrature + 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 The global quadrature rule on the polygon/polyhedron. + */ + template + dealii::Quadrature + compute_quadrature_on_boolean_operation( + const typename dealii::Triangulation::cell_iterator &cell0, + const typename dealii::Triangulation::cell_iterator &cell1, + const unsigned int degree, + const BooleanOperation & bool_op, + const Mapping &mapping0 = + (dealii::ReferenceCells::get_hypercube() + .template get_default_linear_mapping()), + const Mapping &mapping1 = + (dealii::ReferenceCells::get_hypercube() + .template get_default_linear_mapping())); + + /** + * 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 The global quadrature rule on the polygon/polyhedron. + */ + template + dealii::Quadrature + compute_quadrature_on_intersection( + const typename dealii::Triangulation::cell_iterator &cell0, + const typename dealii::Triangulation::cell_iterator &cell1, + const unsigned int degree, + const Mapping &mapping0, + const Mapping &mapping1); } // namespace CGALWrappers # ifndef DOXYGEN @@ -330,6 +408,178 @@ namespace CGALWrappers Assert(res, ExcMessage("The boolean operation was not successfully computed.")); } + + + + template + dealii::Quadrature + 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 quad(degree); + std::vector> pts; + std::vector wts; + std::array, 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( + 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(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(pts, wts); + } + + + + template + dealii::Quadrature + compute_quadrature_on_boolean_operation( + const typename dealii::Triangulation::cell_iterator &cell0, + const typename dealii::Triangulation::cell_iterator &cell1, + const unsigned int degree, + const BooleanOperation & bool_op, + const Mapping &mapping0, + const Mapping &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; + CGAL::Surface_mesh 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>; + using Tr = typename CGAL::Mesh_triangulation_3::type; + using Mesh_criteria = CGAL::Mesh_criteria_3; + using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3; + C3t3 tria; + cgal_surface_mesh_to_cgal_coarse_triangulation(out_surface, tria); + return compute_quadrature(tria, degree); + } + } + + + + template + dealii::Quadrature + compute_quadrature_on_intersection( + const typename dealii::Triangulation::cell_iterator &cell0, + const typename dealii::Triangulation::cell_iterator &cell1, + const unsigned int degree, + const Mapping &mapping0, + const Mapping &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; + using CGALTriangulation = CGAL::Triangulation_3; + + CGAL::Surface_mesh 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 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 diff --git a/tests/cgal/cgal_quadrature_01.cc b/tests/cgal/cgal_quadrature_01.cc new file mode 100644 index 0000000000..3f0ac6472b --- /dev/null +++ b/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 + +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +// Compute volumes by splitting polyhedra into simplices. + +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using CGALPoint = CGAL::Point_3; +using namespace CGALWrappers; +using Mesh_domain = + CGAL::Polyhedral_mesh_domain_with_features_3>; +using Tr = typename CGAL:: + Mesh_triangulation_3::type; + +using Mesh_criteria = CGAL::Mesh_criteria_3; +using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3; + + +void +test() +{ + const std::vector fnames{"input_grids/cube.off", + "input_grids/tetrahedron.off", + "input_grids/hedra.off", + "input_grids/octahedron.off"}; + CGAL::Surface_mesh 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(); +} diff --git a/tests/cgal/cgal_quadrature_01.output b/tests/cgal/cgal_quadrature_01.output new file mode 100644 index 0000000000..f9827d82ce --- /dev/null +++ b/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 diff --git a/tests/cgal/cgal_quadrature_02.cc b/tests/cgal/cgal_quadrature_02.cc new file mode 100644 index 0000000000..f20b133cb1 --- /dev/null +++ b/tests/cgal/cgal_quadrature_02.cc @@ -0,0 +1,86 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#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> 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(); +} diff --git a/tests/cgal/cgal_quadrature_02.output b/tests/cgal/cgal_quadrature_02.output new file mode 100644 index 0000000000..d90cdd607a --- /dev/null +++ b/tests/cgal/cgal_quadrature_02.output @@ -0,0 +1,5 @@ + +DEAL::Volume of poly with Quadrature: 0.430666666667 +DEAL::Volume of poly with Quadrature: 0.166166666667 +DEAL::Volume of poly with Quadrature: 0.121500000000 +DEAL::Volume of poly with Quadrature: 0.972000000000 diff --git a/tests/cgal/input_grids/octahedron.off b/tests/cgal/input_grids/octahedron.off new file mode 100644 index 0000000000..1b38ebfc59 --- /dev/null +++ b/tests/cgal/input_grids/octahedron.off @@ -0,0 +1,16 @@ +OFF +6 8 12 +0.0 0.0 2.0 +2.000000 0.000000 0.000000 +0.000000 2.000000 0.000000 +-2.000000 0.000000 0.000000 +0.000000 -2.000000 0.000000 +0.0 0.0 -2.0 +3 1 0 4 +3 4 0 3 +3 3 0 2 +3 2 0 1 +3 1 5 2 +3 2 5 3 +3 3 5 4 +3 4 5 1 \ No newline at end of file