]> https://gitweb.dealii.org/ - dealii.git/commitdiff
CGAL C3T3 to deal.II Triangulation.
authorLuca Heltai <luca.heltai@sissa.it>
Mon, 16 May 2022 16:59:53 +0000 (19:59 +0300)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 17 May 2022 19:21:29 +0000 (22:21 +0300)
include/deal.II/cgal/triangulation.h
tests/cgal/cgal_triangulation_06.cc [new file with mode: 0644]
tests/cgal/cgal_triangulation_06.output [new file with mode: 0644]

index ffce229d39555b091150486ad85a3b549bfe2b58..53beb4dc742b13a5c6ca38005f350181090caa7c 100644 (file)
@@ -106,6 +106,40 @@ namespace CGALWrappers
     const CGALTriangulation &     cgal_triangulation,
     Triangulation<dim, spacedim> &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 <typename CGALTriangulationType,
+            typename CornerIndexType,
+            typename CurveIndexType>
+  void
+  cgal_triangulation_to_dealii_triangulation(
+    const CGAL::Mesh_complex_3_in_triangulation_3<CGALTriangulationType,
+                                                  CornerIndexType,
+                                                  CurveIndexType>
+      &               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 <typename CGALTriangulationType,
+            typename CornerIndexType,
+            typename CurveIndexType>
+  void
+  cgal_triangulation_to_dealii_triangulation(
+    const CGAL::Mesh_complex_3_in_triangulation_3<CGALTriangulationType,
+                                                  CornerIndexType,
+                                                  CurveIndexType>
+      &               cgal_triangulation,
+    Triangulation<3> &dealii_triangulation)
+  {
+    using C3T3 = CGAL::Mesh_complex_3_in_triangulation_3<CGALTriangulationType,
+                                                         CornerIndexType,
+                                                         CurveIndexType>;
+
+    // Extract all vertices first
+    std::vector<Point<3>> dealii_vertices;
+    std::map<typename C3T3::Vertex_handle, unsigned int>
+      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<CellData<3>> 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<typename C3T3::Surface_patch_index>::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<typename C3T3::Curve_index>::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 <typename CGALTriangulation, int dim, int spacedim>
   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<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.
+              // 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<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.
+              // 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 (file)
index 0000000..aad82fc
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/io.h>
+#include <deal.II/cgal/surface_mesh.h>
+#include <deal.II/cgal/triangulation.h>
+#include <deal.II/cgal/utilities.h>
+
+#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<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>;
+int
+main()
+{
+  initlog();
+  // Build a deal.II triangulation
+  CGAL::Surface_mesh<CGALPoint> 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 (file)
index 0000000..7bd1561
--- /dev/null
@@ -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 
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.