Code review.

doc/news/changes/major/20220426FederHeltai

@@ -1,9 +1,9 @@
 New: Added support for the CGAL library (www.cgal.org). The following features 
 are supported:
 <ul>
-  <li>Conversion from dealii to CGAL point types and viceversa</li>
-  <li>Conversion from dealii cell types to CGAL::Surface_mesh and viceversa</li>
-  <li>Conversion from dealii Triangulation to CGAL::Surface_mesh and viceversa</li>
+  <li>Conversion from deal.II to CGAL point types and viceversa</li>
+  <li>Conversion from deal.II cell types to CGAL::Surface_mesh and viceversa</li>
+  <li>Conversion from deal.II Triangulation to CGAL::Surface_mesh and viceversa</li>
   <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>

include/deal.II/cgal/triangulation.h

@@ -18,15 +18,15 @@
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/grid/tria.h>
+#ifdef DEAL_II_WITH_CGAL
 
-#include <deal.II/cgal/utilities.h>
+#  include <deal.II/grid/tria.h>
 
-#ifdef DEAL_II_WITH_CGAL
 #  include <boost/hana.hpp>
 
 #  include <CGAL/Surface_mesh.h>
 #  include <CGAL/Triangulation_3.h>
+#  include <deal.II/cgal/utilities.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -95,8 +95,8 @@ namespace CGALWrappers
    * input Triangulation `dealii_triangulation`. If this is not the case, an
    * exception is thrown.
    *
-   * @param cgal_triangulation
-   * @param dealii_triangulation
+   * @param[in] cgal_triangulation The input CGAL triangulation.
+   * @param[out] dealii_triangulation The output deal.II triangulation.
    */
   template <typename CGALTriangulation, int dim, int spacedim>
   void
@@ -117,12 +117,15 @@ namespace CGALWrappers
            ExcMessage(
              "The triangulation you pass to this function should be a valid "
              "CGAL triangulation."));
-    using CGALPoint = typename CGALTriangulation::Point;
-    std::vector<CGALPoint> cgal_points(points.size());
-    std::transform(
-      points.begin(), points.end(), cgal_points.begin(), [](const auto &p) {
-        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint>(p);
-      });
+
+    std::vector<typename CGALTriangulation::Point> cgal_points(points.size());
+    std::transform(points.begin(),
+                   points.end(),
+                   cgal_points.begin(),
+                   [](const auto &p) {
+                     return CGALWrappers::dealii_point_to_cgal_point<
+                       typename CGALTriangulation::Point>(p);
+                   });
 
     triangulation.insert(cgal_points.begin(), cgal_points.end());
     Assert(triangulation.is_valid(),
@@ -149,19 +152,18 @@ namespace CGALWrappers
     Assert(dealii_triangulation.n_cells() == 0,
            ExcMessage("The output triangulation object needs to be empty."));
 
-    const auto nv = cgal_triangulation.number_of_vertices();
-
-    // Deal.II storage data structures
-    std::vector<Point<spacedim>> vertices(nv);
-    std::vector<CellData<dim>>   cells;
-    SubCellData                  subcell_data;
+    // deal.II storage data structures
+    std::vector<Point<spacedim>> vertices(
+      cgal_triangulation.number_of_vertices());
+    std::vector<CellData<dim>> cells;
+    SubCellData                subcell_data;
 
     // CGAL storage data structures
     std::map<typename CGALTriangulation::Vertex_handle, unsigned int>
       vertex_map;
     {
       unsigned int i = 0;
-      for (auto v : cgal_triangulation.finite_vertex_handles())
+      for (const auto &v : cgal_triangulation.finite_vertex_handles())
         {
           vertices[i] =
             CGALWrappers::cgal_point_to_dealii_point<spacedim>(v->point());
@@ -182,8 +184,10 @@ namespace CGALWrappers
         if (cgal_triangulation.dimension() == 2)
           for (const auto f : cgal_triangulation.finite_face_handles())
             {
-              CellData<dim> cell(3);
-              for (unsigned int i = 0; i < 3; ++i)
+              CellData<dim> cell(ReferenceCells::Triangle.n_vertices());
+              for (unsigned int i = 0;
+                   i < ReferenceCells::Triangle.n_vertices();
+                   ++i)
                 cell.vertices[i] = vertex_map[f->vertex(i)];
               cells.push_back(cell);
             }
@@ -194,9 +198,16 @@ namespace CGALWrappers
               // An edge is idenfied by a face and a vertex index in the face
               const auto &  f = e.first;
               const auto &  i = e.second;
-              CellData<dim> cell(2);
+              CellData<dim> cell(ReferenceCells::Line.n_vertices());
               unsigned int  id = 0;
-              for (unsigned int j = 0; j < 3; ++j)
+              // 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.
+              for (unsigned int j = 0;
+                   j < ReferenceCells::Triangle.n_vertices();
+                   ++j)
                 if (j != i)
                   cell.vertices[id++] = vertex_map[f->vertex(j)];
               cells.push_back(cell);
@@ -212,8 +223,10 @@ namespace CGALWrappers
         if (cgal_triangulation.dimension() == 3)
           for (const auto c : cgal_triangulation.finite_cell_handles())
             {
-              CellData<dim> cell(4);
-              for (unsigned int i = 0; i < 4; ++i)
+              CellData<dim> cell(ReferenceCells::Tetrahedron.n_vertices());
+              for (unsigned int i = 0;
+                   i < ReferenceCells::Tetrahedron.n_vertices();
+                   ++i)
                 cell.vertices[i] = vertex_map[c->vertex(i)];
               cells.push_back(cell);
             }
@@ -225,9 +238,16 @@ namespace CGALWrappers
               // the face
               const auto &  c = facet.first;
               const auto &  i = facet.second;
-              CellData<dim> cell(3);
+              CellData<dim> cell(ReferenceCells::Triangle.n_vertices());
               unsigned int  id = 0;
-              for (unsigned int j = 0; j < 4; ++j)
+              // 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.
+              for (unsigned int j = 0;
+                   j < ReferenceCells::Tetrahedron.n_vertices();
+                   ++j)
                 if (j != i)
                   cell.vertices[id++] = vertex_map[c->vertex(j)];
               cells.push_back(cell);
@@ -238,7 +258,7 @@ namespace CGALWrappers
             {
               // An edge is idenfied by a cell and its two vertices
               const auto &[c, i, j] = edge;
-              CellData<dim> cell(2);
+              CellData<dim> cell(ReferenceCells::Line.n_vertices());
               cell.vertices[0] = vertex_map[c->vertex(i)];
               cell.vertices[1] = vertex_map[c->vertex(j)];
               cells.push_back(cell);