From 41630d07d67adc1361e9d18d7b3368add10fda16 Mon Sep 17 00:00:00 2001 From: Johannes Heinz Date: Wed, 7 Sep 2022 14:16:59 +0200 Subject: [PATCH] add additional interface and fix subminor version --- cmake/modules/FindCGAL.cmake | 2 +- include/deal.II/cgal/intersections.h | 18 ++ include/deal.II/cgal/utilities.h | 33 ++++ source/cgal/CMakeLists.txt | 1 + source/cgal/intersections.cc | 270 +++++++++++---------------- source/cgal/intersections.inst.in | 31 +++ 6 files changed, 192 insertions(+), 163 deletions(-) create mode 100644 source/cgal/intersections.inst.in diff --git a/cmake/modules/FindCGAL.cmake b/cmake/modules/FindCGAL.cmake index 820f27a57f..34692bee59 100644 --- a/cmake/modules/FindCGAL.cmake +++ b/cmake/modules/FindCGAL.cmake @@ -71,7 +71,7 @@ IF(DEAL_II_HAVE_CXX17 AND NOT FEATURE_BOOST_BUNDLED_CONFIGURED) CGAL_VERSION_SUBMINOR "${CGAL_VERSION}" ) - IF("${CGAL_VERSION_SUBMINOR}" STREQUAL "") + IF("${CGAL_VERSION_SUBMINOR}" MATCHES "^(|${CGAL_VERSION})$") SET(CGAL_VERSION_SUBMINOR "0") ENDIF() ENDIF() diff --git a/include/deal.II/cgal/intersections.h b/include/deal.II/cgal/intersections.h index 5a7f3b4b17..8254d10e35 100644 --- a/include/deal.II/cgal/intersections.h +++ b/include/deal.II/cgal/intersections.h @@ -28,6 +28,24 @@ 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 * vector of simplices, each one identified by an array of deal.II Points. All diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 916a8e4437..4e8bc23a80 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -592,6 +592,39 @@ namespace CGALWrappers Assert(false, ExcInternalError()); } } + + + /** + * 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. + */ + template + std::array, n_vertices> + get_vertices_in_cgal_order( + const typename dealii::Triangulation::cell_iterator &cell, + const Mapping & mapping) + { + // Elements have to be rectangular or simplices + Assert(n_vertices == std::pow(2, dim) || n_vertices == dim + 1, + ExcNotImplemented()); + AssertDimension(mapping.get_vertices(cell).size(), n_vertices); + + std::array, n_vertices> vertices; + + std::copy_n(mapping.get_vertices(cell).begin(), + vertices.size(), + vertices.begin()); + + if (ReferenceCell::n_vertices_to_type(dim, n_vertices) == + ReferenceCells::Quadrilateral) + std::swap(vertices[2], vertices[3]); + + return vertices; + } + } // namespace CGALWrappers # endif diff --git a/source/cgal/CMakeLists.txt b/source/cgal/CMakeLists.txt index 2c951a6b36..3d16483339 100644 --- a/source/cgal/CMakeLists.txt +++ b/source/cgal/CMakeLists.txt @@ -26,6 +26,7 @@ SET(_src SET(_inst surface_mesh.inst.in +intersections.inst.in ) FILE(GLOB _header diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc index a9b23fe29b..7a859849b9 100644 --- a/source/cgal/intersections.cc +++ b/source/cgal/intersections.cc @@ -19,7 +19,6 @@ #include - #ifdef DEAL_II_WITH_CGAL # include @@ -56,8 +55,6 @@ # include # include - - DEAL_II_NAMESPACE_OPEN namespace CGALWrappers @@ -144,8 +141,6 @@ namespace CGALWrappers } } - - void mark_domains(CDT &cdt) { @@ -205,8 +200,6 @@ namespace CGALWrappers return CGAL::intersection(triangle1, triangle2); } - - boost::optional> compute_intersection(const std::array, 3> &first_simplex, const std::array, 2> &second_simplex) @@ -234,8 +227,6 @@ namespace CGALWrappers return CGAL::intersection(segm, triangle); } - - // rectangle-rectangle std::vector compute_intersection(const std::array, 4> &first_simplex, @@ -264,14 +255,11 @@ namespace CGALWrappers return poly_list; } - - boost::optional> compute_intersection(const std::array, 2> &first_simplex, const std::array, 4> &second_simplex) { # if DEAL_II_CGAL_VERSION_GTE(5, 1, 5) - std::array pts0; std::array pts1; std::transform( @@ -294,9 +282,7 @@ namespace CGALWrappers CGALTetra tetra{pts0[0], pts0[1], pts0[2], pts0[3]}; CGALSegment3 segm{pts1[0], pts1[1]}; return CGAL::intersection(segm, tetra); - # else - Assert( false, ExcMessage( @@ -307,8 +293,6 @@ namespace CGALWrappers # endif } - - // tetra, triangle boost::optional std::vector, 3>> - compute_intersection_of_cells<2, 2, 2>( - const typename Triangulation<2, 2>::cell_iterator &cell0, - const typename Triangulation<2, 2>::cell_iterator &cell1, - const Mapping<2, 2> & mapping0, - const Mapping<2, 2> & mapping1, - const double tol) + compute_intersection_of_cells<2, 2, 2, 4, 4>( + const std::array, 4> &vertices0, + const std::array, 4> &vertices1, + const double tol) { - Assert(ReferenceCell::n_vertices_to_type( - 2, mapping0.get_vertices(cell0).size()) == - ReferenceCells::Quadrilateral, - ExcNotImplemented()); - Assert(ReferenceCell::n_vertices_to_type( - 2, mapping1.get_vertices(cell1).size()) == - ReferenceCells::Quadrilateral, - ExcNotImplemented()); - - std::array, 4> vertices0, vertices1; - std::copy_n(mapping0.get_vertices(cell0).begin(), 4, vertices0.begin()); - std::copy_n(mapping1.get_vertices(cell1).begin(), 4, vertices1.begin()); - - std::swap(vertices0[2], vertices0[3]); - std::swap(vertices1[2], vertices1[3]); const auto intersection_test = internal::compute_intersection(vertices0, vertices1); @@ -435,33 +400,14 @@ namespace CGALWrappers } } - - // Specialization for quad \cap line template <> std::vector, 2>> - compute_intersection_of_cells<2, 1, 2>( - const typename Triangulation<2, 2>::cell_iterator &cell0, - const typename Triangulation<1, 2>::cell_iterator &cell1, - const Mapping<2, 2> & mapping0, - const Mapping<1, 2> & mapping1, - const double tol) + compute_intersection_of_cells<2, 1, 2, 4, 2>( + const std::array, 4> &vertices0, + const std::array, 2> &vertices1, + const double tol) { - Assert(ReferenceCell::n_vertices_to_type( - 2, mapping0.get_vertices(cell0).size()) == - ReferenceCells::Quadrilateral, - ExcNotImplemented()); - Assert(ReferenceCell::n_vertices_to_type( - 1, mapping1.get_vertices(cell1).size()) == ReferenceCells::Line, - ExcNotImplemented()); - - std::array, 4> vertices0; - std::array, 2> vertices1; - std::copy_n(mapping0.get_vertices(cell0).begin(), 4, vertices0.begin()); - std::copy_n(mapping1.get_vertices(cell1).begin(), 2, vertices1.begin()); - - std::swap(vertices0[2], vertices0[3]); - std::array pts; std::transform( vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<2> &p) { @@ -496,34 +442,15 @@ namespace CGALWrappers return vertices; } - - // specialization for hex \cap line template <> std::vector, 2>> - compute_intersection_of_cells<3, 1, 3>( - const typename Triangulation<3, 3>::cell_iterator &cell0, - const typename Triangulation<1, 3>::cell_iterator &cell1, - const Mapping<3, 3> & mapping0, - const Mapping<1, 3> & mapping1, - const double tol) + compute_intersection_of_cells<3, 1, 3, 8, 2>( + const std::array, 8> &vertices0, + const std::array, 2> &vertices1, + const double tol) { # if DEAL_II_CGAL_VERSION_GTE(5, 1, 5) - - Assert(ReferenceCell::n_vertices_to_type( - 3, mapping0.get_vertices(cell0).size()) == - ReferenceCells::Hexahedron, - ExcNotImplemented()); - Assert(ReferenceCell::n_vertices_to_type( - 1, mapping1.get_vertices(cell1).size()) == ReferenceCells::Line, - ExcNotImplemented()); - (void)tol; - - std::array, 8> vertices0; // 8 vertices of the hex - std::array, 2> vertices1; // 2 endpoints of the segment - std::copy_n(mapping0.get_vertices(cell0).begin(), 8, vertices0.begin()); - std::copy_n(mapping1.get_vertices(cell1).begin(), 2, vertices1.begin()); - std::array pts; std::transform( vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<3> &p) { @@ -566,43 +493,21 @@ namespace CGALWrappers false, ExcMessage( "This function requires a version of CGAL greater or equal than 5.1.5.")); - (void)cell0; - (void)cell1; - (void)mapping0; - (void)mapping1; + (void)vertices0; + (void)vertices1; (void)tol; return {}; - # endif } - - template <> std::vector, 3>> - compute_intersection_of_cells<3, 2, 3>( - const typename Triangulation<3, 3>::cell_iterator &cell0, - const typename Triangulation<2, 3>::cell_iterator &cell1, - const Mapping<3, 3> & mapping0, - const Mapping<2, 3> & mapping1, - const double tol) + compute_intersection_of_cells<3, 2, 3, 8, 4>( + const std::array, 8> &vertices0, + const std::array, 4> &vertices1, + const double tol) { # if DEAL_II_CGAL_VERSION_GTE(5, 1, 5) - - Assert(ReferenceCell::n_vertices_to_type( - 3, mapping0.get_vertices(cell0).size()) == - ReferenceCells::Hexahedron, - ExcNotImplemented()); - Assert(ReferenceCell::n_vertices_to_type( - 2, mapping1.get_vertices(cell1).size()) == - ReferenceCells::Quadrilateral, - ExcNotImplemented()); - - std::array, 8> vertices0; // 8 vertices of the hex - std::array, 4> vertices1; // 4 vertices of the quad - std::copy_n(mapping0.get_vertices(cell0).begin(), 8, vertices0.begin()); - std::copy_n(mapping1.get_vertices(cell1).begin(), 4, vertices1.begin()); - std::array pts_hex; std::array pts_quad; std::transform( @@ -621,7 +526,6 @@ namespace CGALWrappers return CGALWrappers::dealii_point_to_cgal_point(p); }); - // Subdivide hex into tetrahedrons std::vector, 3>> vertices; Triangulation3_exact tria; @@ -685,71 +589,96 @@ namespace CGALWrappers } return vertices; - # else - Assert( false, ExcMessage( "This function requires a version of CGAL greater or equal than 5.1.5.")); - (void)cell0; - (void)cell1; - (void)mapping0; - (void)mapping1; + (void)vertices0; + (void)vertices1; (void)tol; return {}; # endif } - - template <> std::vector, 4>> - compute_intersection_of_cells<3, 3, 3>( - const typename Triangulation<3, 3>::cell_iterator &cell0, - const typename Triangulation<3, 3>::cell_iterator &cell1, - const Mapping<3, 3> & mapping0, - const Mapping<3, 3> & mapping1, - const double tol) + compute_intersection_of_cells<3, 3, 3, 8, 8>( + const std::array, 8> &vertices0, + const std::array, 8> &vertices1, + const double tol) { - Assert(ReferenceCell::n_vertices_to_type( - 3, mapping0.get_vertices(cell0).size()) == - ReferenceCells::Hexahedron, + // 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 {}; + } + + template + std::vector, dim1 + 1>> + compute_intersection_of_cells( + const typename Triangulation::cell_iterator &cell0, + const typename Triangulation::cell_iterator &cell1, + const Mapping & mapping0, + const Mapping & mapping1, + const double tol) + { + Assert(mapping0.get_vertices(cell0).size() == std::pow(2, dim0), ExcNotImplemented()); - Assert(ReferenceCell::n_vertices_to_type( - 3, mapping1.get_vertices(cell1).size()) == - ReferenceCells::Hexahedron, + Assert(mapping1.get_vertices(cell1).size() == std::pow(2, dim1), ExcNotImplemented()); - 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; - } + + const auto vertices0 = + CGALWrappers::get_vertices_in_cgal_order( + cell0, mapping0); + + const auto vertices1 = + CGALWrappers::get_vertices_in_cgal_order( + cell1, mapping1); + + return compute_intersection_of_cells(vertices0, + vertices1, + tol); } +# include "intersections.inst" + } // namespace CGALWrappers DEAL_II_NAMESPACE_CLOSE @@ -758,6 +687,23 @@ DEAL_II_NAMESPACE_CLOSE DEAL_II_NAMESPACE_OPEN +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; + AssertThrow(false, ExcNeedsCGAL()); +} + template std::vector, dim1 + 1>> compute_intersection_of_cells( diff --git a/source/cgal/intersections.inst.in b/source/cgal/intersections.inst.in new file mode 100644 index 0000000000..d550252e0f --- /dev/null +++ b/source/cgal/intersections.inst.in @@ -0,0 +1,31 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS) + { +#if dim0 <= spacedim && dim0 >= dim1 + + template std::vector, dim1 + 1>> + compute_intersection_of_cells( + const typename Triangulation::cell_iterator &cell0, + const typename Triangulation::cell_iterator &cell1, + const Mapping & mapping0, + const Mapping & mapping1, + const double tol); + +#endif + } -- 2.39.5