CGALSegment2,
CGALTriangle2,
std::vector<CGALPoint2>>>
- compute_intersection(const std::array<Point<2>, 3> &first_simplex,
- const std::array<Point<2>, 3> &second_simplex)
+ compute_intersection_triangle_triangle(
+ const ArrayView<const Point<2>> &triangle0,
+ const ArrayView<const Point<2>> &triangle1)
{
+ AssertDimension(triangle0.size(), 3);
+ AssertDimension(triangle0.size(), triangle1.size());
+
std::array<CGALPoint2, 3> pts0, pts1;
- std::transform(
- first_simplex.begin(),
- first_simplex.end(),
- pts0.begin(),
- [&](const Point<2> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
- });
- std::transform(
- second_simplex.begin(),
- second_simplex.end(),
- pts1.begin(),
- [&](const Point<2> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(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<CGALPoint2, 2>);
+
+ std::transform(triangle1.begin(),
+ triangle1.end(),
+ pts1.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+ 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<std_cxx17::variant<CGALPoint2, CGALSegment2>>
- compute_intersection(const std::array<Point<2>, 3> &first_simplex,
- const std::array<Point<2>, 2> &second_simplex)
+ compute_intersection_triangle_segment(
+ const ArrayView<const Point<2>> &triangle,
+ const ArrayView<const Point<2>> &segment)
{
+ AssertDimension(triangle.size(), 3);
+ AssertDimension(segment.size(), 2);
+
std::array<CGALPoint2, 3> pts0;
std::array<CGALPoint2, 2> pts1;
- std::transform(
- first_simplex.begin(),
- first_simplex.end(),
- pts0.begin(),
- [&](const Point<2> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
- });
- std::transform(
- second_simplex.begin(),
- second_simplex.end(),
- pts1.begin(),
- [&](const Point<2> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(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<CGALPoint2, 2>);
+
+ std::transform(segment.begin(),
+ segment.end(),
+ pts1.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+ 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<Polygon_with_holes_2>
- compute_intersection(const std::array<Point<2>, 4> &first_simplex,
- const std::array<Point<2>, 4> &second_simplex)
+ compute_intersection_rect_rect(const ArrayView<const Point<2>> &rectangle0,
+ const ArrayView<const Point<2>> &rectangle1)
{
+ AssertDimension(rectangle0.size(), 4);
+ AssertDimension(rectangle0.size(), rectangle1.size());
+
std::array<CGALPoint2, 4> pts0, pts1;
- std::transform(
- first_simplex.begin(),
- first_simplex.end(),
- pts0.begin(),
- [&](const Point<2> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
- });
- std::transform(
- second_simplex.begin(),
- second_simplex.end(),
- pts1.begin(),
- [&](const Point<2> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(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<CGALPoint2, 2>);
+
+ std::transform(rectangle1.begin(),
+ rectangle1.end(),
+ pts1.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+ const CGALPolygon first_poly{pts0.begin(), pts0.end()};
+ const CGALPolygon second_poly{pts1.begin(), pts1.end()};
std::vector<Polygon_with_holes_2> 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<std_cxx17::variant<CGALPoint3, CGALSegment3>>
- compute_intersection(const std::array<Point<3>, 2> &first_simplex,
- const std::array<Point<3>, 4> &second_simplex)
+ compute_intersection_tetra_segment(
+ const ArrayView<const Point<3>> &tetrahedron,
+ const ArrayView<const Point<3>> &segment)
{
# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
+
+ AssertDimension(tetrahedron.size(), 4);
+ AssertDimension(segment.size(), 2);
+
std::array<CGALPoint3, 4> pts0;
std::array<CGALPoint3, 2> pts1;
- std::transform(
- first_simplex.begin(),
- first_simplex.end(),
- pts0.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
- });
+ std::transform(tetrahedron.begin(),
+ tetrahedron.end(),
+ pts0.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
- std::transform(
- second_simplex.begin(),
- second_simplex.end(),
- pts1.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(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<CGALPoint3, 3>);
+
+ 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
}
CGALSegment3,
CGALTriangle3,
std::vector<CGALPoint3>>>
- compute_intersection(const std::array<Point<3>, 3> &first_simplex,
- const std::array<Point<3>, 4> &second_simplex)
+ compute_intersection_tetra_triangle(
+ const ArrayView<const Point<3>> &tetrahedron,
+ const ArrayView<const Point<3>> &triangle)
{
# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
+
+ AssertDimension(tetrahedron.size(), 4);
+ AssertDimension(triangle.size(), 3);
+
std::array<CGALPoint3, 4> pts0;
std::array<CGALPoint3, 3> pts1;
- std::transform(
- first_simplex.begin(),
- first_simplex.end(),
- pts0.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(p);
- });
- std::transform(
- second_simplex.begin(),
- second_simplex.end(),
- pts1.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3>(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<CGALPoint3, 3>);
+
+ std::transform(triangle.begin(),
+ triangle.end(),
+ pts1.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3, 3>);
+
+ 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<std::array<Point<2>, 3>>
+ compute_intersection_quad_quad(const ArrayView<const Point<2>> &quad0,
+ const ArrayView<const Point<2>> &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<std::array<Point<2>, 3>>
- compute_intersection_of_cells<2, 2, 2, 4, 4>(
- const std::array<Point<2>, 4> &vertices0,
- const std::array<Point<2>, 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<std::array<Point<2>, 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<std::array<Point<2>, 3>> collection;
-
- CDT cdt;
- cdt.insert_constraint(poly.vertices_begin(),
- poly.vertices_end(),
- true);
-
- internal::mark_domains(cdt);
- std::array<Point<2>, 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<Point<2>, 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<std::array<Point<2>, 2>>
- compute_intersection_of_cells<2, 1, 2, 4, 2>(
- const std::array<Point<2>, 4> &vertices0,
- const std::array<Point<2>, 2> &vertices1,
- const double tol)
- {
- std::array<CGALPoint2, 4> pts;
- std::transform(
- vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<2> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(p);
- });
-
- CGALPolygon poly(pts.begin(), pts.end());
-
- CGALSegment2 segm(
- CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(vertices1[0]),
- CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(vertices1[1]));
- CDT cdt;
- cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true);
- std::vector<std::array<Point<2>, 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<CGALSegment2>(&*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<std::array<Point<2>, 2>>
+ compute_intersection_quad_line(const ArrayView<const Point<2>> &quad,
+ const ArrayView<const Point<2>> &line,
+ const double tol)
+ {
+ AssertDimension(quad.size(), 4);
+ AssertDimension(line.size(), 2);
+
+ std::array<CGALPoint2, 4> pts;
+
+ std::transform(quad.begin(),
+ quad.end(),
+ pts.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint2, 2>);
+
+ CGALPolygon poly(pts.begin(), pts.end());
+
+ CGALSegment2 segm(
+ CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[0]),
+ CGALWrappers::dealii_point_to_cgal_point<CGALPoint2>(line[1]));
+ CDT cdt;
+ cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true);
+ std::vector<std::array<Point<2>, 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<CGALSegment2>(&*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<std::array<Point<3>, 2>>
- compute_intersection_of_cells<3, 1, 3, 8, 2>(
- const std::array<Point<3>, 8> &vertices0,
- const std::array<Point<3>, 2> &vertices1,
- const double tol)
- {
+ return vertices;
+ }
+
+ // specialization for hex \cap line
+ std::vector<std::array<Point<3>, 2>>
+ compute_intersection_hexa_line(const ArrayView<const Point<3>> &hexa,
+ const ArrayView<const Point<3>> &line,
+ const double tol)
+ {
# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
- std::array<CGALPoint3_exact, 8> pts;
- std::transform(
- vertices0.begin(), vertices0.end(), pts.begin(), [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
- });
-
- CGALSegment3_exact segm(
- CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(vertices1[0]),
- CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(vertices1[1]));
-
- // Subdivide the hex into tetrahedrons, and intersect each one of them with
- // the line
- std::vector<std::array<Point<3>, 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<CGALSegment3_exact>(&*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<CGALPoint3_exact, 8> pts;
+
+ std::transform(
+ hexa.begin(),
+ hexa.end(),
+ pts.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
+
+ CGALSegment3_exact cgal_segment(
+ CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[0]),
+ CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(line[1]));
+
+ // Subdivide the hex into tetrahedrons, and intersect each one of them
+ // with the line
+ std::vector<std::array<Point<3>, 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<CGALSegment3_exact>(&*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<std::array<Point<3>, 3>>
- compute_intersection_of_cells<3, 2, 3, 8, 4>(
- const std::array<Point<3>, 8> &vertices0,
- const std::array<Point<3>, 4> &vertices1,
- const double tol)
- {
+ std::vector<std::array<Point<3>, 3>>
+ compute_intersection_hexa_quad(const ArrayView<const Point<3>> &hexa,
+ const ArrayView<const Point<3>> &quad,
+ const double tol)
+ {
# if DEAL_II_CGAL_VERSION_GTE(5, 5, 0)
- std::array<CGALPoint3_exact, 8> pts_hex;
- std::array<CGALPoint3_exact, 4> pts_quad;
- std::transform(
- vertices0.begin(),
- vertices0.end(),
- pts_hex.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
- });
-
- std::transform(
- vertices1.begin(),
- vertices1.end(),
- pts_quad.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact>(p);
- });
-
- // Subdivide hex into tetrahedrons
- std::vector<std::array<Point<3>, 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<CGALTriangle3_exact>(&*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<CGALPoint3_exact> *vps =
- boost::get<std::vector<CGALPoint3_exact>>(&*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<Point<3>, 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<CGALPoint3_exact, 8> pts_hex;
+ std::array<CGALPoint3_exact, 4> pts_quad;
+
+ std::transform(
+ hexa.begin(),
+ hexa.end(),
+ pts_hex.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
+
+ std::transform(
+ quad.begin(),
+ quad.end(),
+ pts_quad.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_exact, 3>);
+
+ // Subdivide hex into tetrahedrons
+ std::vector<std::array<Point<3>, 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<CGALTriangle3_exact>(&*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<CGALPoint3_exact> *vps =
+ boost::get<std::vector<CGALPoint3_exact>>(
+ &*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<Point<3>, 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<std::array<Point<3>, 4>>
+ compute_intersection_hexa_hexa(const ArrayView<const Point<3>> &hexa0,
+ const ArrayView<const Point<3>> &hexa1,
+ const double tol)
+ {
+ AssertDimension(hexa0.size(), 8);
+ AssertDimension(hexa0.size(), hexa1.size());
+ std::array<CGALPoint3_inexact, 8> pts_hex0;
+ std::array<CGALPoint3_inexact, 8> pts_hex1;
+ std::transform(
+ hexa0.begin(),
+ hexa0.end(),
+ pts_hex0.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact, 3>);
+
+ std::transform(
+ hexa1.begin(),
+ hexa1.end(),
+ pts_hex1.begin(),
+ &CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact, 3>);
+
+ Surface_mesh surf0, surf1, sm;
+ // Subdivide hex into tetrahedrons
+ std::vector<std::array<Point<3>, 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<std::array<Point<3>, 4>>
- compute_intersection_of_cells<3, 3, 3, 8, 8>(
- const std::array<Point<3>, 8> &vertices0,
- const std::array<Point<3>, 8> &vertices1,
- const double tol)
+ template <int dim0, int dim1, int spacedim>
+ std::vector<std::array<Point<spacedim>, dim1 + 1>>
+ compute_intersection_of_cells(
+ const ArrayView<const Point<spacedim>> &vertices0,
+ const ArrayView<const Point<spacedim>> &vertices1,
+ const double tol)
{
- std::array<CGALPoint3_inexact, 8> pts_hex0;
- std::array<CGALPoint3_inexact, 8> pts_hex1;
- std::transform(
- vertices0.begin(),
- vertices0.end(),
- pts_hex0.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
- });
-
- std::transform(
- vertices1.begin(),
- vertices1.end(),
- pts_hex1.begin(),
- [&](const Point<3> &p) {
- return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
- });
-
-
- Surface_mesh surf0, surf1, sm;
- // Subdivide hex into tetrahedrons
- std::vector<std::array<Point<3>, 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 <int dim0, int dim1, int spacedim>
std::vector<std::array<Point<spacedim>, dim1 + 1>>
compute_intersection_of_cells(
const Mapping<dim1, spacedim> & mapping1,
const double tol)
{
- Assert(mapping0.get_vertices(cell0).size() == std::pow(2, dim0),
+ Assert(mapping0.get_vertices(cell0).size() ==
+ ReferenceCells::get_hypercube<dim0>().n_vertices(),
ExcNotImplemented());
- Assert(mapping1.get_vertices(cell1).size() == std::pow(2, dim1),
+ Assert(mapping1.get_vertices(cell1).size() ==
+ ReferenceCells::get_hypercube<dim1>().n_vertices(),
ExcNotImplemented());
- const auto vertices0 =
- CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim0)>(
- cell0, mapping0);
- const auto vertices1 =
- CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim1)>(
- cell1, mapping1);
-
- return compute_intersection_of_cells<dim0,
- dim1,
- spacedim,
- Utilities::pow(2, dim0),
- Utilities::pow(2, dim1)>(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<dim0, dim1, spacedim>(vertices0,
+ vertices1,
+ tol);
}
# include "intersections.inst"
projects / dealii.git / commitdiff