From: Marco Feder Date: Fri, 9 Sep 2022 07:57:01 +0000 (+0200) Subject: Add 3D/3D cuts with vertices only and minor fixes X-Git-Tag: v9.5.0-rc1~940^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ac5de68554de654dca261165ae013c19356fc559;p=dealii.git Add 3D/3D cuts with vertices only and minor fixes --- diff --git a/cmake/modules/FindCGAL.cmake b/cmake/modules/FindCGAL.cmake index 34692bee59..7071d5b31e 100644 --- a/cmake/modules/FindCGAL.cmake +++ b/cmake/modules/FindCGAL.cmake @@ -1,17 +1,17 @@ -# # --------------------------------------------------------------------- -# # -# # 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. -# # -# # --------------------------------------------------------------------- +## --------------------------------------------------------------------- +## +## 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. +## +## --------------------------------------------------------------------- # # Try to find the CGAL libraries diff --git a/doc/news/changes/minor/20220920Feder b/doc/news/changes/minor/20220920Feder new file mode 100644 index 0000000000..aaa411a047 --- /dev/null +++ b/doc/news/changes/minor/20220920Feder @@ -0,0 +1,6 @@ +New: The new function CGALWrappers::compute_intersection_of_cells computes +the intersection between two (affine) cells starting from the location of the +vertices. The intersection is described by a vector of arrays, where each array +identifies a simplex of the sub-tassellation of the intersection. +
+(Marco Feder, Johannes Heinz, 2022/09/20) diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 5284bf849c..0cfe1fa4ee 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -644,12 +644,12 @@ public: 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. + * Given a collection of simplices, this function creates a global + * quadrature rule on the area covered by the simplices, by mapping the + * current quadrature 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. diff --git a/include/deal.II/cgal/intersections.h b/include/deal.II/cgal/intersections.h index 8254d10e35..acf05fae36 100644 --- a/include/deal.II/cgal/intersections.h +++ b/include/deal.II/cgal/intersections.h @@ -28,26 +28,8 @@ DEAL_II_NAMESPACE_OPEN namespace CGALWrappers { - template - std::vector, dim1 + 1>> - compute_intersection_of_cells( - const std::array, n_components0> &vertices0, - const std::array, n_components1> &vertices1, - const double tol) - { - (void)vertices0; - (void)vertices1; - (void)tol; - Assert(false, ExcMessage("No explicit template instantiation available")); - return {}; - } - /** - * Given two deal.II affine cells, compute the intersection and return a + * Given two deal.II affine cells, compute their intersection and return a * vector of simplices, each one identified by an array of deal.II Points. All * the simplices together are a subdivision of the intersection. If cells are * non-affine, a geometrical error is introduced. If the @@ -75,6 +57,24 @@ namespace CGALWrappers const Mapping & mapping0, const Mapping & mapping1, const double tol = 1e-9); + + + /** + * Same function as above, but working directly with vertices. + */ + template + std::vector, dim1 + 1>> + compute_intersection_of_cells( + const std::array, n_vertices0> &vertices0, + const std::array, n_vertices1> &vertices1, + const double tol = 1e-9) + { + (void)vertices0; + (void)vertices1; + (void)tol; + Assert(false, ExcMessage("No explicit template instantiation available")); + return {}; + } } // namespace CGALWrappers DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 4e8bc23a80..10fe36523f 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -595,11 +595,11 @@ namespace CGALWrappers /** - * Get vertices of cell in CGAL ordering + * Get vertices of cell in CGAL ordering. * - * @param [in] cell A cell_iterator to a deal.II cell. - * @param [in] mapping Mapping object for the cell. - * @return [out] Array of vertices in CGAL order. + * @param cell A cell_iterator to a deal.II cell. + * @param mapping Mapping object for the cell. + * @return Array of vertices in CGAL order. */ template std::array, n_vertices> diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc index 7a859849b9..fbf5dd76b1 100644 --- a/source/cgal/intersections.cc +++ b/source/cgal/intersections.cc @@ -70,6 +70,7 @@ namespace CGALWrappers using CGALPoint2 = K::Point_2; using CGALPoint3 = K::Point_3; using CGALPoint3_exact = K_exact::Point_3; + using CGALPoint3_inexact = K_inexact::Point_3; using CGALSegment2 = K::Segment_2; using Surface_mesh = CGAL::Surface_mesh; using CGALSegment3 = K::Segment_3; @@ -141,6 +142,8 @@ namespace CGALWrappers } } + + void mark_domains(CDT &cdt) { @@ -200,6 +203,8 @@ namespace CGALWrappers return CGAL::intersection(triangle1, triangle2); } + + boost::optional> compute_intersection(const std::array, 3> &first_simplex, const std::array, 2> &second_simplex) @@ -227,6 +232,8 @@ namespace CGALWrappers return CGAL::intersection(segm, triangle); } + + // rectangle-rectangle std::vector compute_intersection(const std::array, 4> &first_simplex, @@ -255,6 +262,8 @@ namespace CGALWrappers return poly_list; } + + boost::optional> compute_intersection(const std::array, 2> &first_simplex, const std::array, 4> &second_simplex) @@ -335,6 +344,8 @@ namespace CGALWrappers } } // namespace internal + + // Specialization for quads template <> std::vector, 3>> @@ -400,6 +411,8 @@ namespace CGALWrappers } } + + // Specialization for quad \cap line template <> std::vector, 2>> @@ -601,51 +614,93 @@ namespace CGALWrappers # endif } + + template <> std::vector, 4>> compute_intersection_of_cells<3, 3, 3, 8, 8>( const std::array, 8> &vertices0, const std::array, 8> &vertices1, - const double tol) + + const double tol) { - // Surface_mesh surf0, surf1, sm; - // CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surf0); - // CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surf1); - // CGAL::Polygon_mesh_processing::triangulate_faces(surf0); - // CGAL::Polygon_mesh_processing::triangulate_faces(surf1); - // CGALWrappers::compute_boolean_operation( - // surf0, surf1, CGALWrappers::BooleanOperation::compute_intersection, - // sm); - // std::vector, 4>> vertices; - // if (CGAL::Polygon_mesh_processing::volume(sm) > tol) - // { - // // Collect tetrahedrons - // Triangulation3_inexact tria; - // tria.insert(sm.points().begin(), sm.points().end()); - // for (const auto &c : tria.finite_cell_handles()) - // { - // const auto &tet = tria.tetrahedron(c); - // vertices.push_back( - // {{CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(0)), - // CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(1)), - // CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(2)), - // CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(3))}}); - // } - // return vertices; - // } - // else - // { - // return vertices; - // } - - // TODO: implememnt 3d/3d cut - AssertThrow(false, ExcMessage("Not yet implemented")); - (void)vertices0; - (void)vertices1; - (void)tol; - return {}; + std::array pts_hex0; + std::array pts_hex1; + std::transform( + vertices0.begin(), + vertices0.end(), + pts_hex0.begin(), + [&](const Point<3> &p) { + return CGALWrappers::dealii_point_to_cgal_point(p); + }); + + std::transform( + vertices1.begin(), + vertices1.end(), + pts_hex1.begin(), + [&](const Point<3> &p) { + return CGALWrappers::dealii_point_to_cgal_point(p); + }); + + + Surface_mesh surf0, surf1, sm; + // Subdivide hex into tetrahedrons + std::vector, 4>> vertices; + Triangulation3_inexact tria0, tria1; + + tria0.insert(pts_hex0.begin(), pts_hex0.end()); + tria1.insert(pts_hex1.begin(), pts_hex1.end()); + + for (const auto &c0 : tria0.finite_cell_handles()) + { + const auto & tet0 = tria1.tetrahedron(c0); + [[maybe_unused]] const auto &tetg0 = + CGAL::make_tetrahedron(tet0.vertex(0), + tet0.vertex(1), + tet0.vertex(2), + tet0.vertex(3), + surf0); + for (const auto &c1 : tria1.finite_cell_handles()) + { + const auto & tet1 = tria1.tetrahedron(c1); + [[maybe_unused]] const auto &tetg1 = + CGAL::make_tetrahedron(tet1.vertex(0), + tet1.vertex(1), + tet1.vertex(2), + tet1.vertex(3), + surf1); + const bool test_intersection = + CGAL::PMP::corefine_and_compute_intersection(surf0, surf1, sm); + if (CGAL::Polygon_mesh_processing::volume(sm) > tol && + test_intersection) + { + // Collect tetrahedrons + Triangulation3_inexact tria; + tria.insert(sm.points().begin(), sm.points().end()); + for (const auto &c : tria.finite_cell_handles()) + { + const auto &tet = tria.tetrahedron(c); + vertices.push_back( + {{CGALWrappers::cgal_point_to_dealii_point<3>( + tet.vertex(0)), + CGALWrappers::cgal_point_to_dealii_point<3>( + tet.vertex(1)), + CGALWrappers::cgal_point_to_dealii_point<3>( + tet.vertex(2)), + CGALWrappers::cgal_point_to_dealii_point<3>( + tet.vertex(3))}}); + } + } + surf1.clear(); + sm.clear(); + } + surf0.clear(); + } + return vertices; } + + template std::vector, dim1 + 1>> compute_intersection_of_cells( @@ -663,7 +718,6 @@ namespace CGALWrappers const auto vertices0 = CGALWrappers::get_vertices_in_cgal_order( cell0, mapping0); - const auto vertices1 = CGALWrappers::get_vertices_in_cgal_order( cell1, mapping1);