Added CGALWrappers::to_cgal_mesh()

cmake/config/template-arguments.in

@@ -277,3 +277,8 @@ SYM_RANKS := { 2; 4 }
 OUTPUT_FLAG_TYPES := { DXFlags; UcdFlags; GnuplotFlags; PovrayFlags; EpsFlags;
                        GmvFlags; TecplotFlags; VtkFlags; SvgFlags;
                        Deal_II_IntermediateFlags }
+
+// CGAL Kernels
+CGAL_KERNELS := {CGAL::Simple_cartesian<double>; CGAL::Exact_predicates_exact_constructions_kernel;
+                CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt; 
+                CGAL::Exact_predicates_inexact_constructions_kernel }

doc/news/changes/minor/20220428MarcoFederLucaHeltai

@@ -0,0 +1,4 @@
+New: Added function CGALWrappers::to_cgal_mesh() to convert a deal.II cell
+into a CGAL::Surface_mesh.
+<br>
+(Marco Feder, Luca Heltai, 2022/04/28)

include/deal.II/cgal/surface_mesh.h

@@ -0,0 +1,76 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_cgal_surface_mesh_h
+#define dealii_cgal_surface_mesh_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/fe/mapping.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/cgal/utilities.h>
+
+#ifdef DEAL_II_WITH_CGAL
+#  include <CGAL/Surface_mesh.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace CGALWrappers
+{
+  /**
+   * Build a CGAL::Surface_mesh from a deal.II cell.
+   *
+   * The class Surface_mesh implements a halfedge data structure and can be used
+   * to represent polyhedral surfaces. It is an edge-centered data structure
+   * capable of maintaining incidence information of vertices, edges, and faces.
+   * Each edge is represented by two halfedges with opposite orientation. The
+   * orientation of a face is chosen so that the halfedges around a face are
+   * oriented counterclockwise.
+   *
+   * More information on this class is available at
+   * https://doc.cgal.org/latest/Surface_mesh/index.html
+   *
+   * The function will throw an exception in dimension one. In dimension two it
+   * generates a surface mesh of the quadrilateral cell or of the triangle cell,
+   * while in dimension three it will generate the surface mesh of the cell,
+   * i.e., a polyhedral mesh containing the faces of the input cell.
+   *
+   * The generated mesh is useful when performing geometric operations using
+   * CGAL::Polygon_mesh_processing, i.e., to compute boolean operations on
+   * cells, splitting, cutting, slicing, etc.
+   *
+   * For examples on how to use the resulting  CGAL::Surface_mesh see
+   * https://doc.cgal.org/latest/Polygon_mesh_processing/
+   *
+   * @param[in] cell The input deal.II cell iterator
+   * @param[in] mapping The mapping used to map the vertices of the cell
+   * @param[out] mesh The output CGAL::Surface_mesh
+   */
+  template <typename CGALPointType, int dim, int spacedim>
+  void
+  to_cgal_mesh(
+    const typename dealii::Triangulation<dim, spacedim>::cell_iterator &cell,
+    const dealii::Mapping<dim, spacedim> &                              mapping,
+    CGAL::Surface_mesh<CGALPointType> &                                 mesh);
+} // namespace CGALWrappers
+
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+#endif

include/deal.II/cgal/utilities.h

@@ -20,8 +20,12 @@
 
 #ifdef DEAL_II_WITH_CGAL
 #  include <CGAL/Cartesian.h>
+#  include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#  include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
+#  include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #  include <CGAL/Simple_cartesian.h>
 
+
 DEAL_II_NAMESPACE_OPEN
 /**
  * Interface to the Computational Geometry Algorithm Library (CGAL).

source/CMakeLists.txt

@@ -56,6 +56,7 @@ ADD_SUBDIRECTORY(fe)
 ADD_SUBDIRECTORY(dofs)
 ADD_SUBDIRECTORY(lac)
 ADD_SUBDIRECTORY(base)
+ADD_SUBDIRECTORY(cgal)
 ADD_SUBDIRECTORY(gmsh)
 ADD_SUBDIRECTORY(grid)
 ADD_SUBDIRECTORY(hp)

source/cgal/CMakeLists.txt

@@ -0,0 +1,35 @@
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2012 - 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.
+##
+## ---------------------------------------------------------------------
+
+
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
+
+
+SET(_src
+surface_mesh.cc
+)
+
+
+
+SET(_inst
+surface_mesh.inst.in
+)
+
+FILE(GLOB _header
+  ${CMAKE_SOURCE_DIR}/include/deal.II/cgal/*.h
+  )
+
+DEAL_II_ADD_LIBRARY(obj_cgal OBJECT ${_src} ${_header} ${_inst})
+EXPAND_INSTANTIATIONS(obj_cgal "${_inst}")

source/cgal/surface_mesh.cc

@@ -0,0 +1,121 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+#include <deal.II/base/patterns.h>
+
+#include <deal.II/cgal/surface_mesh.h>
+
+#ifdef DEAL_II_WITH_CGAL
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+  template <typename dealiiFace, typename Container, typename CGAL_Mesh>
+  void
+  add_facet(const dealiiFace &face,
+            const Container & deal2cgal,
+            CGAL_Mesh &       mesh,
+            const bool        clockwise_ordering = true)
+  {
+    const unsigned                                nv = face->n_vertices();
+    std::vector<typename CGAL_Mesh::Vertex_index> indices;
+
+
+    switch (nv)
+      {
+        case 2:
+          mesh.add_edge(deal2cgal.at(face->vertex_index(0)),
+                        deal2cgal.at(face->vertex_index(1)));
+          break;
+        case 3:
+          indices = {deal2cgal.at(face->vertex_index(0)),
+                     deal2cgal.at(face->vertex_index(1)),
+                     deal2cgal.at(face->vertex_index(2))};
+          break;
+        case 4:
+          indices = {deal2cgal.at(face->vertex_index(0)),
+                     deal2cgal.at(face->vertex_index(1)),
+                     deal2cgal.at(face->vertex_index(3)),
+                     deal2cgal.at(face->vertex_index(2))};
+          break;
+        default:
+          Assert(false, ExcInternalError());
+          break;
+      }
+    auto f = mesh.null_face();
+    if (clockwise_ordering)
+      f = mesh.add_face(indices);
+    else
+      {
+        std::reverse(indices.begin(), indices.end());
+        f = mesh.add_face(indices);
+      }
+    Assert(f != mesh.null_face(),
+           ExcInternalError("While trying to build a CGAL facet, "
+                            "CGAL encountered a orientation problem that it "
+                            "was not able to solve."));
+  }
+} // namespace
+
+
+
+#  ifndef DOXYGEN
+// Template implementations
+namespace CGALWrappers
+{
+  template <typename CGALPointType, int dim, int spacedim>
+  void
+  to_cgal_mesh(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+               const Mapping<dim, spacedim> &     mapping,
+               CGAL::Surface_mesh<CGALPointType> &mesh)
+  {
+    Assert(dim > 1, ExcImpossibleInDim(dim));
+    using Mesh           = CGAL::Surface_mesh<CGALPointType>;
+    using Vertex_index   = typename Mesh::Vertex_index;
+    const auto &vertices = mapping.get_vertices(cell);
+
+    std::map<unsigned int, Vertex_index> deal2cgal;
+    // Add all vertices to the mesh
+    // Store CGAL ordering
+    for (const auto &i : cell->vertex_indices())
+      deal2cgal[cell->vertex_index(i)] =
+        mesh.add_vertex(CGALWrappers::to_cgal<CGALPointType>(vertices[i]));
+
+    // Add faces
+    if (dim < 3)
+      // simplices and quads
+      add_facet(cell, deal2cgal, mesh);
+    else
+      // faces of 3d cells
+      for (const auto &f : cell->face_indices())
+        add_facet(cell->face(f),
+                  deal2cgal,
+                  mesh,
+                  (f % 2 == 0 || cell->n_vertices() != 8));
+  }
+  // explicit instantiations
+#    include "surface_mesh.inst"
+
+
+} // namespace CGALWrappers
+#  endif
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

source/cgal/surface_mesh.inst.in

@@ -0,0 +1,26 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (dim : DIMENSIONS; spacedim : SPACE_DIMENSIONS; cgal_kernel : CGAL_KERNELS)
+  {
+#if dim <= spacedim
+    template void to_cgal_mesh<typename cgal_kernel::Point_3, dim, spacedim>(
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+      const Mapping<dim, spacedim> &                              mapping,
+      CGAL::Surface_mesh<typename cgal_kernel::Point_3> &         mesh);
+#endif
+  }

tests/cgal/cgal_surface_mesh_01.cc

@@ -0,0 +1,72 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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 deal.II cell to a cgal Surface_mesh
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <CGAL/IO/io.h>
+#include <deal.II/cgal/surface_mesh.h>
+
+#include "../tests.h"
+
+
+using namespace CGALWrappers;
+using CGALPoint = CGAL::Point_3<CGAL::Simple_cartesian<double>>;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  deallog << "dim= " << dim << ",\t spacedim= " << spacedim << std::endl;
+  std::vector<std::vector<unsigned int>> d2t = {{}, {2}, {3, 4}, {4, 5, 6, 8}};
+  for (const auto nv : d2t[dim])
+    {
+      Triangulation<dim, spacedim>  tria;
+      CGAL::Surface_mesh<CGALPoint> mesh;
+
+      const auto ref     = ReferenceCell::n_vertices_to_type(dim, nv);
+      const auto mapping = ref.template get_default_mapping<dim, spacedim>(1);
+      GridGenerator::reference_cell(tria, ref);
+
+      const auto cell = tria.begin_active();
+      to_cgal_mesh(cell, *mapping, mesh);
+
+      Assert(mesh.is_valid(), dealii::ExcMessage("The CGAL mesh is not valid"));
+      deallog << "deal vertices: " << nv << ", cgal vertices "
+              << mesh.num_vertices() << std::endl;
+      deallog << "deal faces: " << cell->n_faces() << ", cgal faces "
+              << mesh.num_faces() << std::endl;
+      deallog << "Valid mesh: " << std::boolalpha << mesh.is_valid()
+              << std::endl;
+      deallog << mesh << std::endl;
+    }
+}
+
+int
+main()
+{
+  initlog();
+  test<2, 2>();
+  test<2, 3>();
+  test<3, 3>();
+}