Skip to content

Commit

Permalink
Conversion from CGAL surface to deal.II tria
Browse files Browse the repository at this point in the history
  • Loading branch information
@fdrmrc @luca-heltai
fdrmrc authored and luca-heltai committed May 8, 2022
1 parent ce3c581 commit 65984f2
Show file tree
Hide file tree
Showing 8 changed files with 662 additions and 1 deletion.
1 change: 1 addition & 0 deletions doc/news/changes/major/20220426FederHeltai
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ are supported:
<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>
<li>Conversion from CGAL::Surface_mesh or CGAL::Polyhedron_3 to Triangulation<2, 3></li>
</ul>
<br>
(Marco Feder, Luca Heltai, 2022/04/26)
146 changes: 145 additions & 1 deletion include/deal.II/cgal/triangulation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,9 @@

# include <boost/hana.hpp>

# include <CGAL/Polyhedron_3.h>
# include <CGAL/Surface_mesh.h>
# include <CGAL/Triangulation_2.h>
# include <CGAL/Triangulation_3.h>
# include <deal.II/cgal/utilities.h>

Expand Down Expand Up @@ -104,6 +106,32 @@ namespace CGALWrappers
const CGALTriangulation & cgal_triangulation,
Triangulation<dim, spacedim> &dealii_triangulation);

/**
* Construct a deal.II surface triangulation starting from a
* CGAL::Surface_mesh or a CGAL::Polyhedron_3.
*
* These types can both represent a polyhedral surface made of general
* polygons. deal.II only supports triangle or quadrilateral triangulations,
* and this function will throw an exception if the input surface mesh
* contains polygonal faces with more than 4 vertices per face.
*
* CGAL orients the faces of a polyhedral surface in a counter clock-wise
* w.r.t. to the material side of the surface, i.e., the normal is directed
* towards the outside of the surface if the surface is closed and represents
* a bounded domain. The opposite orientation is also perfectly valid, and
* would represent an unbounded domain complementary to the region
* represented by the surface. This function does not change the orientation
* of the surface mesh and will produce a triangulation with the same
* orientation of the input mesh.
*
* @param[in] cgal_mesh The input CGAL surface mesh.
* @param[out] triangulation The output deal.II triangulation.
*/
template <typename CGAL_MeshType>
void
cgal_surface_mesh_to_dealii_triangulation(const CGAL_MeshType &cgal_mesh,
Triangulation<2, 3> &triangulation);



# ifndef DOXYGEN
Expand Down Expand Up @@ -195,7 +223,7 @@ namespace CGALWrappers
// 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
// An edge is idendified 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());
Expand Down Expand Up @@ -270,6 +298,122 @@ namespace CGALWrappers
}
dealii_triangulation.create_triangulation(vertices, cells, subcell_data);
}



template <typename CGAL_MeshType>
void
cgal_surface_mesh_to_dealii_triangulation(const CGAL_MeshType &cgal_mesh,
Triangulation<2, 3> &triangulation)
{
Assert(triangulation.n_cells() == 0,
ExcMessage(
"Triangulation must be empty upon calling this function."));

const auto is_surface_mesh =
boost::hana::is_valid([](auto &&obj) -> decltype(obj.faces()) {});

const auto is_polyhedral =
boost::hana::is_valid([](auto &&obj) -> decltype(obj.facets_begin()) {});

// Collect Vertices and cells
std::vector<dealii::Point<3>> vertices;
std::vector<CellData<2>> cells;
SubCellData subcell_data;

// Different loops for Polyhedron or Surface_mesh types
if constexpr (decltype(is_surface_mesh(cgal_mesh)){})
{
AssertThrow(cgal_mesh.num_vertices() > 0,
ExcMessage("CGAL surface mesh is empty."));
vertices.reserve(cgal_mesh.num_vertices());
std::map<typename CGAL_MeshType::Vertex_index, unsigned int> vertex_map;
{
unsigned int i = 0;
for (const auto &v : cgal_mesh.vertices())
{
vertices.emplace_back(CGALWrappers::cgal_point_to_dealii_point<3>(
cgal_mesh.point(v)));
vertex_map[v] = i++;
}
}

// Collect CellData
for (const auto &face : cgal_mesh.faces())
{
const auto face_vertices =
CGAL::vertices_around_face(cgal_mesh.halfedge(face), cgal_mesh);

AssertThrow(face_vertices.size() == 3 || face_vertices.size() == 4,
ExcMessage("Only triangle or quadrilateral surface "
"meshes are supported in deal.II"));

CellData<2> c(face_vertices.size());
auto it_vertex = c.vertices.begin();
for (const auto &v : face_vertices)
*(it_vertex++) = vertex_map[v];

// Make sure the numberfing is consistent with the one in deal.II
if (face_vertices.size() == 4)
std::swap(c.vertices[3], c.vertices[2]);

cells.emplace_back(c);
}
}
else if constexpr (decltype(is_polyhedral(cgal_mesh)){})
{
AssertThrow(cgal_mesh.size_of_vertices() > 0,
ExcMessage("CGAL surface mesh is empty."));
vertices.reserve(cgal_mesh.size_of_vertices());
std::map<decltype(cgal_mesh.vertices_begin()), unsigned int> vertex_map;
{
unsigned int i = 0;
for (auto it = cgal_mesh.vertices_begin();
it != cgal_mesh.vertices_end();
++it)
{
vertices.emplace_back(
CGALWrappers::cgal_point_to_dealii_point<3>(it->point()));
vertex_map[it] = i++;
}
}

// Loop over faces of Polyhedron, fill CellData
for (auto face = cgal_mesh.facets_begin();
face != cgal_mesh.facets_end();
++face)
{
auto j = face->facet_begin();
const unsigned int vertices_per_face = CGAL::circulator_size(j);
AssertThrow(vertices_per_face == 3 || vertices_per_face == 4,
ExcMessage("Only triangle or quadrilateral surface "
"meshes are supported in deal.II. You "
"tried to read a mesh where a face has " +
std::to_string(vertices_per_face) +
" vertices per face."));

CellData<2> c(vertices_per_face);
auto it = c.vertices.begin();
for (unsigned int i = 0; i < vertices_per_face; ++i)
{
*(it++) = vertex_map[j->vertex()];
++j;
}

if (vertices_per_face == 4)
std::swap(c.vertices[3], c.vertices[2]);

cells.emplace_back(c);
}
}
else
{
AssertThrow(false,
ExcInternalError(
"Unsupported CGAL surface triangulation type."));
}
triangulation.create_triangulation(vertices, cells, subcell_data);
}
# endif
} // namespace CGALWrappers

Expand Down
63 changes: 63 additions & 0 deletions tests/cgal/cgal_surface_triangulation_01.cc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
// ---------------------------------------------------------------------
//
// 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 closed CGAL::Surface_mesh or a closed CGAL::Polyhedron_3 to a
// deal.II triangulation. The input mesh is a quad torus.

#include <deal.II/base/config.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>

#include <CGAL/IO/io.h>
#include <deal.II/cgal/triangulation.h>

#include <fstream>

#include "../tests.h"

using namespace CGALWrappers;
using Kernel = CGAL::Simple_cartesian<double>;
using CGALPoint = CGAL::Point_3<Kernel>;
using Mesh = CGAL::Surface_mesh<CGALPoint>;
using Polyhedron = CGAL::Polyhedron_3<Kernel>;

int
main()
{
initlog();
Triangulation<2, 3> tria;
GridOut go;
{
deallog << "Surface mesh" << std::endl;
Mesh sm;
std::ifstream input_sm(SOURCE_DIR "/input_grids/torus_quad.off");
input_sm >> sm;
cgal_surface_mesh_to_dealii_triangulation(sm, tria);
go.write_msh(tria, deallog.get_file_stream());
}
tria.clear();

// Now do the same with a Polyhedron
deallog << "Polyhedron test" << std::endl;
{
Polyhedron P;
std::ifstream input_p(SOURCE_DIR "/input_grids/torus_quad.off");
input_p >> P;
cgal_surface_mesh_to_dealii_triangulation(P, tria);
go.write_msh(tria, deallog.get_file_stream());
}
}
115 changes: 115 additions & 0 deletions tests/cgal/cgal_surface_triangulation_01.output
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@

DEAL::Surface mesh
$NOD
25
1 0.107677 -1.54743e-18 -0.331395
2 -0.281902 -1.54743e-18 -0.204814
3 -0.281902 -1.54743e-18 0.204814
4 0.107676 -1.54743e-18 0.331395
5 0.348450 -1.54743e-18 8.92032e-08
6 0.182080 0.331395 -0.560383
7 -0.476691 0.331395 -0.346336
8 -0.476691 0.331395 0.346336
9 0.182080 0.331395 0.560384
10 0.589223 0.331395 8.92032e-08
11 0.302466 0.204814 -0.930895
12 -0.791866 0.204814 -0.575325
13 -0.791866 0.204814 0.575324
14 0.302466 0.204814 0.930895
15 0.978801 0.204814 8.92032e-08
16 0.302466 -0.204814 -0.930895
17 -0.791866 -0.204814 -0.575325
18 -0.791866 -0.204814 0.575324
19 0.302466 -0.204814 0.930895
20 0.978801 -0.204814 8.92032e-08
21 0.182080 -0.331395 -0.560383
22 -0.476691 -0.331395 -0.346336
23 -0.476691 -0.331395 0.346336
24 0.182080 -0.331395 0.560384
25 0.589222 -0.331395 8.92032e-08
$ENDNOD
$ELM
25
1 3 0 0 4 6 7 2 1
2 3 0 0 4 7 8 3 2
3 3 0 0 4 8 9 4 3
4 3 0 0 4 9 10 5 4
5 3 0 0 4 10 6 1 5
6 3 0 0 4 11 12 7 6
7 3 0 0 4 12 13 8 7
8 3 0 0 4 13 14 9 8
9 3 0 0 4 14 15 10 9
10 3 0 0 4 15 11 6 10
11 3 0 0 4 16 17 12 11
12 3 0 0 4 17 18 13 12
13 3 0 0 4 18 19 14 13
14 3 0 0 4 19 20 15 14
15 3 0 0 4 20 16 11 15
16 3 0 0 4 21 22 17 16
17 3 0 0 4 22 23 18 17
18 3 0 0 4 23 24 19 18
19 3 0 0 4 24 25 20 19
20 3 0 0 4 25 21 16 20
21 3 0 0 4 1 2 22 21
22 3 0 0 4 2 3 23 22
23 3 0 0 4 3 4 24 23
24 3 0 0 4 4 5 25 24
25 3 0 0 4 5 1 21 25
$ENDELM
DEAL::Polyhedron test
$NOD
25
1 0.107677 -1.54743e-18 -0.331395
2 -0.281902 -1.54743e-18 -0.204814
3 -0.281902 -1.54743e-18 0.204814
4 0.107676 -1.54743e-18 0.331395
5 0.348450 -1.54743e-18 8.92032e-08
6 0.182080 0.331395 -0.560383
7 -0.476691 0.331395 -0.346336
8 -0.476691 0.331395 0.346336
9 0.182080 0.331395 0.560384
10 0.589223 0.331395 8.92032e-08
11 0.302466 0.204814 -0.930895
12 -0.791866 0.204814 -0.575325
13 -0.791866 0.204814 0.575324
14 0.302466 0.204814 0.930895
15 0.978801 0.204814 8.92032e-08
16 0.302466 -0.204814 -0.930895
17 -0.791866 -0.204814 -0.575325
18 -0.791866 -0.204814 0.575324
19 0.302466 -0.204814 0.930895
20 0.978801 -0.204814 8.92032e-08
21 0.182080 -0.331395 -0.560383
22 -0.476691 -0.331395 -0.346336
23 -0.476691 -0.331395 0.346336
24 0.182080 -0.331395 0.560384
25 0.589222 -0.331395 8.92032e-08
$ENDNOD
$ELM
25
1 3 0 0 4 6 7 2 1
2 3 0 0 4 7 8 3 2
3 3 0 0 4 8 9 4 3
4 3 0 0 4 9 10 5 4
5 3 0 0 4 10 6 1 5
6 3 0 0 4 11 12 7 6
7 3 0 0 4 12 13 8 7
8 3 0 0 4 13 14 9 8
9 3 0 0 4 14 15 10 9
10 3 0 0 4 15 11 6 10
11 3 0 0 4 16 17 12 11
12 3 0 0 4 17 18 13 12
13 3 0 0 4 18 19 14 13
14 3 0 0 4 19 20 15 14
15 3 0 0 4 20 16 11 15
16 3 0 0 4 21 22 17 16
17 3 0 0 4 22 23 18 17
18 3 0 0 4 23 24 19 18
19 3 0 0 4 24 25 20 19
20 3 0 0 4 25 21 16 20
21 3 0 0 4 1 2 22 21
22 3 0 0 4 2 3 23 22
23 3 0 0 4 3 4 24 23
24 3 0 0 4 4 5 25 24
25 3 0 0 4 5 1 21 25
$ENDELM

0 comments on commit 65984f2

Please sign in to comment.