-# # ---------------------------------------------------------------------
-# #
-# # 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.
-# #
-# # ---------------------------------------------------------------------
+## ---------------------------------------------------------------------
+##
+## 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.
+##
+## ---------------------------------------------------------------------
#
# Try to find the CGAL libraries
DEAL_II_NAMESPACE_OPEN
namespace CGALWrappers
{
- template <int dim0,
- int dim1,
- int spacedim,
- int n_components0,
- int n_components1>
- std::vector<std::array<Point<spacedim>, dim1 + 1>>
- compute_intersection_of_cells(
- const std::array<Point<spacedim>, n_components0> &vertices0,
- const std::array<Point<spacedim>, n_components1> &vertices1,
- const double tol)
- {
- (void)vertices0;
- (void)vertices1;
- (void)tol;
- Assert(false, ExcMessage("No explicit template instantiation available"));
- return {};
- }
-
/**
- * Given two deal.II affine cells, compute the intersection and return a
+ * Given two deal.II affine cells, compute their 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 is introduced. If the
const Mapping<dim0, spacedim> & mapping0,
const Mapping<dim1, spacedim> & mapping1,
const double tol = 1e-9);
+
+
+ /**
+ * Same function as above, but working directly with vertices.
+ */
+ template <int dim0, int dim1, int spacedim, int n_vertices0, int n_vertices1>
+ std::vector<std::array<Point<spacedim>, dim1 + 1>>
+ compute_intersection_of_cells(
+ const std::array<Point<spacedim>, n_vertices0> &vertices0,
+ const std::array<Point<spacedim>, n_vertices1> &vertices1,
+ const double tol = 1e-9)
+ {
+ (void)vertices0;
+ (void)vertices1;
+ (void)tol;
+ Assert(false, ExcMessage("No explicit template instantiation available"));
+ return {};
+ }
} // namespace CGALWrappers
DEAL_II_NAMESPACE_CLOSE
using CGALPoint2 = K::Point_2;
using CGALPoint3 = K::Point_3;
using CGALPoint3_exact = K_exact::Point_3;
+ using CGALPoint3_inexact = K_inexact::Point_3;
using CGALSegment2 = K::Segment_2;
using Surface_mesh = CGAL::Surface_mesh<K_inexact::Point_3>;
using CGALSegment3 = K::Segment_3;
}
}
+
+
void
mark_domains(CDT &cdt)
{
return CGAL::intersection(triangle1, triangle2);
}
+
+
boost::optional<boost::variant<CGALPoint2, CGALSegment2>>
compute_intersection(const std::array<Point<2>, 3> &first_simplex,
const std::array<Point<2>, 2> &second_simplex)
return CGAL::intersection(segm, triangle);
}
+
+
// rectangle-rectangle
std::vector<Polygon_with_holes_2>
compute_intersection(const std::array<Point<2>, 4> &first_simplex,
return poly_list;
}
+
+
boost::optional<boost::variant<CGALPoint3, CGALSegment3>>
compute_intersection(const std::array<Point<3>, 2> &first_simplex,
const std::array<Point<3>, 4> &second_simplex)
}
} // namespace internal
+
+
// Specialization for quads
template <>
std::vector<std::array<Point<2>, 3>>
}
}
+
+
// Specialization for quad \cap line
template <>
std::vector<std::array<Point<2>, 2>>
# endif
}
+
+
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)
+
+ 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);
- // std::vector<std::array<Point<3>, 4>> vertices;
- // if (CGAL::Polygon_mesh_processing::volume(sm) > tol)
- // {
- // // 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))}});
- // }
- // return vertices;
- // }
- // else
- // {
- // return vertices;
- // }
-
- // TODO: implememnt 3d/3d cut
- AssertThrow(false, ExcMessage("Not yet implemented"));
- (void)vertices0;
- (void)vertices1;
- (void)tol;
- return {};
+ 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 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())
+ {
+ 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);
+ const bool test_intersection =
+ CGAL::PMP::corefine_and_compute_intersection(surf0, surf1, sm);
+ if (CGAL::Polygon_mesh_processing::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();
+ }
+ surf0.clear();
+ }
+ return vertices;
}
+
+
template <int dim0, int dim1, int spacedim>
std::vector<std::array<Point<spacedim>, dim1 + 1>>
compute_intersection_of_cells(
const auto vertices0 =
CGALWrappers::get_vertices_in_cgal_order<int(std::pow(2, dim0))>(
cell0, mapping0);
-
const auto vertices1 =
CGALWrappers::get_vertices_in_cgal_order<int(std::pow(2, dim1))>(
cell1, mapping1);
projects / dealii.git / commitdiff