From 9814fb92619852b7b577e7adbda988469e515450 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 16 May 2022 19:59:53 +0300 Subject: [PATCH] CGAL C3T3 to deal.II Triangulation. --- include/deal.II/cgal/triangulation.h | 153 +++++++++++++++++-- tests/cgal/cgal_triangulation_06.cc | 68 +++++++++ tests/cgal/cgal_triangulation_06.output | 193 ++++++++++++++++++++++++ 3 files changed, 403 insertions(+), 11 deletions(-) create mode 100644 tests/cgal/cgal_triangulation_06.cc create mode 100644 tests/cgal/cgal_triangulation_06.output diff --git a/include/deal.II/cgal/triangulation.h b/include/deal.II/cgal/triangulation.h index ffce229d39..53beb4dc74 100644 --- a/include/deal.II/cgal/triangulation.h +++ b/include/deal.II/cgal/triangulation.h @@ -106,6 +106,40 @@ namespace CGALWrappers const CGALTriangulation & cgal_triangulation, Triangulation &dealii_triangulation); + /** + * Specialization of the above function for + * CGAL::Mesh_complex_3_in_triangulation_3 types. + * + * CGAL::Mesh_complex_3_in_triangulation_3 is the class that implements the + * concept of a MeshComplexWithFeatures_3InTriangulation_3 in CGAL. + * + * https://doc.cgal.org/latest/Mesh_3/classMeshComplexWithFeatures__3InTriangulation__3.html + * + * This function translates the information contained in the input mesh + * complex (i.e., a collection of cells, surface patches, and curves + * embedded in a three dimensional simplicial complex) to a deal.II + * Triangulation object. + * + * This function ignores the information about the corners contained in the + * grids, but honors the information about curves, surface patches, and domain + * indices, which are all translated to manifold ids in the output + * Triangulation object. + * + * @param[in] cgal_triangulation An input Mesh_complex_3_in_triangulation_3 + * object + * @param[out] dealii_triangulation The output deal.II Triangulation object + */ + template + void + cgal_triangulation_to_dealii_triangulation( + const CGAL::Mesh_complex_3_in_triangulation_3 + & cgal_triangulation, + Triangulation<3> &dealii_triangulation); + /** * Construct a deal.II surface triangulation starting from a * CGAL::Surface_mesh or a CGAL::Polyhedron_3. @@ -164,6 +198,102 @@ namespace CGALWrappers + template + void + cgal_triangulation_to_dealii_triangulation( + const CGAL::Mesh_complex_3_in_triangulation_3 + & cgal_triangulation, + Triangulation<3> &dealii_triangulation) + { + using C3T3 = CGAL::Mesh_complex_3_in_triangulation_3; + + // Extract all vertices first + std::vector> dealii_vertices; + std::map + cgal_to_dealii_vertex_map; + + std::size_t inum = 0; + for (auto vit = cgal_triangulation.triangulation().finite_vertices_begin(); + vit != cgal_triangulation.triangulation().finite_vertices_end(); + ++vit) + { + if (vit->in_dimension() <= -1) + continue; + cgal_to_dealii_vertex_map[vit] = inum++; + dealii_vertices.emplace_back( + CGALWrappers::cgal_point_to_dealii_point<3>(vit->point())); + } + + // Now build cell connectivity + std::vector> cells; + for (auto cgal_cell = cgal_triangulation.cells_in_complex_begin(); + cgal_cell != cgal_triangulation.cells_in_complex_end(); + ++cgal_cell) + { + CellData<3> cell(4); + for (unsigned int i = 0; i < 4; ++i) + cell.vertices[i] = cgal_to_dealii_vertex_map[cgal_cell->vertex(i)]; + cell.manifold_id = cgal_triangulation.subdomain_index(cgal_cell); + cells.push_back(cell); + } + + // Do the same with surface patches, if possible + SubCellData subcell_data; + if constexpr (std::is_integral::value) + { + for (auto face = cgal_triangulation.facets_in_complex_begin(); + face != cgal_triangulation.facets_in_complex_end(); + ++face) + { + const auto &[cgal_cell, cgal_vertex_face_index] = *face; + CellData<2> dealii_face(3); + // A face is identified by a cell and by the index within the cell + // of the opposite vertex. Loop over vertices, and retain only those + // that belong to this face. + unsigned int j = 0; + for (unsigned int i = 0; i < 4; ++i) + if (i != cgal_vertex_face_index) + dealii_face.vertices[j++] = + cgal_to_dealii_vertex_map[cgal_cell->vertex(i)]; + dealii_face.manifold_id = + cgal_triangulation.surface_patch_index(cgal_cell, + cgal_vertex_face_index); + subcell_data.boundary_quads.emplace_back(dealii_face); + } + } + // and curves + if constexpr (std::is_integral::value) + { + for (auto edge = cgal_triangulation.edges_in_complex_begin(); + edge != cgal_triangulation.edges_in_complex_end(); + ++edge) + { + const auto &[cgal_cell, v1, v2] = *edge; + CellData<1> dealii_edge(2); + dealii_edge.vertices[0] = + cgal_to_dealii_vertex_map[cgal_cell->vertex(v1)]; + dealii_edge.vertices[1] = + cgal_to_dealii_vertex_map[cgal_cell->vertex(v2)]; + dealii_edge.manifold_id = + cgal_triangulation.curve_index(cgal_cell->vertex(v1), + cgal_cell->vertex(v2)); + subcell_data.boundary_lines.emplace_back(dealii_edge); + } + } + + dealii_triangulation.create_triangulation(dealii_vertices, + cells, + subcell_data); + } + + + template void cgal_triangulation_to_dealii_triangulation( @@ -223,16 +353,17 @@ namespace CGALWrappers // This is a degenerate Triangulation_2, made of edges for (const auto &e : cgal_triangulation.finite_edges()) { - // An edge is idendified 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 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. + // 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) @@ -262,17 +393,17 @@ namespace CGALWrappers // 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 + // 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(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. + // 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) diff --git a/tests/cgal/cgal_triangulation_06.cc b/tests/cgal/cgal_triangulation_06.cc new file mode 100644 index 0000000000..aad82fc0c3 --- /dev/null +++ b/tests/cgal/cgal_triangulation_06.cc @@ -0,0 +1,68 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// Read a surface mesh, make a coarse CGAL triangulation out of it, and +// translate the result into a deal.II Triangulation. + +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include "../tests.h" + +using namespace CGALWrappers; + +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using CGALPoint = K::Point_3; + +using Mesh_domain = + CGAL::Polyhedral_mesh_domain_with_features_3>; +using Tr = typename CGAL:: + Mesh_triangulation_3::type; +using Mesh_criteria = CGAL::Mesh_criteria_3; +using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3; +int +main() +{ + initlog(); + // Build a deal.II triangulation + CGAL::Surface_mesh cgal_surface_mesh; + { + std::ifstream input_mesh(SOURCE_DIR "/input_grids/tetrahedron.off"); + input_mesh >> cgal_surface_mesh; + } + + C3t3 cgal_triangulation; + cgal_surface_mesh_to_cgal_coarse_triangulation(cgal_surface_mesh, + cgal_triangulation); + + Triangulation<3> tria; + cgal_triangulation_to_dealii_triangulation(cgal_triangulation, tria); + + { + GridOut go; + std::ofstream out_tria("tria.vtk"); + go.write_vtk(tria, out_tria); + } + cat_file("tria.vtk"); +} diff --git a/tests/cgal/cgal_triangulation_06.output b/tests/cgal/cgal_triangulation_06.output new file mode 100644 index 0000000000..7bd156145c --- /dev/null +++ b/tests/cgal/cgal_triangulation_06.output @@ -0,0 +1,193 @@ + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 34 double +0 0 0 +0 0 1 +0 1 0 +1 0 0 +0.25 0.75 0 +0.333333 0.666667 0 +0.75 0.25 0 +0 0.666667 0 +0 0.333333 0 +0 0.75 0.25 +0 0.333333 0.666667 +0 0.25 0.75 +0.666667 0 0 +0.333333 0 0 +0.75 0 0.25 +0.333333 0 0.666667 +0.25 0 0.75 +0 0 0.333333 +0 0 0.666667 +0.416667 0 0.583333 +0.5 0 0.5 +0.583333 0 0.416667 +0.666667 0 0.333333 +0 0.416667 0.583333 +0 0.5 0.5 +0 0.583333 0.416667 +0 0.666667 0.333333 +0.416667 0.583333 0 +0.5 0.5 0 +0.583333 0.416667 0 +0.666667 0.333333 0 +0.544546 0 0 +0 0 0.544546 +0 0.544546 0 + +CELLS 132 524 +4 25 24 8 27 +4 6 12 3 14 +4 0 17 13 8 +4 14 6 12 31 +4 11 16 18 32 +4 22 13 30 31 +4 16 18 1 11 +4 26 9 4 33 +4 26 25 8 5 +4 15 10 17 32 +4 17 15 19 23 +4 2 9 7 4 +4 29 21 13 30 +4 8 17 13 28 +4 13 21 22 30 +4 8 26 5 33 +4 4 9 7 33 +4 20 13 17 28 +4 11 10 16 32 +4 10 15 17 23 +4 22 6 14 31 +4 20 17 24 28 +4 5 25 8 27 +4 19 17 23 24 +4 17 19 20 24 +4 17 8 24 28 +4 15 16 10 32 +4 27 24 8 28 +4 28 20 13 29 +4 21 13 20 29 +4 6 22 30 31 +4 5 26 4 33 +3 17 0 8 +3 0 17 13 +3 16 1 11 +3 9 2 4 +3 2 9 7 +3 6 3 14 +3 5 4 33 +3 12 6 31 +3 6 12 3 +3 2 7 4 +3 4 7 33 +3 8 5 33 +3 5 8 27 +3 26 8 33 +3 27 8 28 +3 7 9 33 +3 17 10 32 +3 10 17 23 +3 11 10 16 +3 10 11 32 +3 3 12 14 +3 14 12 31 +3 0 13 8 +3 8 13 28 +3 28 13 29 +3 29 13 30 +3 30 13 31 +3 22 14 31 +3 15 10 23 +3 15 16 10 +3 16 15 32 +3 19 15 23 +3 18 16 32 +3 16 18 1 +3 17 8 24 +3 17 15 19 +3 15 17 32 +3 23 17 24 +3 1 18 11 +3 11 18 32 +3 17 19 20 +3 20 19 24 +3 19 23 24 +3 20 13 17 +3 21 20 29 +3 20 28 29 +3 21 13 20 +3 13 21 22 +3 22 21 30 +3 22 6 14 +3 6 22 30 +3 13 22 31 +3 20 24 28 +3 24 27 28 +3 25 5 27 +3 25 24 8 +3 24 25 27 +3 5 26 4 +3 26 9 4 +3 9 26 33 +3 26 25 8 +3 25 26 5 +3 21 29 30 +3 6 30 31 +2 0 17 +2 0 8 +2 13 0 +2 1 16 +2 1 11 +2 2 9 +2 2 4 +2 7 2 +2 3 6 +2 3 14 +2 4 5 +2 12 31 +2 12 3 +2 7 33 +2 8 33 +2 5 27 +2 27 28 +2 17 32 +2 10 23 +2 11 10 +2 28 29 +2 29 30 +2 13 31 +2 14 22 +2 15 16 +2 15 19 +2 18 32 +2 18 1 +2 23 24 +2 19 20 +2 20 21 +2 21 22 +2 30 6 +2 25 24 +2 26 9 +2 26 25 + +CELL_TYPES 132 +10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 +5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 +3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + +CELL_DATA 132 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 +2 3 4 4 2 4 1 1 1 1 1 1 1 2 1 2 2 2 4 2 3 3 1 1 1 1 1 3 4 4 3 4 3 3 2 3 3 2 2 2 3 4 4 3 4 4 3 3 4 4 4 3 4 4 4 2 4 4 4 2 2 4 4 1 +6 2 4 5 3 3 1 2 1 5 1 4 4 2 2 1 1 6 3 3 1 1 4 5 5 5 6 6 3 5 5 5 1 3 3 3 + -- 2.39.5