Conversion from CGAL surface to deal.II tria

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

@@ -18,13 +18,16 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>
 
 #include <deal.II/cgal/utilities.h>
 
 #ifdef DEAL_II_WITH_CGAL
 #  include <boost/hana.hpp>
 
+#  include <CGAL/Polyhedron_3.h>
+#  include <CGAL/Polyhedron_items_with_id_3.h>
 #  include <CGAL/Surface_mesh.h>
 #  include <CGAL/Triangulation_3.h>
 
@@ -100,6 +103,27 @@ namespace CGALWrappers
     const CGALTriangulation &     cgal_triangulation,
     Triangulation<dim, spacedim> &dealii_triangulation);
 
+  /**
+   * @brief Construct a deal.II triangulation starting from a CGAL_MeshType object,
+   * which can either be a CGAL::Surface_mesh or a CGAL::Polyhedron_3. These
+   * types can both represent a polyhedral surface. Therefore, the natural
+   * deal.II object to fill with this function is a Triangulation<2,3>.
+   *
+   * This is created in the same spirit of the function above, i.e. by looping
+   * over all the vertices and cells of the surface.
+   *
+   * @tparam CGAL_MeshType
+   * @tparam dim
+   * @tparam spacedim
+   * @param cgal_mesh
+   * @param triangulation
+   */
+  template <typename CGAL_MeshType, int dim, int spacedim>
+  void
+  cgal_surface_to_dealii_triangulation(
+    const CGAL_MeshType &         cgal_mesh,
+    Triangulation<dim, spacedim> &triangulation);
+
 
 
 #  ifndef DOXYGEN
@@ -246,6 +270,130 @@ namespace CGALWrappers
       }
     dealii_triangulation.create_triangulation(vertices, cells, subcell_data);
   }
+
+
+
+  template <typename CGAL_MeshType, int dim, int spacedim>
+  void
+  cgal_surface_to_dealii_triangulation(
+    const CGAL_MeshType &         cgal_mesh,
+    Triangulation<dim, spacedim> &triangulation)
+  {
+    Assert(triangulation.n_cells() == 0,
+           ExcMessage(
+             "Triangulation must be empty upon calling this function."));
+    Assert(dim == 2 && spacedim == 3,
+           ExcMessage("This function works only with dim==2 and spacedim==3."));
+
+
+    auto is_surface_mesh =
+      boost::hana::is_valid([](auto &&obj) -> decltype(obj.faces()) {});
+
+    auto is_polyhedral =
+      boost::hana::is_valid([](auto &&obj) -> decltype(obj.facets_begin()) {});
+
+    // Collect Vertices
+    std::size_t                          n_cgal_vertices;
+    std::vector<dealii::Point<spacedim>> vertices;
+    std::vector<CellData<dim>>           cells;
+    SubCellData                          subcell_data;
+    if constexpr (is_surface_mesh(cgal_mesh))
+      {
+        n_cgal_vertices = cgal_mesh.num_vertices();
+        vertices.reserve(n_cgal_vertices);
+        for (const auto &p : cgal_mesh.points())
+          {
+            vertices.emplace_back(
+              CGALWrappers::cgal_point_to_dealii_point<spacedim>(p));
+          }
+      }
+    else if constexpr (is_polyhedral(cgal_mesh))
+      {
+        n_cgal_vertices = cgal_mesh.size_of_vertices();
+        vertices.reserve(n_cgal_vertices);
+        for (auto it = cgal_mesh.points_begin(); it != cgal_mesh.points_end();
+             ++it)
+          {
+            vertices.emplace_back(
+              CGALWrappers::cgal_point_to_dealii_point<spacedim>(*it));
+          }
+      }
+
+
+    // Different loops for Polyhedron or Surface_mesh types
+    if constexpr (is_surface_mesh(cgal_mesh))
+      {
+        const unsigned int vertices_per_face =
+          CGAL::vertices_around_face(
+            cgal_mesh.halfedge(*(cgal_mesh.faces().begin())), cgal_mesh)
+            .size();
+
+        // Collect CellData
+        for (const auto &face : cgal_mesh.faces())
+          {
+            CellData<dim> c(vertices_per_face);
+            auto          it_vertex = c.vertices.begin();
+            for (const auto v :
+                 CGAL::vertices_around_face(cgal_mesh.halfedge(face),
+                                            cgal_mesh))
+              {
+                *(it_vertex++) = v;
+              }
+
+            if (vertices_per_face == 4)
+              std::swap(c.vertices[3], c.vertices[2]);
+
+            cells.emplace_back(c);
+          }
+        GridTools::delete_unused_vertices(vertices, cells, subcell_data);
+        GridTools::consistently_order_cells(cells);
+        triangulation.create_triangulation(vertices, cells, subcell_data);
+      }
+    else if constexpr (is_polyhedral(cgal_mesh))
+      {
+        const unsigned int vertices_per_face =
+          cgal_mesh.facets_begin()->facet_degree();
+
+        auto cgal_mesh_nc = const_cast<CGAL_MeshType *>(
+          &cgal_mesh); // Need to change vertex->id() for the Polyhedron
+        std::size_t i = 0;
+        for (auto vertex_it = cgal_mesh_nc->vertices_begin();
+             vertex_it != cgal_mesh_nc->vertices_end();
+             ++vertex_it)
+          {
+            vertex_it->id() = i++;
+          }
+        // Loop over faces of Polyhedron, fill CellData
+        for (auto face_it = cgal_mesh_nc->facets_begin();
+             face_it != cgal_mesh_nc->facets_end();
+             ++face_it)
+          {
+            CellData<dim> c(vertices_per_face);
+            auto          it   = c.vertices.begin();
+            auto          circ = face_it->facet_begin();
+            do
+              {
+                *(it++) = circ->vertex()->id();
+              }
+            while (++circ != face_it->facet_begin());
+
+            if (vertices_per_face == 4)
+              std::swap(c.vertices[3], c.vertices[2]);
+
+            cells.emplace_back(c);
+          }
+        if (vertices_per_face == 4)
+          {
+            GridTools::delete_unused_vertices(vertices, cells, subcell_data);
+            GridTools::consistently_order_cells(cells);
+            triangulation.create_triangulation(vertices, cells, subcell_data);
+          }
+      }
+    else
+      {
+        Assert(false, ExcInternalError());
+      }
+  }
 #  endif
 } // namespace CGALWrappers
 

tests/cgal/cgal_surface_triangulation_06.cc

@@ -0,0 +1,88 @@
+// ---------------------------------------------------------------------
+//
+// 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 CGAL::Surface_mesh or a CGAL::Polyhedron_3 to a deal.II
+// triangulation. CGAL surfaces are generated using an .off file describing
+// the shape of the number '8'.
+
+#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 CGALPoint  = CGAL::Point_3<CGAL::Simple_cartesian<double>>;
+using Mesh       = CGAL::Surface_mesh<CGALPoint>;
+using Polyhedron = CGAL::Polyhedron_3<CGAL::Simple_cartesian<double>,
+                                      CGAL::Polyhedron_items_with_id_3>;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  Mesh          sm;
+  std::ifstream input_sm("eight.off");
+  input_sm >> sm;
+
+  Triangulation<dim, spacedim> tria;
+  cgal_surface_to_dealii_triangulation(sm, tria);
+  tria.refine_global(1);
+  GridOut       go;
+  std::ofstream out("surface_tria.vtk");
+  go.write_vtk(tria, out);
+  go.write_vtk(tria, deallog.get_file_stream());
+
+  // Now do the same with a Polyhedron
+  deallog << "Polyhedron test" << std::endl;
+  Polyhedron    P;
+  std::ifstream input_p("eight.off");
+  input_p >> P;
+
+
+  tria.clear();
+  cgal_surface_to_dealii_triangulation(P, tria);
+  tria.refine_global(1);
+  std::ofstream output("poly_tria.vtk");
+  go.write_vtk(tria, output);
+  go.write_vtk(tria, deallog.get_file_stream());
+
+  // Show enumeration of each facet
+  for (auto face_it = P.facets_begin(); face_it != P.facets_end(); ++face_it)
+    {
+      auto circ = face_it->facet_begin();
+      deallog << "Number of vertices per facet: " << CGAL::circulator_size(circ)
+              << std::endl;
+      do
+        {
+          deallog << "Vertex id: " << circ->vertex()->id() << std::endl;
+        }
+      while (++circ != face_it->facet_begin());
+    }
+}
+
+int
+main()
+{
+  initlog();
+  test<2, 3>();
+}