From 4c8162e4820edfabcccc5927836a7906158f33fa Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sat, 7 May 2022 21:25:49 +0300 Subject: [PATCH] Support for CGAL::Triangulation_2 and CGAL::Triangulation_3. --- doc/news/changes/major/20220426FederHeltai | 11 +- include/deal.II/cgal/triangulation.h | 259 +++++++++++++++++++++ 2 files changed, 269 insertions(+), 1 deletion(-) create mode 100644 include/deal.II/cgal/triangulation.h diff --git a/doc/news/changes/major/20220426FederHeltai b/doc/news/changes/major/20220426FederHeltai index e4d7d88dda..1680193d81 100644 --- a/doc/news/changes/major/20220426FederHeltai +++ b/doc/news/changes/major/20220426FederHeltai @@ -1,3 +1,12 @@ -New: Added support for the CGAL library (www.cgal.org). +New: Added support for the CGAL library (www.cgal.org). The following features +are supported: +
(Marco Feder, Luca Heltai, 2022/04/26) diff --git a/include/deal.II/cgal/triangulation.h b/include/deal.II/cgal/triangulation.h new file mode 100644 index 0000000000..4dffce2a7b --- /dev/null +++ b/include/deal.II/cgal/triangulation.h @@ -0,0 +1,259 @@ +// --------------------------------------------------------------------- +// +// 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_triangulation_h +#define dealii_cgal_triangulation_h + +#include + +#include + +#include + +#ifdef DEAL_II_WITH_CGAL +# include + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace CGALWrappers +{ + /** + * Build any CGAL triangulation from a list of deal.II points. + * + * The CGAL library offers several types of triangulations, including surface + * triangulations, cell complexes, Delaunay triangulations in two and three + * dimensions, and many others. For simple cases, all of these triangulations + * allow you to incrementally build them from a list of points. This function + * provides a convenient way to add a vector of deal.II points to any CGAL + * triangulation that admits the insertion of new points via iterators. + * + * The resulting triangulation will be the convex-hull of the input points. If + * you need to build more complicated triangulations, you can use the + * functions that build a triangulation from a polyhedral surface domain. + * + * More information on the available CGAL triangulation classes is available + * at https://doc.cgal.org/latest/Triangulation_2/index.html for two + * dimensional triangulations, and at + * https://doc.cgal.org/latest/Triangulation_3/index.html for three + * dimensional triangulations. + * + * Notice that CGAL distinguishes between a triangulation and a polygonal or + * polyhedral mesh. Generally speaking, a triangulation is made of simplices, + * whereas a polygonal or polyhedral mesh is made of general polygons or + * general polyhedrons. While CGAL implements the two concepts in a similar + * fashion using half-edge data structures, some optimizations are performed + * for triangulations, where only triangles or tetrahedra are used. + * + * @param[in] points The input points to build the triangulation from. + * @param[out] triangulation The output triangulation. + */ + template + void + add_points_to_cgal_triangulation(const std::vector> &points, + CGALTriangulation &triangulation); + + /** + * Convert any compatible CGAL triangulation type to a deal.II triangulation. + * + * The conversion is done by looping over the finite vertices and finite cells + * of the CGAL triangulation, and building with them a deal.II Triangulation + * object. + * + * CGAL considers triangulations to be a partition of the entire + * spacedim-dimensional space, where every face (including boundary faces) is + * shared between neighboring cells. This is obtained by adding to the list of + * vertices a special vertex (the vertex at infinity, representing a point on + * the spacedim-dimensional sphere with infinite radius), and adding a + * neighboring infinite cell to each boundary face. + * + * Valid conversions require that the capacity of the deal.II triangulation + * matches that of the input CGAL triangulation, i.e., if the input CGAL + * triangulation is a CGAL::Triangulation_2, then spacedim must be greater + * or equal than two. If the input CGAL triangulation is a + * CGAL::Triangulation_3, then spacedim must be greater or equal than three. + * + * Both CGAL::Triangulation_2 and CGAL::Triangulation_3 can store degenerate + * two dimensional triangulations (i.e., a one-dimensional triangulation), or + * a degenerate three dimensional triangulation (i.e., a two-dimensional + * surface triangulation). The dimension of the triangulation, as returned by + * `cgal_triangulation.dimension()`, must match the dimension `dim` of the + * input Triangulation `dealii_triangulation`. If this is not the case, an + * exception is thrown. + * + * @param cgal_triangulation + * @param dealii_triangulation + */ + template + void + cgal_triangulation_to_dealii_triangulation( + const CGALTriangulation & cgal_triangulation, + Triangulation &dealii_triangulation); + + + +# ifndef DOXYGEN + // Template implementation + template + void + add_points_to_cgal_triangulation(const std::vector> &points, + CGALTriangulation &triangulation) + { + Assert(triangulation.is_valid(), + ExcMessage( + "The triangulation you pass to this function should be a valid " + "CGAL triangulation.")); + using CGALPoint = typename CGALTriangulation::Point; + std::vector cgal_points(points.size()); + std::transform( + points.begin(), points.end(), cgal_points.begin(), [](const auto &p) { + return CGALWrappers::dealii_point_to_cgal_point(p); + }); + + triangulation.insert(cgal_points.begin(), cgal_points.end()); + Assert(triangulation.is_valid(), + ExcMessage( + "The Triangulation is no longer valid after inserting the points. " + "Bailing out.")); + } + + + + template + void + cgal_triangulation_to_dealii_triangulation( + const CGALTriangulation & cgal_triangulation, + Triangulation &dealii_triangulation) + { + AssertThrow(cgal_triangulation.dimension() == dim, + ExcMessage("The dimension of the input CGAL triangulation (" + + std::to_string(cgal_triangulation.dimension()) + + ") does not match the dimension of the output " + "deal.II triangulation (" + + std::to_string(dim) + ").")); + + Assert(dealii_triangulation.n_cells() == 0, + ExcMessage("The output triangulation object needs to be empty.")); + + const auto nv = cgal_triangulation.number_of_vertices(); + + // Deal.II storage data structures + std::vector> vertices(nv); + std::vector> cells; + SubCellData subcell_data; + + // CGAL storage data structures + std::map + vertex_map; + { + unsigned int i = 0; + for (auto v : cgal_triangulation.finite_vertex_handles()) + { + vertices[i] = + CGALWrappers::cgal_point_to_dealii_point(v->point()); + vertex_map[v] = i++; + } + } + + auto has_faces = boost::hana::is_valid( + [](auto &&obj) -> decltype(obj.finite_face_handles()) {}); + + auto has_cells = boost::hana::is_valid( + [](auto &&obj) -> decltype(obj.finite_cell_handles()) {}); + + // Different loops for Triangulation_2 and Triangulation_3 types. + if constexpr (decltype(has_faces(cgal_triangulation)){}) + { + // This is a non-degenerate Triangulation_2 + if (cgal_triangulation.dimension() == 2) + for (const auto f : cgal_triangulation.finite_face_handles()) + { + CellData cell(3); + for (unsigned int i = 0; i < 3; ++i) + cell.vertices[i] = vertex_map[f->vertex(i)]; + cells.push_back(cell); + } + else if (cgal_triangulation.dimension() == 1) + // This is a degenerate Triangulation_2, made of edges + for (const auto e : cgal_triangulation.finite_edges()) + { + // An edge is idenfied by a face and a vertex index in the face + const auto & f = e.first; + const auto & i = e.second; + CellData cell(2); + unsigned int id = 0; + for (unsigned int j = 0; j < 3; ++j) + if (j != i) + cell.vertices[id++] = vertex_map[f->vertex(j)]; + cells.push_back(cell); + } + else + { + Assert(false, ExcInternalError()); + } + } + else if constexpr (decltype(has_cells(cgal_triangulation)){}) + { + // This is a non-degenerate Triangulation_3 + if (cgal_triangulation.dimension() == 3) + for (const auto c : cgal_triangulation.finite_cell_handles()) + { + CellData cell(4); + for (unsigned int i = 0; i < 4; ++i) + cell.vertices[i] = vertex_map[c->vertex(i)]; + cells.push_back(cell); + } + else if (cgal_triangulation.dimension() == 2) + // This is a degenerate Triangulation_3, made of triangles + for (const auto facet : cgal_triangulation.finite_facets()) + { + // A facet is idenfied by a cell and the opposite vertex index in + // the face + const auto & c = facet.first; + const auto & i = facet.second; + CellData cell(3); + unsigned int id = 0; + for (unsigned int j = 0; j < 4; ++j) + if (j != i) + cell.vertices[id++] = vertex_map[c->vertex(j)]; + cells.push_back(cell); + } + else if (cgal_triangulation.dimension() == 1) + // This is a degenerate Triangulation_3, made of edges + for (const auto edge : cgal_triangulation.finite_edges()) + { + // An edge is idenfied by a cell and its two vertices + const auto &[c, i, j] = edge; + CellData cell(2); + cell.vertices[0] = vertex_map[c->vertex(i)]; + cell.vertices[1] = vertex_map[c->vertex(j)]; + cells.push_back(cell); + } + else + { + Assert(false, ExcInternalError()); + } + } + dealii_triangulation.create_triangulation(vertices, cells, subcell_data); + } +# endif +} // namespace CGALWrappers + +DEAL_II_NAMESPACE_CLOSE + +#endif +#endif -- 2.39.5