constexpr int spacedim =
CGALTriangulationType::Point::Ambient_dimension::value;
- QGaussSimplex<spacedim> quad(degree);
- std::vector<dealii::Point<spacedim>> pts;
- std::vector<double> wts;
- std::array<dealii::Point<spacedim>, spacedim + 1> vertices; // tets
+ std::vector<std::array<dealii::Point<spacedim>, spacedim + 1>>
+ vec_of_simplices; // tets
- const auto is_c3t3 = boost::hana::is_valid(
- [](auto &&obj) -> decltype(obj.cells_in_complex_begin()) {});
-
- const auto is_tria3 = boost::hana::is_valid(
- [](auto &&obj) -> decltype(obj.finite_cell_handles()) {});
-
- if constexpr (decltype(is_c3t3(tria)){})
- {
- for (auto iter = tria.cells_in_complex_begin();
- iter != tria.cells_in_complex_end();
- ++iter)
- {
- for (unsigned int i = 0; i < (spacedim + 1); ++i)
- {
- vertices[i] = cgal_point_to_dealii_point<spacedim>(
- iter->vertex(i)->point());
- }
-
- auto local_quad = quad.compute_affine_transformation(vertices);
- std::transform(local_quad.get_points().begin(),
- local_quad.get_points().end(),
- std::back_inserter(pts),
- [&pts](const auto &p) { return p; });
- std::transform(local_quad.get_weights().begin(),
- local_quad.get_weights().end(),
- std::back_inserter(wts),
- [&wts](const double w) { return w; });
- }
- }
- else if constexpr (decltype(is_tria3(tria)){})
+ std::array<dealii::Point<spacedim>, spacedim + 1> simplex;
+ for (const auto &f : tria.finite_cell_handles())
{
- for (const auto &f : tria.finite_cell_handles())
+ for (unsigned int i = 0; i < (spacedim + 1); ++i)
{
- for (unsigned int i = 0; i < (spacedim + 1); ++i)
- {
- vertices[i] =
- cgal_point_to_dealii_point<spacedim>(f->vertex(i)->point());
- }
-
- auto local_quad = quad.compute_affine_transformation(vertices);
- std::transform(local_quad.get_points().begin(),
- local_quad.get_points().end(),
- std::back_inserter(pts),
- [&pts](const auto &p) { return p; });
- std::transform(local_quad.get_weights().begin(),
- local_quad.get_weights().end(),
- std::back_inserter(wts),
- [&wts](const double w) { return w; });
+ simplex[i] =
+ cgal_point_to_dealii_point<spacedim>(f->vertex(i)->point());
}
+
+ vec_of_simplices.push_back(simplex);
}
- else
- {
- Assert(false, ExcMessage("Not a valid CGAL Triangulation."));
- }
- return Quadrature<spacedim>(pts, wts);
+
+ return QGaussSimplex<spacedim>(degree).mapped_quadrature(vec_of_simplices);
}
{
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using CGALPoint = CGAL::Point_3<K>;
+ using CGALTriangulation = CGAL::Triangulation_3<K>;
CGAL::Surface_mesh<CGALPoint> surface_1, surface_2, out_surface;
dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surface_1);
dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surface_2);
// They have to be triangle meshes
CGAL::Polygon_mesh_processing::triangulate_faces(surface_1);
CGAL::Polygon_mesh_processing::triangulate_faces(surface_2);
- Assert(CGAL::is_triangle_mesh(surface_1),
+ Assert(CGAL::is_triangle_mesh(surface_1) &&
+ CGAL::is_triangle_mesh(surface_2),
ExcMessage("The surface must be a triangle mesh."));
compute_boolean_operation(surface_1, surface_2, bool_op, out_surface);
- using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<
- K,
- CGAL::Surface_mesh<CGALPoint>>;
- using Tr = typename CGAL::Mesh_triangulation_3<Mesh_domain,
- CGAL::Default,
- ConcurrencyTag>::type;
- using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
- using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
- C3t3 tria;
- cgal_surface_mesh_to_cgal_triangulation(out_surface, tria);
+ CGALTriangulation tria;
+ tria.insert(out_surface.points().begin(), out_surface.points().end());
+
return compute_quadrature(tria, degree);
}
}
+template <int dim>
+template <int spacedim>
+Quadrature<spacedim>
+QSimplex<dim>::mapped_quadrature(
+ const std::vector<std::array<Point<spacedim>, dim + 1>> &simplices) const
+{
+ Assert(!(dim == 1 && spacedim == 1),
+ ExcMessage("This function is not supposed to work in 1D-1D case."));
+ Assert(dim <= spacedim,
+ ExcMessage("Invalid combination of dim and spacedim ."));
+
+ std::vector<Point<spacedim>> qp;
+ std::vector<double> ws;
+ for (const auto &simplex : simplices)
+ {
+ const auto rule = this->compute_affine_transformation(simplex);
+ std::transform(rule.get_points().begin(),
+ rule.get_points().end(),
+ std::back_inserter(qp),
+ [&](const Point<spacedim> &p) { return p; });
+ std::transform(rule.get_weights().begin(),
+ rule.get_weights().end(),
+ std::back_inserter(ws),
+ [&](const double w) { return w; });
+ }
+ return Quadrature<spacedim>(qp, ws);
+}
+
+
+
QTrianglePolar::QTrianglePolar(const Quadrature<1> &radial_quadrature,
const Quadrature<1> &angular_quadrature)
: QSimplex<2>(Quadrature<2>())
QSimplex<2>::compute_affine_transformation(
const std::array<Point<3>, 2 + 1> &vertices) const;
} // namespace dealii
+
+namespace dealii
+{
+ template Quadrature<2>
+ QSimplex<1>::mapped_quadrature(
+ const std::vector<std::array<Point<2>, 1 + 1>> &simplices) const;
+
+ template Quadrature<3>
+ QSimplex<1>::mapped_quadrature(
+ const std::vector<std::array<Point<3>, 1 + 1>> &simplices) const;
+
+ template Quadrature<2>
+ QSimplex<2>::mapped_quadrature(
+ const std::vector<std::array<Point<2>, 2 + 1>> &simplices) const;
+
+ template Quadrature<3>
+ QSimplex<2>::mapped_quadrature(
+ const std::vector<std::array<Point<3>, 2 + 1>> &simplices) const;
+
+ template Quadrature<3>
+ QSimplex<3>::mapped_quadrature(
+ const std::vector<std::array<Point<3>, 3 + 1>> &simplices) const;
+} // namespace dealii
// Compute volumes by splitting polyhedra into simplices.
-using K = CGAL::Exact_predicates_inexact_constructions_kernel;
-using CGALPoint = CGAL::Point_3<K>;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using CGALPoint = CGAL::Point_3<K>;
+using CGALTriangulation = CGAL::Triangulation_3<K>;
using namespace CGALWrappers;
-using Mesh_domain =
- CGAL::Polyhedral_mesh_domain_with_features_3<K,
- CGAL::Surface_mesh<CGALPoint>>;
-using Tr = typename CGAL::
- Mesh_triangulation_3<Mesh_domain, CGAL::Default, ConcurrencyTag>::type;
-using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
-using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr,
- Mesh_domain::Corner_index,
- Mesh_domain::Curve_index>;
void
SOURCE_DIR
"/input_grids/octahedron.off"};
CGAL::Surface_mesh<CGALPoint> sm;
- C3t3 tria;
+ CGALTriangulation tria;
constexpr int degree = 3;
for (const auto &name : fnames)
{
std::ifstream input(name);
input >> sm;
- cgal_surface_mesh_to_cgal_triangulation(sm, tria);
+ tria.insert(sm.points().begin(), sm.points().end());
auto b = compute_quadrature(tria, degree);
deallog << "Volume of poly with Quadrature: " << std::setprecision(12)
<< std::accumulate(b.get_weights().begin(),
projects / dealii.git / commitdiff