]> https://gitweb.dealii.org/ - dealii.git/commitdiff
code review
authorMarco Feder <marco.feder@sissa.it>
Wed, 7 Sep 2022 09:41:20 +0000 (11:41 +0200)
committerMarco Feder <marco.feder@sissa.it>
Wed, 7 Sep 2022 11:29:54 +0000 (13:29 +0200)
14 files changed:
cmake/modules/FindCGAL.cmake
doc/doxygen/options.dox.in
include/deal.II/base/config.h.in
include/deal.II/base/exceptions.h
include/deal.II/cgal/intersections.h
source/cgal/intersections.cc
tests/cgal/cgal_intersect_simplices_1d_2d.cc
tests/cgal/cgal_intersect_simplices_1d_2d.output
tests/cgal/cgal_intersect_simplices_1d_3d.cc
tests/cgal/cgal_intersect_simplices_1d_3d.output
tests/cgal/cgal_intersect_simplices_2d_2d.cc
tests/cgal/cgal_intersect_simplices_2d_2d.output
tests/cgal/cgal_intersect_simplices_3d_3d.cc
tests/cgal/cgal_intersect_simplices_3d_3d.output

index f3aa3e963197499bb57ac6fecc9a2cb57b7dbc89..820f27a57f2d4da8342f75afe2752ab94fd076ec 100644 (file)
@@ -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
+)
index ac58e22a2ed8e45116ea37d7d7f20d37e83f4251..9e0271628ef12cfe7345cbebf41c040c5c10d8ef 100644 (file)
@@ -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 \
index 30d2e672ff179f9e15935656706ea978e8c48c48..a201f301ded359061aca74ba3f76c8f4a12cc971 100644 (file)
     (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
  */
index e4e62af343fee2b73e7e8f30488ec11c557b2a56..f33ba5fbc369ae0553a8236cacca0892d406cebf 100644 (file)
@@ -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
index 96883afda5a96a3d14784ec8b133364e35bcc6ec..5a7f3b4b17e3fe282297c6762ab3f7f0fb0b278b 100644 (file)
@@ -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.
index 9537cae7453d5a3946b7c8e3396f753160932a63..a9b23fe29ba94334b8682e93b6d00af5c77ceb0f 100644 (file)
 #  include <fstream>
 #  include <type_traits>
 
-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<K>;
-using Polygon_with_holes_2   = CGAL::Polygon_with_holes_2<K>;
-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<K_inexact::Point_3>;
-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<K>;
-using Triangulation3         = CGAL::Triangulation_3<K>;
-using Triangulation3_exact   = CGAL::Triangulation_3<K_exact>;
-using Triangulation3_inexact = CGAL::Triangulation_3<K_inexact>;
-
-struct FaceInfo2
-{
-  FaceInfo2()
-  {}
-  int nesting_level;
-  bool
-  in_domain()
-  {
-    return nesting_level % 2 == 1;
-  }
-};
-
-using Vb            = CGAL::Triangulation_vertex_base_2<K>;
-using Fbb           = CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>;
-using CFb           = CGAL::Constrained_triangulation_face_base_2<K, Fbb>;
-using Fb            = CGAL::Delaunay_mesh_face_base_2<K, CFb>;
-using Tds           = CGAL::Triangulation_data_structure_2<Vb, Fb>;
-using Itag          = CGAL::Exact_predicates_tag;
-using CDT           = CGAL::Constrained_Delaunay_triangulation_2<K, Tds, Itag>;
-using Criteria      = CGAL::Delaunay_mesh_size_criteria_2<CDT>;
-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<K>;
+  using Polygon_with_holes_2   = CGAL::Polygon_with_holes_2<K>;
+  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<K_inexact::Point_3>;
+  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<K>;
+  using Triangulation3         = CGAL::Triangulation_3<K>;
+  using Triangulation3_exact   = CGAL::Triangulation_3<K_exact>;
+  using Triangulation3_inexact = CGAL::Triangulation_3<K_inexact>;
+
+  struct FaceInfo2
+  {
+    FaceInfo2()
+    {}
+    int nesting_level;
+    bool
+    in_domain()
+    {
+      return nesting_level % 2 == 1;
+    }
+  };
+
+  using Vb       = CGAL::Triangulation_vertex_base_2<K>;
+  using Fbb      = CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>;
+  using CFb      = CGAL::Constrained_triangulation_face_base_2<K, Fbb>;
+  using Fb       = CGAL::Delaunay_mesh_face_base_2<K, CFb>;
+  using Tds      = CGAL::Triangulation_data_structure_2<Vb, Fb>;
+  using Itag     = CGAL::Exact_predicates_tag;
+  using CDT      = CGAL::Constrained_Delaunay_triangulation_2<K, Tds, Itag>;
+  using Criteria = CGAL::Delaunay_mesh_size_criteria_2<CDT>;
+  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<Polygon_with_holes_2>
     compute_intersection(const std::array<Point<2>, 4> &first_simplex,
                          const std::array<Point<2>, 4> &second_simplex)
     {
@@ -270,7 +270,7 @@ namespace CGALWrappers
     compute_intersection(const std::array<Point<3>, 2> &first_simplex,
                          const std::array<Point<3>, 4> &second_simplex)
     {
-#  if defined(CGAL_GEQ_515)
+#  if DEAL_II_CGAL_VERSION_GTE(5, 1, 5)
 
       std::array<CGALPoint3, 4> pts0;
       std::array<CGALPoint3, 2> pts1;
@@ -317,7 +317,7 @@ namespace CGALWrappers
     compute_intersection(const std::array<Point<3>, 3> &first_simplex,
                          const std::array<Point<3>, 4> &second_simplex)
     {
-#  if defined(CGAL_GEQ_515)
+#  if DEAL_II_CGAL_VERSION_GTE(5, 1, 5)
       std::array<CGALPoint3, 4> pts0;
       std::array<CGALPoint3, 3> 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<Point<2>, 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<Point<2>, 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<Point<2>, 4> vertices0;
     std::array<Point<2>, 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<Point<3>, 8> vertices0; // 8 vertices of the hex
     std::array<Point<3>, 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<Point<3>, 8> vertices0; // 8 vertices of the hex
     std::array<Point<3>, 4> vertices1; // 4 vertices of the quad
@@ -668,6 +702,7 @@ namespace CGALWrappers
   }
 
 
+
   template <>
   std::vector<std::array<Point<3>, 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<std::array<Point<3>, 4>> vertices;
     if (CGAL::Polygon_mesh_processing::volume(sm) > tol)
       {
         // Collect tetrahedrons
-        std::vector<std::array<Point<3>, 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<dim1, spacedim>::cell_iterator &cell1,
   const Mapping<dim0, spacedim> &                              mapping0,
   const Mapping<dim1, spacedim> &                              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
index 5d01660fdf42d270f036c60a19496415438162d3..b0d00ba33cb8cf32331b19bda60c03059ad9ac08 100644 (file)
@@ -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);
 }
index 6dd05e653700ccadb3ddb0ddb4217c88b0ffa5f5..fb71de2867cdd56077fa628bd1498d2d19382019 100644 (file)
@@ -1,5 +1,4 @@
 
-DEAL::Inside: 
 DEAL::OK
-DEAL::Intersecting: 
+DEAL::OK
 DEAL::OK
index 69c8505f7e262077d601154a9e5f51e0911cfdc2..d7bd153ac43bb9cf94ba1dd87d3ec7e0a60d6332 100644 (file)
@@ -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);
 }
index 8b3b075900ce960fcb04b78ce0e95f33d6db4ed2..fb71de2867cdd56077fa628bd1498d2d19382019 100644 (file)
@@ -1,3 +1,4 @@
 
 DEAL::OK
 DEAL::OK
+DEAL::OK
index 0e9bf03bba5fd9d97bc6d063c78903ef256b3520..c4e99b60b6845849aa65be58ec66b88cf7038065 100644 (file)
@@ -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);
 }
index 8b3b075900ce960fcb04b78ce0e95f33d6db4ed2..fb71de2867cdd56077fa628bd1498d2d19382019 100644 (file)
@@ -1,3 +1,4 @@
 
 DEAL::OK
 DEAL::OK
+DEAL::OK
index 639d30a128042b3cfdd9bb71282870d6db8a4e44..d19f3f3ef9c00f701b5e9504170ddfb397b229aa 100644 (file)
@@ -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);
 }
index 0fd8fc12f0b442283edd8868867114c242b04d11..8b3b075900ce960fcb04b78ce0e95f33d6db4ed2 100644 (file)
@@ -1,2 +1,3 @@
 
 DEAL::OK
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.