From 3e2491b0c2635178d473f42d4286321fb0cf1c89 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Tue, 5 Jul 2022 18:55:37 +0200 Subject: [PATCH] Clean cgal utilities --- include/deal.II/base/quadrature_lib.h | 17 ++++++ include/deal.II/cgal/utilities.h | 83 ++++++--------------------- source/base/quadrature_lib.cc | 53 +++++++++++++++++ tests/cgal/cgal_quadrature_01.cc | 18 ++---- 4 files changed, 92 insertions(+), 79 deletions(-) diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index c4dd2edc20..5284bf849c 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -642,6 +642,23 @@ public: Quadrature compute_affine_transformation( const std::array, 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 + Quadrature + mapped_quadrature( + const std::vector, dim + 1>> &simplices) const; }; /** diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 0bb57f6580..93ecacacae 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -410,66 +410,22 @@ namespace CGALWrappers constexpr int spacedim = CGALTriangulationType::Point::Ambient_dimension::value; - QGaussSimplex quad(degree); - std::vector> pts; - std::vector wts; - std::array, spacedim + 1> vertices; // tets + std::vector, 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( - 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, 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(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(f->vertex(i)->point()); } + + vec_of_simplices.push_back(simplex); } - else - { - Assert(false, ExcMessage("Not a valid CGAL Triangulation.")); - } - return Quadrature(pts, wts); + + return QGaussSimplex(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; + 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); - 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>; - 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_triangulation(out_surface, tria); + CGALTriangulation tria; + tria.insert(out_surface.points().begin(), out_surface.points().end()); + return compute_quadrature(tria, degree); } } diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index bfe1d52b23..7803cb6196 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1260,6 +1260,36 @@ QSimplex::compute_affine_transformation( +template +template +Quadrature +QSimplex::mapped_quadrature( + const std::vector, 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> qp; + std::vector 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 &p) { return p; }); + std::transform(rule.get_weights().begin(), + rule.get_weights().end(), + std::back_inserter(ws), + [&](const double w) { return w; }); + } + return Quadrature(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, 2 + 1> &vertices) const; } // namespace dealii + +namespace dealii +{ + template Quadrature<2> + QSimplex<1>::mapped_quadrature( + const std::vector, 1 + 1>> &simplices) const; + + template Quadrature<3> + QSimplex<1>::mapped_quadrature( + const std::vector, 1 + 1>> &simplices) const; + + template Quadrature<2> + QSimplex<2>::mapped_quadrature( + const std::vector, 2 + 1>> &simplices) const; + + template Quadrature<3> + QSimplex<2>::mapped_quadrature( + const std::vector, 2 + 1>> &simplices) const; + + template Quadrature<3> + QSimplex<3>::mapped_quadrature( + const std::vector, 3 + 1>> &simplices) const; +} // namespace dealii diff --git a/tests/cgal/cgal_quadrature_01.cc b/tests/cgal/cgal_quadrature_01.cc index e20674cd68..71df579487 100644 --- a/tests/cgal/cgal_quadrature_01.cc +++ b/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; +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using CGALPoint = CGAL::Point_3; +using CGALTriangulation = CGAL::Triangulation_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 @@ -55,13 +47,13 @@ test() SOURCE_DIR "/input_grids/octahedron.off"}; CGAL::Surface_mesh 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(), -- 2.39.5