From 6598e76e282292bb67c3d4c0060df27e707263cd Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Wed, 7 Sep 2022 11:41:20 +0200 Subject: [PATCH] code review --- cmake/modules/FindCGAL.cmake | 74 +++++--- doc/doxygen/options.dox.in | 1 + include/deal.II/base/config.h.in | 15 ++ include/deal.II/base/exceptions.h | 9 + include/deal.II/cgal/intersections.h | 7 +- source/cgal/intersections.cc | 159 +++++++++++------- tests/cgal/cgal_intersect_simplices_1d_2d.cc | 31 +++- .../cgal_intersect_simplices_1d_2d.output | 3 +- tests/cgal/cgal_intersect_simplices_1d_3d.cc | 30 ++++ .../cgal_intersect_simplices_1d_3d.output | 1 + tests/cgal/cgal_intersect_simplices_2d_2d.cc | 35 +++- .../cgal_intersect_simplices_2d_2d.output | 1 + tests/cgal/cgal_intersect_simplices_3d_3d.cc | 32 +++- .../cgal_intersect_simplices_3d_3d.output | 1 + 14 files changed, 304 insertions(+), 95 deletions(-) diff --git a/cmake/modules/FindCGAL.cmake b/cmake/modules/FindCGAL.cmake index f3aa3e9631..820f27a57f 100644 --- a/cmake/modules/FindCGAL.cmake +++ b/cmake/modules/FindCGAL.cmake @@ -1,24 +1,27 @@ -## --------------------------------------------------------------------- -## -## 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 # -# This module exports +# This module exports: # -# CGAL_INCLUDE_DIRS +# CGAL_INCLUDE_DIRS +# CGAL_VERSION_MAJOR +# CGAL_VERSION_MINOR +# CGAL_VERSION_SUBMINOR # SET(CGAL_DIR "" CACHE PATH "An optional hint to a CGAL installation") @@ -38,6 +41,7 @@ IF(DEAL_II_HAVE_CXX17 AND NOT FEATURE_BOOST_BUNDLED_CONFIGURED) LIST(REMOVE_ITEM CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake/modules/) SET(CGAL_DO_NOT_WARN_ABOUT_CMAKE_BUILD_TYPE ON) FIND_PACKAGE(CGAL QUIET) + # Check version manually. Older binary distros don't do this properly. IF(CGAL_MAJOR_VERSION LESS 5) SET(CGAL_FOUND FALSE) @@ -45,13 +49,33 @@ IF(DEAL_II_HAVE_CXX17 AND NOT FEATURE_BOOST_BUNDLED_CONFIGURED) SET(CGAL_LIBRARIES "-NOTFOUND") MESSAGE(STATUS "CGAL wrappers require CGAL version 5 and above.") ENDIF() + LIST(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake/modules/) + IF(CGAL_FOUND) GET_TARGET_PROPERTY(CGAL_LIBRARIES CGAL::CGAL INTERFACE_LINK_LIBRARIES) - IF(CGAL_VERSION GREATER_EQUAL 5.1.5) - ADD_COMPILE_DEFINITIONS(CGAL_GEQ_515) + IF(DEFINED CGAL_VERSION) + SET(CGAL_VERSION "${CGAL_VERSION}") + + STRING(REGEX REPLACE + "^([0-9]+).*$" "\\1" + CGAL_VERSION_MAJOR "${CGAL_VERSION}") + + STRING(REGEX REPLACE + "^[0-9]+\\.([0-9]+).*$" "\\1" + CGAL_VERSION_MINOR "${CGAL_VERSION}") + + STRING(REGEX REPLACE + "^[0-9]+\\.[0-9]+\\.([0-9]+).*$" "\\1" + CGAL_VERSION_SUBMINOR "${CGAL_VERSION}" + ) + + IF("${CGAL_VERSION_SUBMINOR}" STREQUAL "") + SET(CGAL_VERSION_SUBMINOR "0") + ENDIF() ENDIF() + # Make sure we dont' pass Boost::Boost over to deal.II. LIST(FILTER CGAL_LIBRARIES EXCLUDE REGEX "::") ENDIF() @@ -64,12 +88,12 @@ ENDIF() DEAL_II_PACKAGE_HANDLE(CGAL INCLUDE_DIRS - REQUIRED CGAL_INCLUDE_DIRS + REQUIRED CGAL_INCLUDE_DIRS LIBRARIES - OPTIONAL CGAL_LIBRARIES + OPTIONAL CGAL_LIBRARIES USER_INCLUDE_DIRS - REQUIRED CGAL_INCLUDE_DIRS - CLEAR - CGAL_INCLUDE_DIRS - CGAL_LIBRARIES - ) + REQUIRED CGAL_INCLUDE_DIRS + CLEAR + CGAL_INCLUDE_DIRS + CGAL_LIBRARIES +) diff --git a/doc/doxygen/options.dox.in b/doc/doxygen/options.dox.in index ac58e22a2e..9e0271628e 100644 --- a/doc/doxygen/options.dox.in +++ b/doc/doxygen/options.dox.in @@ -195,6 +195,7 @@ PREDEFINED = DOXYGEN=1 \ DEAL_II_WITH_ASSIMP=1 \ DEAL_II_WITH_BOOST=1 \ DEAL_II_WITH_CGAL=1 \ + DEAL_II_CGAL_VERSION_GTE=1 \ DEAL_II_WITH_TASKFLOW=1 \ DEAL_II_WITH_COMPLEX_VALUES=1 \ DEAL_II_WITH_CUDA=1 \ diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index 30d2e672ff..a201f301de 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -415,6 +415,21 @@ (major)*10000 + (minor)*100 + (subminor)) #endif +/* + * CGAL: + */ + +#ifdef DEAL_II_WITH_CGAL +# define DEAL_II_CGAL_VERSION_MAJOR @CGAL_VERSION_MAJOR@ +# define DEAL_II_CGAL_VERSION_MINOR @CGAL_VERSION_MINOR@ +# define DEAL_II_CGAL_VERSION_SUBMINOR @CGAL_VERSION_SUBMINOR@ + +# define DEAL_II_CGAL_VERSION_GTE(major, minor, subminor) \ + ((DEAL_II_CGAL_VERSION_MAJOR * 10000 + DEAL_II_CGAL_VERSION_MINOR * 100 + \ + DEAL_II_CGAL_VERSION_SUBMINOR) >= \ + (major)*10000 + (minor)*100 + (subminor)) +#endif + /* * MPI */ diff --git a/include/deal.II/base/exceptions.h b/include/deal.II/base/exceptions.h index e4e62af343..f33ba5fbc3 100644 --- a/include/deal.II/base/exceptions.h +++ b/include/deal.II/base/exceptions.h @@ -1160,6 +1160,15 @@ namespace StandardExceptions "was configured to use Trilinos' SEACAS library (which provides ExodusII), " "but cmake did not find a valid SEACAS library."); + /** + * This function requires support for the CGAL library. + */ + DeclExceptionMsg( + ExcNeedsCGAL, + "You are attempting to use functionality that is only available " + "if deal.II was configured to use CGAL, but cmake did not " + "find a valid CGAL library."); + #ifdef DEAL_II_WITH_MPI /** * Exception for MPI errors. This exception is only defined if diff --git a/include/deal.II/cgal/intersections.h b/include/deal.II/cgal/intersections.h index 96883afda5..5a7f3b4b17 100644 --- a/include/deal.II/cgal/intersections.h +++ b/include/deal.II/cgal/intersections.h @@ -32,9 +32,10 @@ namespace CGALWrappers * 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 * the simplices together are a subdivision of the intersection. If cells are - * non-affine, a geometrical error will be necessarily introduced. If the - * measure of one of the simplices is below a certain treshold, than it is - * discarded. + * non-affine, a geometrical error is introduced. If the + * measure of one of the simplices is below a certain treshold which defaults + * to 1e-9, then it is discarded. In case the two cells are disjoint, an empty + * array is returned. * * @note If the two cells have different intrinsic dimensions, than the * first iterator is assumed to have with higher intrinsic dimension. diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc index 9537cae745..a9b23fe29b 100644 --- a/source/cgal/intersections.cc +++ b/source/cgal/intersections.cc @@ -56,57 +56,57 @@ # include # include -using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt; -using K_exact = CGAL::Exact_predicates_exact_constructions_kernel; -using K_inexact = CGAL::Exact_predicates_inexact_constructions_kernel; -using CGALPolygon = CGAL::Polygon_2; -using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; -using CGALTriangle2 = K::Triangle_2; -using CGALTriangle3 = K::Triangle_3; -using CGALTriangle3_exact = K_exact::Triangle_3; -using CGALPoint2 = K::Point_2; -using CGALPoint3 = K::Point_3; -using CGALPoint3_exact = K_exact::Point_3; -using CGALSegment2 = K::Segment_2; -using Surface_mesh = CGAL::Surface_mesh; -using CGALSegment3 = K::Segment_3; -using CGALSegment3_exact = K_exact::Segment_3; -using CGALTetra = K::Tetrahedron_3; -using CGALTetra_exact = K_exact::Tetrahedron_3; -using Triangulation2 = CGAL::Triangulation_2; -using Triangulation3 = CGAL::Triangulation_3; -using Triangulation3_exact = CGAL::Triangulation_3; -using Triangulation3_inexact = CGAL::Triangulation_3; - -struct FaceInfo2 -{ - FaceInfo2() - {} - int nesting_level; - bool - in_domain() - { - return nesting_level % 2 == 1; - } -}; - -using Vb = CGAL::Triangulation_vertex_base_2; -using Fbb = CGAL::Triangulation_face_base_with_info_2; -using CFb = CGAL::Constrained_triangulation_face_base_2; -using Fb = CGAL::Delaunay_mesh_face_base_2; -using Tds = CGAL::Triangulation_data_structure_2; -using Itag = CGAL::Exact_predicates_tag; -using CDT = CGAL::Constrained_Delaunay_triangulation_2; -using Criteria = CGAL::Delaunay_mesh_size_criteria_2; -using Vertex_handle = CDT::Vertex_handle; -using Face_handle = CDT::Face_handle; - DEAL_II_NAMESPACE_OPEN namespace CGALWrappers { + using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt; + using K_exact = CGAL::Exact_predicates_exact_constructions_kernel; + using K_inexact = CGAL::Exact_predicates_inexact_constructions_kernel; + using CGALPolygon = CGAL::Polygon_2; + using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; + using CGALTriangle2 = K::Triangle_2; + using CGALTriangle3 = K::Triangle_3; + using CGALTriangle3_exact = K_exact::Triangle_3; + using CGALPoint2 = K::Point_2; + using CGALPoint3 = K::Point_3; + using CGALPoint3_exact = K_exact::Point_3; + using CGALSegment2 = K::Segment_2; + using Surface_mesh = CGAL::Surface_mesh; + using CGALSegment3 = K::Segment_3; + using CGALSegment3_exact = K_exact::Segment_3; + using CGALTetra = K::Tetrahedron_3; + using CGALTetra_exact = K_exact::Tetrahedron_3; + using Triangulation2 = CGAL::Triangulation_2; + using Triangulation3 = CGAL::Triangulation_3; + using Triangulation3_exact = CGAL::Triangulation_3; + using Triangulation3_inexact = CGAL::Triangulation_3; + + struct FaceInfo2 + { + FaceInfo2() + {} + int nesting_level; + bool + in_domain() + { + return nesting_level % 2 == 1; + } + }; + + using Vb = CGAL::Triangulation_vertex_base_2; + using Fbb = CGAL::Triangulation_face_base_with_info_2; + using CFb = CGAL::Constrained_triangulation_face_base_2; + using Fb = CGAL::Delaunay_mesh_face_base_2; + using Tds = CGAL::Triangulation_data_structure_2; + using Itag = CGAL::Exact_predicates_tag; + using CDT = CGAL::Constrained_Delaunay_triangulation_2; + using Criteria = CGAL::Delaunay_mesh_size_criteria_2; + using Vertex_handle = CDT::Vertex_handle; + using Face_handle = CDT::Face_handle; + namespace internal { void @@ -237,7 +237,7 @@ namespace CGALWrappers // rectangle-rectangle - decltype(auto) + std::vector compute_intersection(const std::array, 4> &first_simplex, const std::array, 4> &second_simplex) { @@ -270,7 +270,7 @@ namespace CGALWrappers compute_intersection(const std::array, 2> &first_simplex, const std::array, 4> &second_simplex) { -# if defined(CGAL_GEQ_515) +# if DEAL_II_CGAL_VERSION_GTE(5, 1, 5) std::array pts0; std::array pts1; @@ -317,7 +317,7 @@ namespace CGALWrappers compute_intersection(const std::array, 3> &first_simplex, const std::array, 4> &second_simplex) { -# if defined(CGAL_GEQ_515) +# if DEAL_II_CGAL_VERSION_GTE(5, 1, 5) std::array pts0; std::array pts1; std::transform( @@ -363,6 +363,15 @@ namespace CGALWrappers const Mapping<2, 2> & mapping1, 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()); @@ -398,7 +407,7 @@ namespace CGALWrappers internal::mark_domains(cdt); std::array, 3> vertices; - for (Face_handle f : cdt.finite_face_handles()) + for (const Face_handle &f : cdt.finite_face_handles()) { if (f->info().in_domain() && CGAL::to_double(cdt.triangle(f).area()) > tol) @@ -422,7 +431,6 @@ namespace CGALWrappers } else { - Assert(false, ExcMessage("Cells do not intersect.")); return {}; } } @@ -439,6 +447,14 @@ namespace CGALWrappers const Mapping<1, 2> & mapping1, 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()); @@ -492,8 +508,17 @@ namespace CGALWrappers const Mapping<1, 3> & mapping1, const double tol) { -# if defined(CGAL_GEQ_515) - (void)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()); @@ -562,7 +587,16 @@ namespace CGALWrappers const Mapping<2, 3> & mapping1, const double tol) { -# if defined(CGAL_GEQ_515) +# 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 @@ -668,6 +702,7 @@ namespace CGALWrappers } + template <> std::vector, 4>> compute_intersection_of_cells<3, 3, 3>( @@ -677,6 +712,14 @@ namespace CGALWrappers const Mapping<3, 3> & mapping1, const double tol) { + Assert(ReferenceCell::n_vertices_to_type( + 3, mapping0.get_vertices(cell0).size()) == + ReferenceCells::Hexahedron, + ExcNotImplemented()); + Assert(ReferenceCell::n_vertices_to_type( + 3, mapping1.get_vertices(cell1).size()) == + ReferenceCells::Hexahedron, + 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); @@ -684,11 +727,11 @@ namespace CGALWrappers 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 - std::vector, 4>> vertices; - Triangulation3_inexact tria; + Triangulation3_inexact tria; tria.insert(sm.points().begin(), sm.points().end()); for (const auto &c : tria.finite_cell_handles()) { @@ -703,7 +746,7 @@ namespace CGALWrappers } else { - return {}; + return vertices; } } @@ -722,14 +765,14 @@ compute_intersection_of_cells( const typename Triangulation::cell_iterator &cell1, const Mapping & mapping0, const Mapping & mapping1, - const double tol = 1e-9) + const double tol) { (void)cell0; (void)cell1; (void)mapping0; (void)mapping1; (void)tol; - AssertThrow(false, ExcMessage("This function needs CGAL to be installed.")); + AssertThrow(false, ExcNeedsCGAL()); } DEAL_II_NAMESPACE_CLOSE diff --git a/tests/cgal/cgal_intersect_simplices_1d_2d.cc b/tests/cgal/cgal_intersect_simplices_1d_2d.cc index 5d01660fdf..b0d00ba33c 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_2d.cc +++ b/tests/cgal/cgal_intersect_simplices_1d_2d.cc @@ -82,18 +82,45 @@ test_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1) +void +test_failing_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, 1.5, 2.); + const double expected_area = 0.; + const auto cell0 = tria0.begin_active(); + const auto cell1 = tria1.begin_active(); + + const auto vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<2, 1, 2>(cell0, + cell1, + MappingQ1<2>(), + MappingQ1<1, 2>()); + + const auto quad = qgauss.mapped_quadrature(vec_of_arrays); + const double sum = + std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.); + assert(std::abs(expected_area - sum) < 1e-15); + deallog << "OK" << std::endl; +} + + + int main() { initlog(); Triangulation<2> tria0; // ambient space Triangulation<1, 2> tria1; // immersed grid - deallog << "Inside: " << std::endl; test_inside_intersection(tria0, tria1); tria0.clear(); tria1.clear(); - deallog << "Intersecting: " << std::endl; test_intersection(tria0, tria1); + + tria0.clear(); + tria1.clear(); + + test_failing_intersection(tria0, tria1); } diff --git a/tests/cgal/cgal_intersect_simplices_1d_2d.output b/tests/cgal/cgal_intersect_simplices_1d_2d.output index 6dd05e6537..fb71de2867 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_2d.output +++ b/tests/cgal/cgal_intersect_simplices_1d_2d.output @@ -1,5 +1,4 @@ -DEAL::Inside: DEAL::OK -DEAL::Intersecting: +DEAL::OK DEAL::OK diff --git a/tests/cgal/cgal_intersect_simplices_1d_3d.cc b/tests/cgal/cgal_intersect_simplices_1d_3d.cc index 69c8505f7e..d7bd153ac4 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_1d_3d.cc @@ -80,6 +80,31 @@ test_intersection(Triangulation<3> &tria0, Triangulation<1, 3> &tria1) +void +test_failing_intersection(Triangulation<3> &tria0, Triangulation<1, 3> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, 2., 3.); + + const auto cell0 = tria0.begin_active(); + const auto cell1 = tria1.begin_active(); + const double expected_measure = 0.; + + const auto vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<3, 1, 3>(cell0, + cell1, + MappingQ1<3>(), + MappingQ1<1, 3>()); + + const auto quad = QGaussSimplex<1>(1).mapped_quadrature(vec_of_arrays); + const double sum = + std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.); + assert(std::abs(sum - expected_measure) < 1e-15); + deallog << "OK" << std::endl; +} + + + int main() { @@ -88,7 +113,12 @@ main() Triangulation<1, 3> tria1; test_intersection_inside(tria0, tria1); + tria0.clear(); tria1.clear(); test_intersection(tria0, tria1); + + tria0.clear(); + tria1.clear(); + test_failing_intersection(tria0, tria1); } diff --git a/tests/cgal/cgal_intersect_simplices_1d_3d.output b/tests/cgal/cgal_intersect_simplices_1d_3d.output index 8b3b075900..fb71de2867 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_3d.output +++ b/tests/cgal/cgal_intersect_simplices_1d_3d.output @@ -1,3 +1,4 @@ DEAL::OK DEAL::OK +DEAL::OK diff --git a/tests/cgal/cgal_intersect_simplices_2d_2d.cc b/tests/cgal/cgal_intersect_simplices_2d_2d.cc index 0e9bf03bba..c4e99b60b6 100644 --- a/tests/cgal/cgal_intersect_simplices_2d_2d.cc +++ b/tests/cgal/cgal_intersect_simplices_2d_2d.cc @@ -30,7 +30,6 @@ #include "../tests.h" using namespace CGALWrappers; -static QGaussSimplex<2> qgauss(1); // use a degree 1 QGaussSimplex<2> void test_inside_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) @@ -49,7 +48,7 @@ test_inside_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) MappingQ1<2>()); - const auto quad = qgauss.mapped_quadrature(vec_of_arrays); + const auto quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays); const double sum = std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.); assert(std::abs(expected_area - sum) < 1e-15); @@ -74,7 +73,32 @@ test_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) MappingQ1<2>(), MappingQ1<2>()); - const auto quad = qgauss.mapped_quadrature(vec_of_arrays); + const auto quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays); + const double sum = + std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.); + assert(std::abs(expected_area - sum) < 1e-15); + deallog << "OK" << std::endl; +} + + + +void +test_failing_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, 1.5, 2.5); + const double expected_area = 0.; + + const auto cell0 = tria0.begin_active(); + const auto cell1 = tria1.begin_active(); + + const auto vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<2, 2, 2>(cell0, + cell1, + MappingQ1<2>(), + MappingQ1<2>()); + + const auto quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays); const double sum = std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.); assert(std::abs(expected_area - sum) < 1e-15); @@ -96,4 +120,9 @@ main() tria1.clear(); test_intersection(tria0, tria1); + + tria0.clear(); + tria1.clear(); + + test_failing_intersection(tria0, tria1); } diff --git a/tests/cgal/cgal_intersect_simplices_2d_2d.output b/tests/cgal/cgal_intersect_simplices_2d_2d.output index 8b3b075900..fb71de2867 100644 --- a/tests/cgal/cgal_intersect_simplices_2d_2d.output +++ b/tests/cgal/cgal_intersect_simplices_2d_2d.output @@ -1,3 +1,4 @@ DEAL::OK DEAL::OK +DEAL::OK diff --git a/tests/cgal/cgal_intersect_simplices_3d_3d.cc b/tests/cgal/cgal_intersect_simplices_3d_3d.cc index 639d30a128..d19f3f3ef9 100644 --- a/tests/cgal/cgal_intersect_simplices_3d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_3d_3d.cc @@ -27,7 +27,6 @@ #include "../tests.h" -static QGaussSimplex<3> qgauss(1); void test_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1) { @@ -43,7 +42,31 @@ test_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1) MappingQ1<3>(), MappingQ1<3>()); - const auto quad = qgauss.mapped_quadrature(vec_of_arrays); + const auto quad = QGaussSimplex<3>(1).mapped_quadrature(vec_of_arrays); + const double sum = + std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.); + assert(std::abs(sum - expected_measure) < 1e-15); + deallog << "OK" << std::endl; +} + + + +void +test_failing_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, 2.5, 3.); + static constexpr double expected_measure = 0.; + const auto cell0 = tria0.begin_active(); + const auto cell1 = tria1.begin_active(); + + const auto vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<3, 3, 3>(cell0, + cell1, + MappingQ1<3>(), + MappingQ1<3>()); + + const auto quad = QGaussSimplex<3>(1).mapped_quadrature(vec_of_arrays); const double sum = std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.); assert(std::abs(sum - expected_measure) < 1e-15); @@ -60,4 +83,9 @@ main() Triangulation<3> tria1; test_intersection(tria0, tria1); + + tria0.clear(); + tria1.clear(); + + test_failing_intersection(tria0, tria1); } diff --git a/tests/cgal/cgal_intersect_simplices_3d_3d.output b/tests/cgal/cgal_intersect_simplices_3d_3d.output index 0fd8fc12f0..8b3b075900 100644 --- a/tests/cgal/cgal_intersect_simplices_3d_3d.output +++ b/tests/cgal/cgal_intersect_simplices_3d_3d.output @@ -1,2 +1,3 @@ DEAL::OK +DEAL::OK -- 2.39.5