From: Johannes Heinz Date: Wed, 19 Oct 2022 07:47:18 +0000 (+0200) Subject: Fix linker error, due to compute_intersection_of_cells() X-Git-Tag: v9.5.0-rc1~460^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14362%2Fhead;p=dealii.git Fix linker error, due to compute_intersection_of_cells() Co-authored-by: Marco Feder --- diff --git a/include/deal.II/cgal/intersections.h b/include/deal.II/cgal/intersections.h index 0435b04af9..1162a25ec9 100644 --- a/include/deal.II/cgal/intersections.h +++ b/include/deal.II/cgal/intersections.h @@ -18,6 +18,7 @@ #include +#include #include #include @@ -60,21 +61,20 @@ namespace CGALWrappers /** - * Same function as above, but working directly with vertices. + * Same function as above, but working directly with vertices. It's assumed + * that the first vertices are coming from a Triangulation, + * while the other vertices are from a Triangulation, with + * @p dim0 > @p dim1. + * + * @note The vertices have to be given in CGAL order. */ - template + template std::vector, dim1 + 1>> compute_intersection_of_cells( - const std::array, n_vertices0> &vertices0, - const std::array, n_vertices1> &vertices1, - const double tol = 1e-9) - { - (void)vertices0; - (void)vertices1; - (void)tol; - Assert(false, ExcMessage("No explicit template instantiation available")); - return {}; - } + const ArrayView> &vertices0, + const ArrayView> &vertices1, + const double tol = 1e-9); + } // namespace CGALWrappers DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index b6f4634012..041b10c7b8 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -601,28 +601,28 @@ namespace CGALWrappers * @param mapping Mapping object for the cell. * @return Array of vertices in CGAL order. */ - template - std::array, n_vertices> + template + std::vector> 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, + const unsigned int n_vertices = cell->n_vertices(); + Assert((n_vertices == ReferenceCells::get_hypercube().n_vertices()) || + (n_vertices == ReferenceCells::get_simplex().n_vertices()), ExcNotImplemented()); - AssertDimension(mapping.get_vertices(cell).size(), n_vertices); - - std::array, n_vertices> vertices; + std::vector> ordered_vertices(n_vertices); std::copy_n(mapping.get_vertices(cell).begin(), - vertices.size(), - vertices.begin()); + n_vertices, + ordered_vertices.begin()); if (ReferenceCell::n_vertices_to_type(dim, n_vertices) == ReferenceCells::Quadrilateral) - std::swap(vertices[2], vertices[3]); + std::swap(ordered_vertices[2], ordered_vertices[3]); - return vertices; + return ordered_vertices; } /** diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc index a1634706f3..245103a748 100644 --- a/source/cgal/intersections.cc +++ b/source/cgal/intersections.cc @@ -229,126 +229,127 @@ namespace CGALWrappers CGALSegment2, CGALTriangle2, std::vector>> - compute_intersection(const std::array, 3> &first_simplex, - const std::array, 3> &second_simplex) + compute_intersection_triangle_triangle( + const ArrayView> &triangle0, + const ArrayView> &triangle1) { + AssertDimension(triangle0.size(), 3); + AssertDimension(triangle0.size(), triangle1.size()); + std::array pts0, pts1; - std::transform( - first_simplex.begin(), - first_simplex.end(), - pts0.begin(), - [&](const Point<2> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - std::transform( - second_simplex.begin(), - second_simplex.end(), - pts1.begin(), - [&](const Point<2> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - CGALTriangle2 triangle1{pts0[0], pts0[1], pts0[2]}; - CGALTriangle2 triangle2{pts1[0], pts1[1], pts1[2]}; - return convert_boost_to_std(CGAL::intersection(triangle1, triangle2)); - } + std::transform(triangle0.begin(), + triangle0.end(), + pts0.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + std::transform(triangle1.begin(), + triangle1.end(), + pts1.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + CGALTriangle2 cgal_triangle0{pts0[0], pts0[1], pts0[2]}; + CGALTriangle2 cgal_triangle1{pts1[0], pts1[1], pts1[2]}; + return convert_boost_to_std( + CGAL::intersection(cgal_triangle0, cgal_triangle1)); + } std_cxx17::optional> - compute_intersection(const std::array, 3> &first_simplex, - const std::array, 2> &second_simplex) + compute_intersection_triangle_segment( + const ArrayView> &triangle, + const ArrayView> &segment) { + AssertDimension(triangle.size(), 3); + AssertDimension(segment.size(), 2); + std::array pts0; std::array pts1; - std::transform( - first_simplex.begin(), - first_simplex.end(), - pts0.begin(), - [&](const Point<2> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - std::transform( - second_simplex.begin(), - second_simplex.end(), - pts1.begin(), - [&](const Point<2> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - CGALTriangle2 triangle{pts0[0], pts0[1], pts0[2]}; - CGALSegment2 segm{pts1[0], pts1[1]}; - return convert_boost_to_std(CGAL::intersection(segm, triangle)); + std::transform(triangle.begin(), + triangle.end(), + pts0.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + std::transform(segment.begin(), + segment.end(), + pts1.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + CGALTriangle2 cgal_triangle{pts0[0], pts0[1], pts0[2]}; + CGALSegment2 cgal_segment{pts1[0], pts1[1]}; + return convert_boost_to_std( + CGAL::intersection(cgal_segment, cgal_triangle)); } // rectangle-rectangle std::vector - compute_intersection(const std::array, 4> &first_simplex, - const std::array, 4> &second_simplex) + compute_intersection_rect_rect(const ArrayView> &rectangle0, + const ArrayView> &rectangle1) { + AssertDimension(rectangle0.size(), 4); + AssertDimension(rectangle0.size(), rectangle1.size()); + std::array pts0, pts1; - std::transform( - first_simplex.begin(), - first_simplex.end(), - pts0.begin(), - [&](const Point<2> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - std::transform( - second_simplex.begin(), - second_simplex.end(), - pts1.begin(), - [&](const Point<2> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - const CGALPolygon first{pts0.begin(), pts0.end()}; - const CGALPolygon second{pts1.begin(), pts1.end()}; + + std::transform(rectangle0.begin(), + rectangle0.end(), + pts0.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + std::transform(rectangle1.begin(), + rectangle1.end(), + pts1.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + const CGALPolygon first_poly{pts0.begin(), pts0.end()}; + const CGALPolygon second_poly{pts1.begin(), pts1.end()}; std::vector poly_list; - CGAL::intersection(first, second, std::back_inserter(poly_list)); + CGAL::intersection(first_poly, + second_poly, + std::back_inserter(poly_list)); return poly_list; } std_cxx17::optional> - compute_intersection(const std::array, 2> &first_simplex, - const std::array, 4> &second_simplex) + compute_intersection_tetra_segment( + const ArrayView> &tetrahedron, + const ArrayView> &segment) { # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0) + + AssertDimension(tetrahedron.size(), 4); + AssertDimension(segment.size(), 2); + std::array pts0; std::array pts1; - std::transform( - first_simplex.begin(), - first_simplex.end(), - pts0.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); + std::transform(tetrahedron.begin(), + tetrahedron.end(), + pts0.begin(), + &CGALWrappers::dealii_point_to_cgal_point); - std::transform( - second_simplex.begin(), - second_simplex.end(), - pts1.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - CGALTetra tetra{pts0[0], pts0[1], pts0[2], pts0[3]}; - CGALSegment3 segm{pts1[0], pts1[1]}; - return convert_boost_to_std(CGAL::intersection(segm, tetra)); + std::transform(segment.begin(), + segment.end(), + pts1.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]}; + CGALSegment3 cgal_segment{pts1[0], pts1[1]}; + return convert_boost_to_std( + CGAL::intersection(cgal_segment, cgal_tetrahedron)); # else Assert( false, ExcMessage( "This function requires a version of CGAL greater or equal than 5.5.")); - (void)first_simplex; - (void)second_simplex; + (void)tetrahedron; + (void)segment; return {}; # endif } @@ -359,400 +360,473 @@ namespace CGALWrappers CGALSegment3, CGALTriangle3, std::vector>> - compute_intersection(const std::array, 3> &first_simplex, - const std::array, 4> &second_simplex) + compute_intersection_tetra_triangle( + const ArrayView> &tetrahedron, + const ArrayView> &triangle) { # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0) + + AssertDimension(tetrahedron.size(), 4); + AssertDimension(triangle.size(), 3); + std::array pts0; std::array pts1; - std::transform( - first_simplex.begin(), - first_simplex.end(), - pts0.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - std::transform( - second_simplex.begin(), - second_simplex.end(), - pts1.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - CGALTetra tetra{pts0[0], pts0[1], pts0[2], pts0[3]}; - CGALTriangle3 triangle{pts1[0], pts1[1], pts1[2]}; - return convert_boost_to_std(CGAL::intersection(triangle, tetra)); + std::transform(tetrahedron.begin(), + tetrahedron.end(), + pts0.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + std::transform(triangle.begin(), + triangle.end(), + pts1.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + CGALTetra cgal_tetrahedron{pts0[0], pts0[1], pts0[2], pts0[3]}; + CGALTriangle3 cgal_triangle{pts1[0], pts1[1], pts1[2]}; + return convert_boost_to_std( + CGAL::intersection(cgal_triangle, cgal_tetrahedron)); # else Assert( false, ExcMessage( "This function requires a version of CGAL greater or equal than 5.5.")); - (void)first_simplex; - (void)second_simplex; + (void)tetrahedron; + (void)triangle; return {}; # endif } - } // namespace internal + // quad-quad + std::vector, 3>> + compute_intersection_quad_quad(const ArrayView> &quad0, + const ArrayView> &quad1, + const double tol) + { + AssertDimension(quad0.size(), 4); + AssertDimension(quad0.size(), quad1.size()); + const auto intersection_test = + internal::compute_intersection_rect_rect(quad0, quad1); - // Specialization for quads - template <> - std::vector, 3>> - compute_intersection_of_cells<2, 2, 2, 4, 4>( - const std::array, 4> &vertices0, - const std::array, 4> &vertices1, - const double tol) - { - const auto intersection_test = - internal::compute_intersection(vertices0, vertices1); + if (!intersection_test.empty()) + { + const auto & poly = intersection_test[0].outer_boundary(); + const unsigned int size_poly = poly.size(); + if (size_poly == 3) + { + // intersection is a triangle itself, so directly return its + // vertices. + return { + {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)), + CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)), + CGALWrappers::cgal_point_to_dealii_point<2>( + poly.vertex(2))}}}; + } + else if (size_poly >= 4) + { + // intersection is a polygon, need to triangulate it. + std::vector, 3>> collection; - if (!intersection_test.empty()) - { - const auto & poly = intersection_test[0].outer_boundary(); - const unsigned int size_poly = poly.size(); - if (size_poly == 3) - { - // intersection is a triangle itself, so directly return its - // vertices. - return { - {{CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(0)), - CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(1)), - CGALWrappers::cgal_point_to_dealii_point<2>(poly.vertex(2))}}}; - } - else if (size_poly >= 4) - { - // intersection is a polygon, need to triangulate it. - std::vector, 3>> collection; - - CDT cdt; - cdt.insert_constraint(poly.vertices_begin(), - poly.vertices_end(), - true); - - internal::mark_domains(cdt); - std::array, 3> vertices; - - for (const Face_handle f : cdt.finite_face_handles()) - { - if (f->info().in_domain() && - CGAL::to_double(cdt.triangle(f).area()) > tol) - { - collection.push_back( - {{CGALWrappers::cgal_point_to_dealii_point<2>( - cdt.triangle(f).vertex(0)), - CGALWrappers::cgal_point_to_dealii_point<2>( - cdt.triangle(f).vertex(1)), - CGALWrappers::cgal_point_to_dealii_point<2>( - cdt.triangle(f).vertex(2))}}); - } - } - return collection; - } - else - { - Assert(false, ExcMessage("The polygon is degenerate.")); - return {}; - } - } - else - { - return {}; - } - } + CDT cdt; + cdt.insert_constraint(poly.vertices_begin(), + poly.vertices_end(), + true); + internal::mark_domains(cdt); + std::array, 3> vertices; + for (Face_handle f : cdt.finite_face_handles()) + { + if (f->info().in_domain() && + CGAL::to_double(cdt.triangle(f).area()) > tol) + { + collection.push_back( + {{CGALWrappers::cgal_point_to_dealii_point<2>( + cdt.triangle(f).vertex(0)), + CGALWrappers::cgal_point_to_dealii_point<2>( + cdt.triangle(f).vertex(1)), + CGALWrappers::cgal_point_to_dealii_point<2>( + cdt.triangle(f).vertex(2))}}); + } + } + return collection; + } + else + { + Assert(false, ExcMessage("The polygon is degenerate.")); + return {}; + } + } + else + { + return {}; + } + } - // Specialization for quad \cap line - template <> - std::vector, 2>> - compute_intersection_of_cells<2, 1, 2, 4, 2>( - const std::array, 4> &vertices0, - const std::array, 2> &vertices1, - const double tol) - { - std::array pts; - std::transform( - vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<2> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - CGALPolygon poly(pts.begin(), pts.end()); - - CGALSegment2 segm( - CGALWrappers::dealii_point_to_cgal_point(vertices1[0]), - CGALWrappers::dealii_point_to_cgal_point(vertices1[1])); - CDT cdt; - cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true); - std::vector, 2>> vertices; - internal::mark_domains(cdt); - for (Face_handle f : cdt.finite_face_handles()) - { - if (f->info().in_domain() && - CGAL::to_double(cdt.triangle(f).area()) > tol && - CGAL::do_intersect(segm, cdt.triangle(f))) - { - const auto intersection = CGAL::intersection(segm, cdt.triangle(f)); - if (const CGALSegment2 *s = - boost::get(&*intersection)) - { - vertices.push_back( - {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]), - CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}}); - } - } - } - return vertices; - } + // Specialization for quad \cap line + std::vector, 2>> + compute_intersection_quad_line(const ArrayView> &quad, + const ArrayView> &line, + const double tol) + { + AssertDimension(quad.size(), 4); + AssertDimension(line.size(), 2); + + std::array pts; + + std::transform(quad.begin(), + quad.end(), + pts.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + CGALPolygon poly(pts.begin(), pts.end()); + + CGALSegment2 segm( + CGALWrappers::dealii_point_to_cgal_point(line[0]), + CGALWrappers::dealii_point_to_cgal_point(line[1])); + CDT cdt; + cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true); + std::vector, 2>> vertices; + internal::mark_domains(cdt); + for (Face_handle f : cdt.finite_face_handles()) + { + if (f->info().in_domain() && + CGAL::to_double(cdt.triangle(f).area()) > tol && + CGAL::do_intersect(segm, cdt.triangle(f))) + { + const auto intersection = + CGAL::intersection(segm, cdt.triangle(f)); + if (const CGALSegment2 *s = + boost::get(&*intersection)) + { + vertices.push_back( + {{CGALWrappers::cgal_point_to_dealii_point<2>((*s)[0]), + CGALWrappers::cgal_point_to_dealii_point<2>((*s)[1])}}); + } + } + } - // specialization for hex \cap line - template <> - std::vector, 2>> - compute_intersection_of_cells<3, 1, 3, 8, 2>( - const std::array, 8> &vertices0, - const std::array, 2> &vertices1, - const double tol) - { + return vertices; + } + + // specialization for hex \cap line + std::vector, 2>> + compute_intersection_hexa_line(const ArrayView> &hexa, + const ArrayView> &line, + const double tol) + { # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0) - std::array pts; - std::transform( - vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - CGALSegment3_exact segm( - CGALWrappers::dealii_point_to_cgal_point(vertices1[0]), - CGALWrappers::dealii_point_to_cgal_point(vertices1[1])); - - // Subdivide the hex into tetrahedrons, and intersect each one of them with - // the line - std::vector, 2>> vertices; - Triangulation3_exact tria; - tria.insert(pts.begin(), pts.end()); - for (const auto &c : tria.finite_cell_handles()) - { - const auto &tet = tria.tetrahedron(c); - if (CGAL::do_intersect(segm, tet)) - { - const auto intersection = CGAL::intersection(segm, tet); - if (const CGALSegment3_exact *s = - boost::get(&*intersection)) - { - if (s->squared_length() > tol * tol) - { - vertices.push_back( - {{CGALWrappers::cgal_point_to_dealii_point<3>( - s->vertex(0)), - CGALWrappers::cgal_point_to_dealii_point<3>( - s->vertex(1))}}); - } - } - } - } - return vertices; + AssertDimension(hexa.size(), 8); + AssertDimension(line.size(), 2); + + std::array pts; + + std::transform( + hexa.begin(), + hexa.end(), + pts.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + CGALSegment3_exact cgal_segment( + CGALWrappers::dealii_point_to_cgal_point(line[0]), + CGALWrappers::dealii_point_to_cgal_point(line[1])); + + // Subdivide the hex into tetrahedrons, and intersect each one of them + // with the line + std::vector, 2>> vertices; + Triangulation3_exact cgal_triangulation; + cgal_triangulation.insert(pts.begin(), pts.end()); + for (const auto &c : cgal_triangulation.finite_cell_handles()) + { + const auto &cgal_tetrahedron = cgal_triangulation.tetrahedron(c); + if (CGAL::do_intersect(cgal_segment, cgal_tetrahedron)) + { + const auto intersection = + CGAL::intersection(cgal_segment, cgal_tetrahedron); + if (const CGALSegment3_exact *s = + boost::get(&*intersection)) + { + if (s->squared_length() > tol * tol) + { + vertices.push_back( + {{CGALWrappers::cgal_point_to_dealii_point<3>( + s->vertex(0)), + CGALWrappers::cgal_point_to_dealii_point<3>( + s->vertex(1))}}); + } + } + } + } + return vertices; # else - Assert( - false, - ExcMessage( - "This function requires a version of CGAL greater or equal than 5.5.")); - (void)vertices0; - (void)vertices1; - (void)tol; - return {}; + Assert( + false, + ExcMessage( + "This function requires a version of CGAL greater or equal than 5.5.")); + (void)hexa; + (void)line; + (void)tol; + return {}; # endif - } + } - template <> - std::vector, 3>> - compute_intersection_of_cells<3, 2, 3, 8, 4>( - const std::array, 8> &vertices0, - const std::array, 4> &vertices1, - const double tol) - { + std::vector, 3>> + compute_intersection_hexa_quad(const ArrayView> &hexa, + const ArrayView> &quad, + const double tol) + { # if DEAL_II_CGAL_VERSION_GTE(5, 5, 0) - std::array pts_hex; - std::array pts_quad; - std::transform( - vertices0.begin(), - vertices0.end(), - pts_hex.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - std::transform( - vertices1.begin(), - vertices1.end(), - pts_quad.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - // Subdivide hex into tetrahedrons - std::vector, 3>> vertices; - Triangulation3_exact tria; - tria.insert(pts_hex.begin(), pts_hex.end()); - - // Subdivide quad into triangles - 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_facets()) - { - if (CGAL::do_intersect(tet, tria_quad.triangle(f))) - { - const auto intersection = - CGAL::intersection(tria_quad.triangle(f), tet); - - if (const CGALTriangle3_exact *t = - boost::get(&*intersection)) - { - if (CGAL::to_double(t->squared_area()) > tol * tol) - { - vertices.push_back( - {{cgal_point_to_dealii_point<3>((*t)[0]), - cgal_point_to_dealii_point<3>((*t)[1]), - cgal_point_to_dealii_point<3>((*t)[2])}}); - } - } - - if (const std::vector *vps = - boost::get>(&*intersection)) - { - 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) - { - 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>( - triangle[0]), - CGALWrappers::cgal_point_to_dealii_point<3>( - triangle[1]), - CGALWrappers::cgal_point_to_dealii_point<3>( - triangle[2])}}; - - vertices.push_back(verts); - } - } - } - } - } - } + AssertDimension(hexa.size(), 8); + AssertDimension(quad.size(), 4); + + std::array pts_hex; + std::array pts_quad; + + std::transform( + hexa.begin(), + hexa.end(), + pts_hex.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + std::transform( + quad.begin(), + quad.end(), + pts_quad.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + // Subdivide hex into tetrahedrons + std::vector, 3>> vertices; + Triangulation3_exact triangulation_hexa; + triangulation_hexa.insert(pts_hex.begin(), pts_hex.end()); + + // Subdivide quad into triangles + Triangulation3_exact triangulation_quad; + triangulation_quad.insert(pts_quad.begin(), pts_quad.end()); + + for (const auto &c : triangulation_hexa.finite_cell_handles()) + { + const auto &tet = triangulation_hexa.tetrahedron(c); + + for (const auto &f : triangulation_quad.finite_facets()) + { + if (CGAL::do_intersect(tet, triangulation_quad.triangle(f))) + { + const auto intersection = + CGAL::intersection(triangulation_quad.triangle(f), tet); + + if (const CGALTriangle3_exact *t = + boost::get(&*intersection)) + { + if (CGAL::to_double(t->squared_area()) > tol * tol) + { + vertices.push_back( + {{cgal_point_to_dealii_point<3>((*t)[0]), + cgal_point_to_dealii_point<3>((*t)[1]), + cgal_point_to_dealii_point<3>((*t)[2])}}); + } + } + + if (const std::vector *vps = + boost::get>( + &*intersection)) + { + 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) + { + 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>( + triangle[0]), + CGALWrappers::cgal_point_to_dealii_point<3>( + triangle[1]), + CGALWrappers::cgal_point_to_dealii_point<3>( + triangle[2])}}; + + vertices.push_back(verts); + } + } + } + } + } + } - return vertices; + return vertices; # else - Assert( - false, - ExcMessage( - "This function requires a version of CGAL greater or equal than 5.5.")); - (void)vertices0; - (void)vertices1; - (void)tol; - return {}; + Assert( + false, + ExcMessage( + "This function requires a version of CGAL greater or equal than 5.5.")); + (void)hexa; + (void)quad; + (void)tol; + return {}; # endif - } + } + + std::vector, 4>> + compute_intersection_hexa_hexa(const ArrayView> &hexa0, + const ArrayView> &hexa1, + const double tol) + { + AssertDimension(hexa0.size(), 8); + AssertDimension(hexa0.size(), hexa1.size()); + std::array pts_hex0; + std::array pts_hex1; + std::transform( + hexa0.begin(), + hexa0.end(), + pts_hex0.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + std::transform( + hexa1.begin(), + hexa1.end(), + pts_hex1.begin(), + &CGALWrappers::dealii_point_to_cgal_point); + + Surface_mesh surf0, surf1, sm; + // Subdivide hex into tetrahedrons + std::vector, 4>> vertices; + Triangulation3_inexact tria0, tria1; + + tria0.insert(pts_hex0.begin(), pts_hex0.end()); + tria1.insert(pts_hex1.begin(), pts_hex1.end()); + + for (const auto &c0 : tria0.finite_cell_handles()) + { + const auto &tet0 = tria1.tetrahedron(c0); + const auto &tetg0 = CGAL::make_tetrahedron(tet0.vertex(0), + tet0.vertex(1), + tet0.vertex(2), + tet0.vertex(3), + surf0); + (void)tetg0; // instead of C++ 17s [[maybe unused]] + for (const auto &c1 : tria1.finite_cell_handles()) + { + const auto &tet1 = tria1.tetrahedron(c1); + const auto &tetg1 = CGAL::make_tetrahedron(tet1.vertex(0), + tet1.vertex(1), + tet1.vertex(2), + tet1.vertex(3), + surf1); + (void)tetg1; // instead of C++ 17s [[maybe unused]] + namespace PMP = CGAL::Polygon_mesh_processing; + const bool test_intersection = + PMP::corefine_and_compute_intersection(surf0, surf1, sm); + if (PMP::volume(sm) > tol && test_intersection) + { + // Collect tetrahedrons + Triangulation3_inexact triangulation_hexa; + triangulation_hexa.insert(sm.points().begin(), + sm.points().end()); + for (const auto &c : triangulation_hexa.finite_cell_handles()) + { + const auto &tet = triangulation_hexa.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))}}); + } + } + surf1.clear(); + sm.clear(); + } + surf0.clear(); + } + return vertices; + } + + } // namespace internal - template <> - std::vector, 4>> - compute_intersection_of_cells<3, 3, 3, 8, 8>( - const std::array, 8> &vertices0, - const std::array, 8> &vertices1, - const double tol) + template + std::vector, dim1 + 1>> + compute_intersection_of_cells( + const ArrayView> &vertices0, + const ArrayView> &vertices1, + const double tol) { - std::array pts_hex0; - std::array pts_hex1; - std::transform( - vertices0.begin(), - vertices0.end(), - pts_hex0.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - std::transform( - vertices1.begin(), - vertices1.end(), - pts_hex1.begin(), - [&](const Point<3> &p) { - return CGALWrappers::dealii_point_to_cgal_point(p); - }); - - - Surface_mesh surf0, surf1, sm; - // Subdivide hex into tetrahedrons - std::vector, 4>> vertices; - Triangulation3_inexact tria0, tria1; - - tria0.insert(pts_hex0.begin(), pts_hex0.end()); - tria1.insert(pts_hex1.begin(), pts_hex1.end()); - - for (const auto &c0 : tria0.finite_cell_handles()) + const unsigned int n_vertices0 = vertices0.size(); + const unsigned int n_vertices1 = vertices1.size(); + + Assert( + n_vertices0 > 0 || n_vertices1 > 0, + ExcMessage( + "The intersection cannot be computed as at least one of the two cells has no vertices.")); + + if constexpr (dim0 == 2 && dim1 == 2 && spacedim == 2) { - const auto & tet0 = tria1.tetrahedron(c0); - [[maybe_unused]] const auto &tetg0 = - CGAL::make_tetrahedron(tet0.vertex(0), - tet0.vertex(1), - tet0.vertex(2), - tet0.vertex(3), - surf0); - for (const auto &c1 : tria1.finite_cell_handles()) + if (n_vertices0 == 4 && n_vertices1 == 4) { - const auto & tet1 = tria1.tetrahedron(c1); - [[maybe_unused]] const auto &tetg1 = - CGAL::make_tetrahedron(tet1.vertex(0), - tet1.vertex(1), - tet1.vertex(2), - tet1.vertex(3), - surf1); - namespace PMP = CGAL::Polygon_mesh_processing; - const bool test_intersection = - PMP::corefine_and_compute_intersection(surf0, surf1, sm); - if (PMP::volume(sm) > tol && test_intersection) - { - // 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))}}); - } - } - surf1.clear(); - sm.clear(); + return internal::compute_intersection_quad_quad(vertices0, + vertices1, + tol); } - surf0.clear(); } - return vertices; + else if constexpr (dim0 == 2 && dim1 == 1 && spacedim == 2) + { + if (n_vertices0 == 4 && n_vertices1 == 2) + { + return internal::compute_intersection_quad_line(vertices0, + vertices1, + tol); + } + } + else if constexpr (dim0 == 3 && dim1 == 1 && spacedim == 3) + { + if (n_vertices0 == 8 && n_vertices1 == 2) + { + return internal::compute_intersection_hexa_line(vertices0, + vertices1, + tol); + } + } + else if constexpr (dim0 == 3 && dim1 == 2 && spacedim == 3) + { + if (n_vertices0 == 8 && n_vertices1 == 4) + { + return internal::compute_intersection_hexa_quad(vertices0, + vertices1, + tol); + } + } + else if constexpr (dim0 == 3 && dim1 == 3 && spacedim == 3) + { + if (n_vertices0 == 8 && n_vertices1 == 8) + { + return internal::compute_intersection_hexa_hexa(vertices0, + vertices1, + tol); + } + } + else + { + Assert(false, ExcNotImplemented()); + return {}; + } + (void)tol; + return {}; } - template std::vector, dim1 + 1>> compute_intersection_of_cells( @@ -762,25 +836,21 @@ namespace CGALWrappers const Mapping & mapping1, const double tol) { - Assert(mapping0.get_vertices(cell0).size() == std::pow(2, dim0), + Assert(mapping0.get_vertices(cell0).size() == + ReferenceCells::get_hypercube().n_vertices(), ExcNotImplemented()); - Assert(mapping1.get_vertices(cell1).size() == std::pow(2, dim1), + Assert(mapping1.get_vertices(cell1).size() == + ReferenceCells::get_hypercube().n_vertices(), ExcNotImplemented()); - 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); + 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" diff --git a/source/cgal/intersections.inst.in b/source/cgal/intersections.inst.in index d550252e0f..3024147791 100644 --- a/source/cgal/intersections.inst.in +++ b/source/cgal/intersections.inst.in @@ -27,5 +27,11 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS) const Mapping & mapping1, const double tol); + template std::vector, dim1 + 1>> + compute_intersection_of_cells( + const ArrayView> &vertices0, + const ArrayView> &vertices1, + const double tol); + #endif } diff --git a/tests/cgal/cgal_intersect_simplices_1d_2d.cc b/tests/cgal/cgal_intersect_simplices_1d_2d.cc index b0d00ba33c..fdb33fa109 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_2d.cc +++ b/tests/cgal/cgal_intersect_simplices_1d_2d.cc @@ -42,11 +42,8 @@ test_inside_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1) const auto cell1 = tria1.begin_active(); // cell1->vertex(1)+=Point<2>(1.0,1.); - const auto vec_of_arrays = - CGALWrappers::compute_intersection_of_cells<2, 1, 2>(cell0, - cell1, - MappingQ1<2>(), - MappingQ1<1, 2>()); + const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<2>(), MappingQ1<1, 2>()); const auto quad = qgauss.mapped_quadrature(vec_of_arrays); const double sum = @@ -67,11 +64,8 @@ test_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<2>(), MappingQ1<1, 2>()); const auto quad = qgauss.mapped_quadrature(vec_of_arrays); const double sum = @@ -91,11 +85,8 @@ test_failing_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<2>(), MappingQ1<1, 2>()); const auto quad = qgauss.mapped_quadrature(vec_of_arrays); const double sum = diff --git a/tests/cgal/cgal_intersect_simplices_1d_3d.cc b/tests/cgal/cgal_intersect_simplices_1d_3d.cc index d7bd153ac4..7e64cc8426 100644 --- a/tests/cgal/cgal_intersect_simplices_1d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_1d_3d.cc @@ -41,10 +41,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<1, 3> &tria1) (tria1.begin_active()->vertex(1) - tria1.begin_active()->vertex(0)).norm(); const auto vec_of_arrays = - CGALWrappers::compute_intersection_of_cells<3, 1, 3>(tria0.begin_active(), - tria1.begin_active(), - MappingQ1<3>(), - MappingQ1<1, 3>()); + CGALWrappers::compute_intersection_of_cells(tria0.begin_active(), + tria1.begin_active(), + MappingQ1<3>(), + MappingQ1<1, 3>()); const auto quad = QGaussSimplex<1>(1).mapped_quadrature(vec_of_arrays); const double sum = @@ -65,11 +65,8 @@ test_intersection(Triangulation<3> &tria0, Triangulation<1, 3> &tria1) cell1->vertex(1) = Point<3>(1.5, 1.5, 1.5); const double expected_measure = std::sqrt(3.); - const auto vec_of_arrays = - CGALWrappers::compute_intersection_of_cells<3, 1, 3>(tria0.begin_active(), - cell1, - MappingQ1<3>(), - MappingQ1<1, 3>()); + const auto vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + tria0.begin_active(), cell1, MappingQ1<3>(), MappingQ1<1, 3>()); const auto quad = QGaussSimplex<1>(1).mapped_quadrature(vec_of_arrays); const double sum = @@ -90,11 +87,8 @@ test_failing_intersection(Triangulation<3> &tria0, Triangulation<1, 3> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<3>(), MappingQ1<1, 3>()); const auto quad = QGaussSimplex<1>(1).mapped_quadrature(vec_of_arrays); const double sum = diff --git a/tests/cgal/cgal_intersect_simplices_2d_2d.cc b/tests/cgal/cgal_intersect_simplices_2d_2d.cc index c4e99b60b6..1d77622d51 100644 --- a/tests/cgal/cgal_intersect_simplices_2d_2d.cc +++ b/tests/cgal/cgal_intersect_simplices_2d_2d.cc @@ -41,11 +41,8 @@ test_inside_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<2>(), MappingQ1<2>()); const auto quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays); @@ -67,11 +64,8 @@ test_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<2>(), MappingQ1<2>()); const auto quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays); const double sum = @@ -92,11 +86,8 @@ test_failing_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<2>(), MappingQ1<2>()); const auto quad = QGaussSimplex<2>(1).mapped_quadrature(vec_of_arrays); const double sum = diff --git a/tests/cgal/cgal_intersect_simplices_2d_3d.cc b/tests/cgal/cgal_intersect_simplices_2d_3d.cc index be43f35fd7..225637f71c 100644 --- a/tests/cgal/cgal_intersect_simplices_2d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_2d_3d.cc @@ -58,10 +58,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<2, 3> &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>()); + CGALWrappers::compute_intersection_of_cells(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 = @@ -73,10 +73,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<2, 3> &tria1) 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>()); + CGALWrappers::compute_intersection_of_cells(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 = @@ -88,10 +88,10 @@ test_intersection_inside(Triangulation<3> &tria0, Triangulation<2, 3> &tria1) 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>()); + CGALWrappers::compute_intersection_of_cells(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 = @@ -112,10 +112,10 @@ test_intersection(Triangulation<3> &tria0, Triangulation<2, 3> &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>()); + CGALWrappers::compute_intersection_of_cells(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 = @@ -127,10 +127,10 @@ test_intersection(Triangulation<3> &tria0, Triangulation<2, 3> &tria1) 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>()); + CGALWrappers::compute_intersection_of_cells(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 = @@ -142,10 +142,10 @@ test_intersection(Triangulation<3> &tria0, Triangulation<2, 3> &tria1) 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>()); + CGALWrappers::compute_intersection_of_cells(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 = diff --git a/tests/cgal/cgal_intersect_simplices_3d_3d.cc b/tests/cgal/cgal_intersect_simplices_3d_3d.cc index d19f3f3ef9..e181eba7f9 100644 --- a/tests/cgal/cgal_intersect_simplices_3d_3d.cc +++ b/tests/cgal/cgal_intersect_simplices_3d_3d.cc @@ -36,11 +36,8 @@ test_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<3>(), MappingQ1<3>()); const auto quad = QGaussSimplex<3>(1).mapped_quadrature(vec_of_arrays); const double sum = @@ -60,11 +57,8 @@ test_failing_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1) 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 vec_of_arrays = CGALWrappers::compute_intersection_of_cells( + cell0, cell1, MappingQ1<3>(), MappingQ1<3>()); const auto quad = QGaussSimplex<3>(1).mapped_quadrature(vec_of_arrays); const double sum = diff --git a/tests/cgal/cgal_intersect_simplices_with_vertices.cc b/tests/cgal/cgal_intersect_simplices_with_vertices.cc new file mode 100644 index 0000000000..a494e9103a --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_with_vertices.cc @@ -0,0 +1,136 @@ +// --------------------------------------------------------------------- + +// 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. +// +// --------------------------------------------------------------------- + +// This test is the same as cgal_intersection_simplices_2d_2d() but directly +// working with vertices While the functionality of the intersections is tested +// with the iterator interface this test ensures there are no linker errors if +// the vertex version is called directly. + +#include + +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace CGALWrappers; + +void +test_inside_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, -0.45, 0.45); + GridTools::rotate(numbers::PI_4, tria1); + const double expected_area = 0.81; + const auto cell0 = tria0.begin_active(); + const auto cell1 = tria1.begin_active(); + + const auto &vertices0 = + CGALWrappers::get_vertices_in_cgal_order(cell0, MappingQ1<2>()); + const auto &vertices1 = + CGALWrappers::get_vertices_in_cgal_order(cell1, MappingQ1<2>()); + + const auto &vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<2, 2, 2>(vertices0, vertices1); + + 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_intersection(Triangulation<2> &tria0, Triangulation<2> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, .5, 1.45); + const double expected_area = 0.25; + + const auto cell0 = tria0.begin_active(); + const auto cell1 = tria1.begin_active(); + + const auto &vertices0 = + CGALWrappers::get_vertices_in_cgal_order(cell0, MappingQ1<2>()); + const auto &vertices1 = + CGALWrappers::get_vertices_in_cgal_order(cell1, MappingQ1<2>()); + + const auto &vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<2, 2, 2>(vertices0, vertices1); + + + 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 &vertices0 = + CGALWrappers::get_vertices_in_cgal_order(cell0, MappingQ1<2>()); + const auto &vertices1 = + CGALWrappers::get_vertices_in_cgal_order(cell1, MappingQ1<2>()); + + const auto &vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<2, 2, 2>(vertices0, vertices1); + + 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; +} + + + +int +main() +{ + initlog(); + Triangulation<2> tria0; + Triangulation<2> tria1; + + test_inside_intersection(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_with_vertices.output b/tests/cgal/cgal_intersect_simplices_with_vertices.output new file mode 100644 index 0000000000..fb71de2867 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_with_vertices.output @@ -0,0 +1,4 @@ + +DEAL::OK +DEAL::OK +DEAL::OK