-
Notifications
You must be signed in to change notification settings - Fork 708
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Cgal triangulation #13661
Merged
Merged
Cgal triangulation #13661
Changes from all commits
Commits
Show all changes
3 commits
Select commit
Hold shift + click to select a range
File filter
Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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: | ||
<ul> | ||
<li>Conversion from deal.II to CGAL point types and viceversa</li> | ||
<li>Conversion from deal.II cell types to CGAL::Surface_mesh and viceversa</li> | ||
<li>Conversion from deal.II Triangulation to CGAL::Surface_mesh and viceversa</li> | ||
<li>Insertion of deal.II points in CGAL triangulation types</li> | ||
<li>Conversion from CGAL::Triangulation_2 to Triangulation<dim, 2></li> | ||
<li>Conversion from CGAL::Triangulation_3 to Triangulation<dim, 3></li> | ||
</ul> | ||
<br> | ||
(Marco Feder, Luca Heltai, 2022/04/26) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,279 @@ | ||
// --------------------------------------------------------------------- | ||
// | ||
// 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 <deal.II/base/config.h> | ||
|
||
#ifdef DEAL_II_WITH_CGAL | ||
|
||
# include <deal.II/grid/tria.h> | ||
|
||
# include <boost/hana.hpp> | ||
|
||
# include <CGAL/Surface_mesh.h> | ||
# include <CGAL/Triangulation_3.h> | ||
# include <deal.II/cgal/utilities.h> | ||
|
||
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 <int spacedim, typename CGALTriangulation> | ||
void | ||
add_points_to_cgal_triangulation(const std::vector<Point<spacedim>> &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[in] cgal_triangulation The input CGAL triangulation. | ||
* @param[out] dealii_triangulation The output deal.II triangulation. | ||
*/ | ||
template <typename CGALTriangulation, int dim, int spacedim> | ||
void | ||
cgal_triangulation_to_dealii_triangulation( | ||
const CGALTriangulation & cgal_triangulation, | ||
Triangulation<dim, spacedim> &dealii_triangulation); | ||
|
||
|
||
|
||
# ifndef DOXYGEN | ||
// Template implementation | ||
template <int spacedim, typename CGALTriangulation> | ||
void | ||
add_points_to_cgal_triangulation(const std::vector<Point<spacedim>> &points, | ||
CGALTriangulation &triangulation) | ||
{ | ||
Assert(triangulation.is_valid(), | ||
ExcMessage( | ||
"The triangulation you pass to this function should be a valid " | ||
"CGAL triangulation.")); | ||
|
||
std::vector<typename CGALTriangulation::Point> cgal_points(points.size()); | ||
std::transform(points.begin(), | ||
points.end(), | ||
cgal_points.begin(), | ||
[](const auto &p) { | ||
return CGALWrappers::dealii_point_to_cgal_point< | ||
typename CGALTriangulation::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 <typename CGALTriangulation, int dim, int spacedim> | ||
void | ||
cgal_triangulation_to_dealii_triangulation( | ||
const CGALTriangulation & cgal_triangulation, | ||
Triangulation<dim, spacedim> &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.")); | ||
|
||
// deal.II storage data structures | ||
std::vector<Point<spacedim>> vertices( | ||
cgal_triangulation.number_of_vertices()); | ||
std::vector<CellData<dim>> cells; | ||
SubCellData subcell_data; | ||
|
||
// CGAL storage data structures | ||
std::map<typename CGALTriangulation::Vertex_handle, unsigned int> | ||
vertex_map; | ||
{ | ||
unsigned int i = 0; | ||
for (const auto &v : cgal_triangulation.finite_vertex_handles()) | ||
{ | ||
vertices[i] = | ||
CGALWrappers::cgal_point_to_dealii_point<spacedim>(v->point()); | ||
vertex_map[v] = i++; | ||
} | ||
} | ||
|
||
const auto has_faces = boost::hana::is_valid( | ||
[](auto &&obj) -> decltype(obj.finite_face_handles()) {}); | ||
|
||
const 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<dim> cell(ReferenceCells::Triangle.n_vertices()); | ||
for (unsigned int i = 0; | ||
i < ReferenceCells::Triangle.n_vertices(); | ||
++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<dim> cell(ReferenceCells::Line.n_vertices()); | ||
unsigned int id = 0; | ||
// Since an edge is identified by a face (a triangle) and the | ||
// index of the opposite vertex to the edge, we can use this logic | ||
// to infer the indices of the vertices of the edge: loop over all | ||
// vertices, and keep only those that are not the opposite vertex | ||
// of the edge. | ||
for (unsigned int j = 0; | ||
j < ReferenceCells::Triangle.n_vertices(); | ||
++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<dim> cell(ReferenceCells::Tetrahedron.n_vertices()); | ||
for (unsigned int i = 0; | ||
i < ReferenceCells::Tetrahedron.n_vertices(); | ||
++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<dim> cell(ReferenceCells::Triangle.n_vertices()); | ||
unsigned int id = 0; | ||
// Since a face is identified by a cell (a tetrahedron) and the | ||
// index of the opposite vertex to the face, we can use this logic | ||
// to infer the indices of the vertices of the face: loop over all | ||
// vertices, and keep only those that are not the opposite vertex | ||
// of the face. | ||
for (unsigned int j = 0; | ||
j < ReferenceCells::Tetrahedron.n_vertices(); | ||
++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; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This is also a C++17 feature, isn't it? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes. The wrappers are only enabled if c++17 is enabled. |
||
CellData<dim> cell(ReferenceCells::Line.n_vertices()); | ||
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,56 @@ | ||
// --------------------------------------------------------------------- | ||
// | ||
// 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. | ||
// | ||
// --------------------------------------------------------------------- | ||
|
||
// Convert a vector of deal.II points to a cgal Triangulation | ||
|
||
#include <deal.II/base/config.h> | ||
|
||
#include <deal.II/grid/grid_generator.h> | ||
#include <deal.II/grid/tria.h> | ||
|
||
#include <CGAL/IO/io.h> | ||
#include <CGAL/Triangulation_3.h> | ||
#include <deal.II/cgal/triangulation.h> | ||
|
||
#include "../tests.h" | ||
|
||
using namespace CGALWrappers; | ||
|
||
using K = CGAL::Simple_cartesian<double>; | ||
using CGALTriangulation = CGAL::Triangulation_3<K>; | ||
|
||
template <int dim, int spacedim> | ||
void | ||
test() | ||
{ | ||
Triangulation<dim, spacedim> tria; | ||
GridGenerator::hyper_cube(tria); | ||
CGALTriangulation tr; | ||
add_points_to_cgal_triangulation(tria.get_vertices(), tr); | ||
deallog << "dim " << dim << ", spacedim " << spacedim << std::endl | ||
<< tr << std::endl; | ||
} | ||
|
||
int | ||
main() | ||
{ | ||
initlog(); | ||
test<1, 1>(); | ||
test<1, 2>(); | ||
test<1, 3>(); | ||
test<2, 2>(); | ||
test<2, 3>(); | ||
test<3, 3>(); | ||
} |
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We cannot use
if constexp
yet. That is only available once we switch to C++17. You will need to do template specialization.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
CGALWrappers are enabled only if c++17 is enabled.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK. I see.