Merge pull request #13664 from luca-heltai/cgal-surface-to-tria

Cgal surface to deal.II triangulation

doc/news/changes/major/20220426FederHeltai

@@ -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)

include/deal.II/cgal/triangulation.h

@@ -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>
 
@@ -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
@@ -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());
@@ -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
 

source/grid/tria.cc

@@ -11641,35 +11641,35 @@ Triangulation<dim, spacedim>::create_triangulation(
 
 
   /*
-      When the triangulation is a manifold (dim < spacedim), the normal field
-      provided from the map class depends on the order of the vertices.
-      It may happen that this normal field is discontinuous.
-      The following code takes care that this is not the case by setting the
-      cell direction flag on those cell that produce the wrong orientation.
-
-      To determine if 2 neighbours have the same or opposite orientation
-      we use a table of truth.
-      Its entries are indexes by the local indices of the common face.
-      For example if two elements share a face, and this face is
-      face 0 for element 0 and face 1 for element 1, then
-      table(0,1) will tell whether the orientation are the same (true) or
-      opposite (false).
-
-      Even though there may be a combinatorial/graph theory argument to get
-      this table in any dimension, I tested by hand all the different possible
-      cases in 1D and 2D to generate the table.
+      When the triangulation is a manifold (dim < spacedim) and made of
+      quadrilaterals, the normal field provided from the map class depends on
+      the order of the vertices. It may happen that this normal field is
+      discontinuous. The following code takes care that this is not the case by
+      setting the cell direction flag on those cell that produce the wrong
+      orientation.
+
+      To determine if 2 neighbours have the same or opposite orientation we use
+      a table of truth. Its entries are indexes by the local indices of the
+      common face. For example if two elements share a face, and this face is
+      face 0 for element 0 and face 1 for element 1, then table(0,1) will tell
+      whether the orientation are the same (true) or opposite (false).
+
+      Even though there may be a combinatorial/graph theory argument to get this
+      table in any dimension, I tested by hand all the different possible cases
+      in 1D and 2D to generate the table.
 
       Assuming that a surface respects the standard orientation for 2d meshes,
       the tables of truth are symmetric and their true values are the following
-      1D curves:  (0,1)
-      2D surface: (0,1),(0,2),(1,3),(2,3)
+
+      - 1D curves:  (0,1)
+      - 2D surface: (0,1),(0,2),(1,3),(2,3)
 
       We store this data using an n_faces x n_faces full matrix, which is
      actually much bigger than the minimal data required, but it makes the code
      more readable.
 
     */
-  if (dim < spacedim)
+  if (dim < spacedim && all_reference_cells_are_hyper_cube())
     {
       Table<2, bool> correct(GeometryInfo<dim>::faces_per_cell,
                              GeometryInfo<dim>::faces_per_cell);

tests/cgal/cgal_surface_triangulation_01.cc

@@ -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());
+  }
+}