From: Marco Feder Date: Tue, 5 Jul 2022 16:55:37 +0000 (+0200) Subject: Provide dimension-independent intersection between deal.II cells X-Git-Tag: v9.5.0-rc1~940^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e49a3823654a44264bbc1796c7d60df3beae3d7f;p=dealii.git Provide dimension-independent intersection between deal.II cells --- diff --git a/cmake/modules/FindCGAL.cmake b/cmake/modules/FindCGAL.cmake index 20ae44e93c..f3aa3e9631 100644 --- a/cmake/modules/FindCGAL.cmake +++ b/cmake/modules/FindCGAL.cmake @@ -48,6 +48,10 @@ IF(DEAL_II_HAVE_CXX17 AND NOT FEATURE_BOOST_BUNDLED_CONFIGURED) 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) + ENDIF() # Make sure we dont' pass Boost::Boost over to deal.II. LIST(FILTER CGAL_LIBRARIES EXCLUDE REGEX "::") ENDIF() diff --git a/include/deal.II/cgal/intersections.h b/include/deal.II/cgal/intersections.h new file mode 100644 index 0000000000..96883afda5 --- /dev/null +++ b/include/deal.II/cgal/intersections.h @@ -0,0 +1,63 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_cgal_intersections_h +#define dealii_cgal_intersections_h + +#include + +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN +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. + * + * @note If the two cells have different intrinsic dimensions, than the + * first iterator is assumed to have with higher intrinsic dimension. + * The same holds for the two Mapping objects. + * + * + * @param cell0 Iterator to the first cell. + * @param cell1 Iterator to the second cell. + * @param mapping0 Mapping for the first cell. + * @param mapping1 Mapping for the second cell. + * @param tol Treshold to decide whether or not a simplex is included. + * @return Vector of arrays, where each array identify a simplex by its vertices. + */ + template + std::vector, dim1 + 1>> + compute_intersection_of_cells( + const typename Triangulation::cell_iterator &cell0, + const typename Triangulation::cell_iterator &cell1, + const Mapping & mapping0, + const Mapping & mapping1, + const double tol = 1e-9); +} // namespace CGALWrappers + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 93ecacacae..916a8e4437 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -476,7 +476,6 @@ namespace CGALWrappers CGALTriangulation tria; tria.insert(out_surface.points().begin(), out_surface.points().end()); - return compute_quadrature(tria, degree); } } diff --git a/source/cgal/CMakeLists.txt b/source/cgal/CMakeLists.txt index 0e93c1c761..2c951a6b36 100644 --- a/source/cgal/CMakeLists.txt +++ b/source/cgal/CMakeLists.txt @@ -18,7 +18,8 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src -surface_mesh.cc + surface_mesh.cc + intersections.cc ) diff --git a/source/cgal/intersections.cc b/source/cgal/intersections.cc new file mode 100644 index 0000000000..d9b618c784 --- /dev/null +++ b/source/cgal/intersections.cc @@ -0,0 +1,736 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include + + +#ifdef DEAL_II_WITH_CGAL + +# include + +# include + +# include + +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include + +# include +# include + +using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt; +using K_exact = CGAL::Exact_predicates_exact_constructions_kernel; +using K_inexact = CGAL::Exact_predicates_inexact_constructions_kernel; +using CGALPolygon = CGAL::Polygon_2; +using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; +using CGALTriangle2 = K::Triangle_2; +using CGALTriangle3 = K::Triangle_3; +using CGALTriangle3_exact = K_exact::Triangle_3; +using CGALPoint2 = K::Point_2; +using CGALPoint3 = K::Point_3; +using CGALPoint3_exact = K_exact::Point_3; +using CGALSegment2 = K::Segment_2; +using Surface_mesh = CGAL::Surface_mesh; +using CGALSegment3 = K::Segment_3; +using CGALSegment3_exact = K_exact::Segment_3; +using CGALTetra = K::Tetrahedron_3; +using CGALTetra_exact = K_exact::Tetrahedron_3; +using Triangulation2 = CGAL::Triangulation_2; +using Triangulation3 = CGAL::Triangulation_3; +using Triangulation3_exact = CGAL::Triangulation_3; +using Triangulation3_inexact = CGAL::Triangulation_3; +typedef CGAL::Projection_traits_xy_3 Gt; +typedef CGAL::Delaunay_triangulation_2 Delaunay; + +struct FaceInfo2 +{ + FaceInfo2() + {} + int nesting_level; + bool + in_domain() + { + return nesting_level % 2 == 1; + } +}; + +using Vb = CGAL::Triangulation_vertex_base_2; +using Fbb = CGAL::Triangulation_face_base_with_info_2; +using CFb = CGAL::Constrained_triangulation_face_base_2; +using Fb = CGAL::Delaunay_mesh_face_base_2; +using Tds = CGAL::Triangulation_data_structure_2; +using Itag = CGAL::Exact_predicates_tag; +using CDT = CGAL::Constrained_Delaunay_triangulation_2; +using Criteria = CGAL::Delaunay_mesh_size_criteria_2; +using Vertex_handle = CDT::Vertex_handle; +using Face_handle = CDT::Face_handle; + + + +DEAL_II_NAMESPACE_OPEN + +namespace CGALWrappers +{ + namespace internal + { + void + mark_domains(CDT & ct, + Face_handle start, + int index, + std::list &border) + { + if (start->info().nesting_level != -1) + { + return; + } + std::list queue; + queue.push_back(start); + while (!queue.empty()) + { + Face_handle fh = queue.front(); + queue.pop_front(); + if (fh->info().nesting_level == -1) + { + fh->info().nesting_level = index; + for (int i = 0; i < 3; i++) + { + CDT::Edge e(fh, i); + Face_handle n = fh->neighbor(i); + if (n->info().nesting_level == -1) + { + if (ct.is_constrained(e)) + border.push_back(e); + else + queue.push_back(n); + } + } + } + } + } + + + + void + mark_domains(CDT &cdt) + { + for (CDT::Face_handle f : cdt.all_face_handles()) + { + f->info().nesting_level = -1; + } + std::list border; + mark_domains(cdt, cdt.infinite_face(), 0, border); + while (!border.empty()) + { + CDT::Edge e = border.front(); + border.pop_front(); + Face_handle n = e.first->neighbor(e.second); + if (n->info().nesting_level == -1) + { + mark_domains(cdt, n, e.first->info().nesting_level + 1, border); + } + } + } + + // Collection of utilities that compute intersection between simplices + // identified by array of points. The return type is the one of + // CGAL::intersection(), i.e. a boost::optional>. + // Intersection between 2D and 3D objects and 1D/3D objects are available + // only with CGAL versions greater or equal than 5.1.5, hence the + // corresponding functions are guarded by #ifdef directives. All the + // signatures follow the convection that the first entity has an intrinsic + // dimension higher than the second one. + + boost::optional>> + compute_intersection(const std::array, 3> &first_simplex, + const std::array, 3> &second_simplex) + { + 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 CGAL::intersection(triangle1, triangle2); + } + + + + boost::optional> + compute_intersection(const std::array, 3> &first_simplex, + const std::array, 2> &second_simplex) + { + 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 CGAL::intersection(segm, triangle); + } + + + + // rectangle-rectangle + decltype(auto) + compute_intersection(const std::array, 4> &first_simplex, + const std::array, 4> &second_simplex) + { + 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::vector poly_list; + CGAL::intersection(first, second, std::back_inserter(poly_list)); + return poly_list; + } + + + + boost::optional> + compute_intersection(const std::array, 2> &first_simplex, + const std::array, 4> &second_simplex) + { +# if defined(CGAL_GEQ_515) + + 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]}; + CGALSegment3 segm{pts1[0], pts1[1]}; + return CGAL::intersection(segm, tetra); + +# else + + Assert( + false, + ExcMessage( + "This function requires a version of CGAL greater or equal than 5.1.5.")); + (void)first_simplex; + (void)second_simplex; + return {}; +# endif + } + + + + // tetra, triangle + boost::optional>> + compute_intersection(const std::array, 3> &first_simplex, + const std::array, 4> &second_simplex) + { +# if defined(CGAL_GEQ_515) + 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 CGAL::intersection(triangle, tetra); +# else + + Assert( + false, + ExcMessage( + "This function requires a version of CGAL greater or equal than 5.1.5.")); + (void)first_simplex; + (void)second_simplex; + return {}; +# endif + } + } // namespace internal + + + + // Specialization for quads + template <> + std::vector, 3>> + compute_intersection_of_cells<2, 2, 2>( + const typename Triangulation<2, 2>::cell_iterator &cell0, + const typename Triangulation<2, 2>::cell_iterator &cell1, + const Mapping<2, 2> & mapping0, + const Mapping<2, 2> & mapping1, + const double tol) + { + std::array, 4> vertices0, vertices1; + std::copy_n(mapping0.get_vertices(cell0).begin(), 4, vertices0.begin()); + std::copy_n(mapping1.get_vertices(cell1).begin(), 4, vertices1.begin()); + + std::swap(vertices0[2], vertices0[3]); + std::swap(vertices1[2], vertices1[3]); + 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; + + 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 + { + Assert(false, ExcMessage("Cells do not intersect.")); + return {}; + } + } + + + + // Specialization for quad \cap line + template <> + std::vector, 2>> + compute_intersection_of_cells<2, 1, 2>( + const typename Triangulation<2, 2>::cell_iterator &cell0, + const typename Triangulation<1, 2>::cell_iterator &cell1, + const Mapping<2, 2> & mapping0, + const Mapping<1, 2> & mapping1, + const double tol) + { + std::array, 4> vertices0; + std::array, 2> vertices1; + std::copy_n(mapping0.get_vertices(cell0).begin(), 4, vertices0.begin()); + std::copy_n(mapping1.get_vertices(cell1).begin(), 2, vertices1.begin()); + + std::swap(vertices0[2], vertices0[3]); + + 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 hex \cap line + template <> + std::vector, 2>> + compute_intersection_of_cells<3, 1, 3>( + const typename Triangulation<3, 3>::cell_iterator &cell0, + const typename Triangulation<1, 3>::cell_iterator &cell1, + const Mapping<3, 3> & mapping0, + const Mapping<1, 3> & mapping1, + const double tol) + { +# if defined(CGAL_GEQ_515) + (void)tol; // + std::array, 8> vertices0; // 8 vertices of the hex + std::array, 2> vertices1; // 2 endpoints of the segment + std::copy_n(mapping0.get_vertices(cell0).begin(), 8, vertices0.begin()); + std::copy_n(mapping1.get_vertices(cell1).begin(), 2, vertices1.begin()); + + 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; +# else + Assert( + false, + ExcMessage( + "This function requires a version of CGAL greater or equal than 5.1.5.")); + (void)cell0; + (void)cell1; + (void)mapping0; + (void)mapping1; + (void)tol; + return {}; + +# endif + } + + + + template <> + std::vector, 3>> + compute_intersection_of_cells<3, 2, 3>( + const typename Triangulation<3, 3>::cell_iterator &cell0, + const typename Triangulation<2, 3>::cell_iterator &cell1, + const Mapping<3, 3> & mapping0, + const Mapping<2, 3> & mapping1, + const double tol) + { +# if defined(CGAL_GEQ_515) + + std::array, 8> vertices0; // 8 vertices of the hex + std::array, 4> vertices1; // 4 vertices of the quad + std::copy_n(mapping0.get_vertices(cell0).begin(), 8, vertices0.begin()); + std::copy_n(mapping1.get_vertices(cell1).begin(), 4, vertices1.begin()); + + 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 + Delaunay tria_quad(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()) + { + 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)) + { + Delaunay tria_inter(vps->begin(), vps->end()); + for (auto it = tria_inter.finite_faces_begin(); + it != tria_inter.finite_faces_end(); + ++it) + { + if (CGAL::to_double( + tria_inter.triangle(it).squared_area()) > + tol * tol) + { + std::array, 3> verts = { + {CGALWrappers::cgal_point_to_dealii_point<3>( + it->vertex(0)->point()), + CGALWrappers::cgal_point_to_dealii_point<3>( + it->vertex(1)->point()), + CGALWrappers::cgal_point_to_dealii_point<3>( + it->vertex(2)->point())}}; + + vertices.push_back(verts); + } + } + } + } + } + } + + return vertices; + +# else + + Assert( + false, + ExcMessage( + "This function requires a version of CGAL greater or equal than 5.1.5.")); + (void)cell0; + (void)cell1; + (void)mapping0; + (void)mapping1; + (void)tol; + return {}; +# endif + } + + + template <> + std::vector, 4>> + compute_intersection_of_cells<3, 3, 3>( + const typename Triangulation<3, 3>::cell_iterator &cell0, + const typename Triangulation<3, 3>::cell_iterator &cell1, + const Mapping<3, 3> & mapping0, + const Mapping<3, 3> & mapping1, + const double tol) + { + 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); + CGAL::Polygon_mesh_processing::triangulate_faces(surf0); + CGAL::Polygon_mesh_processing::triangulate_faces(surf1); + CGALWrappers::compute_boolean_operation( + surf0, surf1, CGALWrappers::BooleanOperation::compute_intersection, sm); + if (CGAL::Polygon_mesh_processing::volume(sm) > tol) + { + // Collect tetrahedrons + std::vector, 4>> vertices; + 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))}}); + } + return vertices; + } + else + { + return {}; + } + } + +} // namespace CGALWrappers + +DEAL_II_NAMESPACE_CLOSE + +#else + +DEAL_II_NAMESPACE_OPEN + +template +std::vector, dim1 + 1>> +compute_intersection_of_cells( + const typename Triangulation::cell_iterator &cell0, + const typename Triangulation::cell_iterator &cell1, + const Mapping & mapping0, + const Mapping & mapping1, + const double tol = 1e-9) +{ + (void)cell0; + (void)cell1; + (void)mapping0; + (void)mapping1; + (void)tol; + AssertThrow(false, ExcMessage("This function needs CGAL to be installed.")); +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/cgal/cgal_intersect_simplices_1d_2d.cc b/tests/cgal/cgal_intersect_simplices_1d_2d.cc new file mode 100644 index 0000000000..435bbcab08 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_1d_2d.cc @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- + +// 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. +// +// --------------------------------------------------------------------- + +// Compute intersection of a segment embedded in 2D and a deal.II cell in 2D, +// and return a vector of arrays where you can build Quadrature rules. Then +// check that the sum of weights give the correct area for each region. + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +using namespace CGALWrappers; +static QGaussSimplex<1> qgauss(1); // use a degree 1 QGaussSimplex<2> + +void +test_inside_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, -0.45, 0.45); + const double expected_area = 0.9; + const auto cell0 = tria0.begin_active(); + 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 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; +} + + + +void +test_intersection(Triangulation<2> &tria0, Triangulation<1, 2> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, .5, 1.45); + const double expected_area = .5; + 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); +} diff --git a/tests/cgal/cgal_intersect_simplices_1d_2d.output b/tests/cgal/cgal_intersect_simplices_1d_2d.output new file mode 100644 index 0000000000..6dd05e6537 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_1d_2d.output @@ -0,0 +1,5 @@ + +DEAL::Inside: +DEAL::OK +DEAL::Intersecting: +DEAL::OK diff --git a/tests/cgal/cgal_intersect_simplices_1d_3d.cc b/tests/cgal/cgal_intersect_simplices_1d_3d.cc new file mode 100644 index 0000000000..bfd00f9d7a --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_1d_3d.cc @@ -0,0 +1,92 @@ +// --------------------------------------------------------------------- + +// 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. +// +// --------------------------------------------------------------------- + +// Compute intersection of a line with a cube and return a Quadrature formula +// over it. + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + + +void +test_intersection_inside(Triangulation<3> &tria0, Triangulation<1, 3> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, -0.4, 0.2); + tria1.begin_active()->vertex(0) = Point<3>(-.4, 0., 0.); + tria1.begin_active()->vertex(1) = Point<3>(.5, 0., 0.2); + const double expected_measure = + (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>()); + + 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; +} + + + +void +test_intersection(Triangulation<3> &tria0, Triangulation<1, 3> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1); + auto cell1 = tria1.begin_active(); + cell1->vertex(0) = Point<3>(); + 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 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() +{ + initlog(); + Triangulation<3> tria0; + Triangulation<1, 3> tria1; + + test_intersection_inside(tria0, tria1); + tria0.clear(); + tria1.clear(); + test_intersection(tria0, tria1); +} diff --git a/tests/cgal/cgal_intersect_simplices_1d_3d.output b/tests/cgal/cgal_intersect_simplices_1d_3d.output new file mode 100644 index 0000000000..8b3b075900 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_1d_3d.output @@ -0,0 +1,3 @@ + +DEAL::OK +DEAL::OK diff --git a/tests/cgal/cgal_intersect_simplices_2d_2d.cc b/tests/cgal/cgal_intersect_simplices_2d_2d.cc new file mode 100644 index 0000000000..d2d0cbf36b --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_2d_2d.cc @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- + +// 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. +// +// --------------------------------------------------------------------- + +// Compute intersection of simplices in 2D, and return a vector of arrays where +// you can build Quadrature rules. Then check that the sum of weights give the +// correct area for each region. + +#include + +#include +#include +#include + +#include + +#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) +{ + 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 vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<2, 2, 2>(cell0, + cell1, + MappingQ1<2>(), + MappingQ1<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; +} + + + +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 vec_of_arrays = + CGALWrappers::compute_intersection_of_cells<2, 2, 2>(cell0, + cell1, + MappingQ1<2>(), + MappingQ1<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; + Triangulation<2> tria1; + + test_inside_intersection(tria0, tria1); + + tria0.clear(); + tria1.clear(); + + test_intersection(tria0, tria1); +} diff --git a/tests/cgal/cgal_intersect_simplices_2d_2d.output b/tests/cgal/cgal_intersect_simplices_2d_2d.output new file mode 100644 index 0000000000..8b3b075900 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_2d_2d.output @@ -0,0 +1,3 @@ + +DEAL::OK +DEAL::OK diff --git a/tests/cgal/cgal_intersect_simplices_2d_3d.cc b/tests/cgal/cgal_intersect_simplices_2d_3d.cc new file mode 100644 index 0000000000..38070da4a6 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_2d_3d.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- + +// 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. +// +// --------------------------------------------------------------------- + +// Compute intersection of a square with a cube and return a Quadrature rule +// over its intersection. + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +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); + deallog << "OK" << std::endl; +} + + + +void +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); + deallog << "OK" << std::endl; +} + + +// +int +main() +{ + initlog(); + Triangulation<3> tria0; + Triangulation<2, 3> tria1; + + test_intersection_inside(tria0, tria1); + tria0.clear(); + tria1.clear(); + test_intersection(tria0, tria1); +} diff --git a/tests/cgal/cgal_intersect_simplices_2d_3d.output b/tests/cgal/cgal_intersect_simplices_2d_3d.output new file mode 100644 index 0000000000..8b3b075900 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_2d_3d.output @@ -0,0 +1,3 @@ + +DEAL::OK +DEAL::OK diff --git a/tests/cgal/cgal_intersect_simplices_3d_3d.cc b/tests/cgal/cgal_intersect_simplices_3d_3d.cc new file mode 100644 index 0000000000..fd61a9f319 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_3d_3d.cc @@ -0,0 +1,61 @@ +// --------------------------------------------------------------------- + +// 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. +// +// --------------------------------------------------------------------- + +// Compute intersection of two 3D cells, and return a Quadrature rule over it. + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +static QGaussSimplex<3> qgauss(1); +void +test_intersection(Triangulation<3> &tria0, Triangulation<3> &tria1) +{ + GridGenerator::hyper_cube(tria0, -1., 1.); + GridGenerator::hyper_cube(tria1, .5, 1.5); + static constexpr double expected_measure = 0.5 * 0.5 * 0.5; + 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 = qgauss.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() +{ + initlog(); + Triangulation<3> tria0; + Triangulation<3> tria1; + + test_intersection(tria0, tria1); +} diff --git a/tests/cgal/cgal_intersect_simplices_3d_3d.output b/tests/cgal/cgal_intersect_simplices_3d_3d.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/cgal/cgal_intersect_simplices_3d_3d.output @@ -0,0 +1,2 @@ + +DEAL::OK