fix 2d/3d if not in xy plane

diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc
index d9b618c7842d43abd044d561f9d5d7a5a7e50822..9537cae7453d5a3946b7c8e3396f753160932a63 100644 (file)
--- a/source/cgal/intersections.cc
+++ b/source/cgal/intersections.cc
@@ -77,8 +77,6 @@ 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>;
-typedef CGAL::Projection_traits_xy_3<K_exact> Gt;
-typedef CGAL::Delaunay_triangulation_2<Gt>    Delaunay;
 
 struct FaceInfo2
 {
@@ -596,13 +594,14 @@ namespace CGALWrappers
     tria.insert(pts_hex.begin(), pts_hex.end());
 
     // Subdivide quad into triangles
-    Delaunay tria_quad(pts_quad.begin(), pts_quad.end());
+    Triangulation3_exact tria_quad;
+    tria_quad.insert(pts_quad.begin(), pts_quad.end());
 
     for (const auto &c : tria.finite_cell_handles())
       {
         const auto &tet = tria.tetrahedron(c);
 
-        for (const auto f : tria_quad.finite_face_handles())
+        for (const auto f : tria_quad.finite_facets())
           {
             if (CGAL::do_intersect(tet, tria_quad.triangle(f)))
               {
@@ -624,22 +623,24 @@ namespace CGALWrappers
                 if (const std::vector<CGALPoint3_exact> *vps =
                       boost::get<std::vector<CGALPoint3_exact>>(&*intersection))
                   {
-                    Delaunay tria_inter(vps->begin(), vps->end());
-                    for (auto it = tria_inter.finite_faces_begin();
-                         it != tria_inter.finite_faces_end();
+                    Triangulation3_exact tria_inter;
+                    tria_inter.insert(vps->begin(), vps->end());
+
+                    for (auto it = tria_inter.finite_facets_begin();
+                         it != tria_inter.finite_facets_end();
                          ++it)
                       {
-                        if (CGAL::to_double(
-                              tria_inter.triangle(it).squared_area()) >
+                        const auto triangle = tria_inter.triangle(*it);
+                        if (CGAL::to_double(triangle.squared_area()) >
                             tol * tol)
                           {
                             std::array<Point<3>, 3> verts = {
                               {CGALWrappers::cgal_point_to_dealii_point<3>(
-                                 it->vertex(0)->point()),
+                                 triangle[0]),
                                CGALWrappers::cgal_point_to_dealii_point<3>(
-                                 it->vertex(1)->point()),
+                                 triangle[1]),
                                CGALWrappers::cgal_point_to_dealii_point<3>(
-                                 it->vertex(2)->point())}};
+                                 triangle[2])}};
 
                             vertices.push_back(verts);
                           }

diff --git a/tests/cgal/cgal_intersect_simplices_1d_2d.cc b/tests/cgal/cgal_intersect_simplices_1d_2d.cc
index 435bbcab0807f700e86fd7d3d8702fdddd250413..5d01660fdf42d270f036c60a19496415438162d3 100644 (file)
--- a/tests/cgal/cgal_intersect_simplices_1d_2d.cc
+++ b/tests/cgal/cgal_intersect_simplices_1d_2d.cc
@@ -19,6 +19,8 @@
 
 #include <deal.II/base/quadrature_lib.h>
 
+#include <deal.II/fe/mapping_q1.h>
+
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>

diff --git a/tests/cgal/cgal_intersect_simplices_1d_3d.cc b/tests/cgal/cgal_intersect_simplices_1d_3d.cc
index bfd00f9d7a91c24fc889282619988c087787c421..69c8505f7e262077d601154a9e5f51e0911cfdc2 100644 (file)
--- a/tests/cgal/cgal_intersect_simplices_1d_3d.cc
+++ b/tests/cgal/cgal_intersect_simplices_1d_3d.cc
@@ -18,6 +18,8 @@
 
 #include <deal.II/base/quadrature_lib.h>
 
+#include <deal.II/fe/mapping_q1.h>
+
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_out.h>
 #include <deal.II/grid/grid_tools.h>

diff --git a/tests/cgal/cgal_intersect_simplices_2d_2d.cc b/tests/cgal/cgal_intersect_simplices_2d_2d.cc
index d2d0cbf36bebba820e9016640482923f1d50d8e4..0e9bf03bba5fd9d97bc6d063c78903ef256b3520 100644 (file)
--- a/tests/cgal/cgal_intersect_simplices_2d_2d.cc
+++ b/tests/cgal/cgal_intersect_simplices_2d_2d.cc
@@ -19,6 +19,8 @@
 
 #include <deal.II/base/quadrature_lib.h>
 
+#include <deal.II/fe/mapping_q1.h>
+
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>

diff --git a/tests/cgal/cgal_intersect_simplices_2d_3d.cc b/tests/cgal/cgal_intersect_simplices_2d_3d.cc
index 38070da4a641d4ff1c95fa81fb17d77ebc73b1e7..be43f35fd7ff47ce7e8c2f32abe6e691beae168a 100644 (file)
--- a/tests/cgal/cgal_intersect_simplices_2d_3d.cc
+++ b/tests/cgal/cgal_intersect_simplices_2d_3d.cc
@@ -18,6 +18,8 @@
 
 #include <deal.II/base/quadrature_lib.h>
 
+#include <deal.II/fe/mapping_q1.h>
+
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_out.h>
 #include <deal.II/grid/grid_tools.h>
@@ -27,22 +29,75 @@
 
 #include "../tests.h"
 
+struct swap_coordinates
+{
+  swap_coordinates(const unsigned int i, const unsigned int j)
+    : i(i)
+    , j(j)
+  {}
+
+  Point<3>
+  operator()(const Point<3> &p_init) const
+  {
+    Point<3> p = p_init;
+    std::swap(p[i], p[j]);
+    return p;
+  }
+
+private:
+  const unsigned int i;
+  const unsigned int j;
+};
+
 void
 test_intersection_inside(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
 {
   GridGenerator::hyper_cube(tria0, -1., 1.);
   GridGenerator::hyper_cube(tria1, -0.4, 0.2);
-  const double expected_measure = 0.36;
-  const auto   vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                         tria1.begin_active(),
-                                                         MappingQ1<3>(),
-                                                         MappingQ1<2, 3>());
-
-  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(sum - expected_measure) < 1e-15);
+
+  {
+    const double expected_measure = 0.36;
+    const auto   vec_of_arrays =
+      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
+                                                           tria1.begin_active(),
+                                                           MappingQ1<3>(),
+                                                           MappingQ1<2, 3>());
+
+    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(sum - expected_measure) < 1e-15);
+  }
+
+  {
+    GridTools::transform(swap_coordinates(1, 2), tria1);
+    const double expected_measure = 0.36;
+    const auto   vec_of_arrays =
+      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
+                                                           tria1.begin_active(),
+                                                           MappingQ1<3>(),
+                                                           MappingQ1<2, 3>());
+
+    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(sum - expected_measure) < 1e-15);
+  }
+
+  {
+    GridTools::transform(swap_coordinates(0, 1), tria1);
+    const double expected_measure = 0.36;
+    const auto   vec_of_arrays =
+      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
+                                                           tria1.begin_active(),
+                                                           MappingQ1<3>(),
+                                                           MappingQ1<2, 3>());
+
+    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(sum - expected_measure) < 1e-15);
+  }
   deallog << "OK" << std::endl;
 }
 
@@ -53,17 +108,50 @@ test_intersection(Triangulation<3> &tria0, Triangulation<2, 3> &tria1)
 {
   GridGenerator::hyper_cube(tria0, -1., 1.);
   GridGenerator::hyper_cube(tria1, .5, 1.5);
-  const double expected_measure = 0.25;
-  const auto   vec_of_arrays =
-    CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
-                                                         tria1.begin_active(),
-                                                         MappingQ1<3>(),
-                                                         MappingQ1<2, 3>());
-
-  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(sum - expected_measure) < 1e-15);
+
+  {
+    const double expected_measure = 0.25;
+    const auto   vec_of_arrays =
+      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
+                                                           tria1.begin_active(),
+                                                           MappingQ1<3>(),
+                                                           MappingQ1<2, 3>());
+
+    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(sum - expected_measure) < 1e-15);
+  }
+
+  {
+    GridTools::transform(swap_coordinates(1, 2), tria1);
+    const double expected_measure = 0.25;
+    const auto   vec_of_arrays =
+      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
+                                                           tria1.begin_active(),
+                                                           MappingQ1<3>(),
+                                                           MappingQ1<2, 3>());
+
+    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(sum - expected_measure) < 1e-15);
+  }
+
+  {
+    GridTools::transform(swap_coordinates(0, 1), tria1);
+    const double expected_measure = 0.25;
+    const auto   vec_of_arrays =
+      CGALWrappers::compute_intersection_of_cells<3, 2, 3>(tria0.begin_active(),
+                                                           tria1.begin_active(),
+                                                           MappingQ1<3>(),
+                                                           MappingQ1<2, 3>());
+
+    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(sum - expected_measure) < 1e-15);
+  }
   deallog << "OK" << std::endl;
 }
 

diff --git a/tests/cgal/cgal_intersect_simplices_3d_3d.cc b/tests/cgal/cgal_intersect_simplices_3d_3d.cc
index fd61a9f31918928be756952c832c3ece749b677b..639d30a128042b3cfdd9bb71282870d6db8a4e44 100644 (file)
--- a/tests/cgal/cgal_intersect_simplices_3d_3d.cc
+++ b/tests/cgal/cgal_intersect_simplices_3d_3d.cc
@@ -17,6 +17,8 @@
 
 #include <deal.II/base/quadrature_lib.h>
 
+#include <deal.II/fe/mapping_q1.h>
+
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>