From: Johannes Heinz Date: Fri, 2 Sep 2022 11:46:54 +0000 (+0200) Subject: fix 2d/3d if not in xy plane X-Git-Tag: v9.5.0-rc1~940^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=021ae5403f4f751ae1110848b06b80efcdfa7433;p=dealii.git fix 2d/3d if not in xy plane --- diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc index d9b618c784..9537cae745 100644 --- a/source/cgal/intersections.cc +++ b/source/cgal/intersections.cc @@ -77,8 +77,6 @@ using Triangulation2 = CGAL::Triangulation_2; using Triangulation3 = CGAL::Triangulation_3; using Triangulation3_exact = CGAL::Triangulation_3; using Triangulation3_inexact = CGAL::Triangulation_3; -typedef CGAL::Projection_traits_xy_3 Gt; -typedef CGAL::Delaunay_triangulation_2 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 *vps = boost::get>(&*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, 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 435bbcab08..5d01660fdf 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_2d.cc +++ b/tests/cgal/cgal_intersect_simplices_1d_2d.cc @@ -19,6 +19,8 @@ #include +#include + #include #include #include diff --git a/tests/cgal/cgal_intersect_simplices_1d_3d.cc b/tests/cgal/cgal_intersect_simplices_1d_3d.cc index bfd00f9d7a..69c8505f7e 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_1d_3d.cc @@ -18,6 +18,8 @@ #include +#include + #include #include #include diff --git a/tests/cgal/cgal_intersect_simplices_2d_2d.cc b/tests/cgal/cgal_intersect_simplices_2d_2d.cc index d2d0cbf36b..0e9bf03bba 100644 --- a/tests/cgal/cgal_intersect_simplices_2d_2d.cc +++ b/tests/cgal/cgal_intersect_simplices_2d_2d.cc @@ -19,6 +19,8 @@ #include +#include + #include #include #include diff --git a/tests/cgal/cgal_intersect_simplices_2d_3d.cc b/tests/cgal/cgal_intersect_simplices_2d_3d.cc index 38070da4a6..be43f35fd7 100644 --- a/tests/cgal/cgal_intersect_simplices_2d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_2d_3d.cc @@ -18,6 +18,8 @@ #include +#include + #include #include #include @@ -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 fd61a9f319..639d30a128 100644 --- a/tests/cgal/cgal_intersect_simplices_3d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_3d_3d.cc @@ -17,6 +17,8 @@ #include +#include + #include #include #include