Add treatments for diamond cells

Note that VTK output does not support diamond cells and only deals
with pyramids based on quadrangles

src/mesh/CellType.hpp

@@ -12,10 +12,11 @@ enum class CellType : unsigned short
   Triangle,
   Quadrangle,
 
-  Tetrahedron,
-  Pyramid,
+  Diamond,
+  Hexahedron,
   Prism,
-  Hexahedron
+  Pyramid,
+  Tetrahedron
 };
 
 PUGS_INLINE
@@ -29,14 +30,16 @@ name(CellType cell_type)
     return "triangle";
   case CellType::Quadrangle:
     return "quadrangle";
-  case CellType::Tetrahedron:
-    return "tetrahedron";
-  case CellType::Pyramid:
-    return "pyramid";
-  case CellType::Prism:
-    return "prism";
+  case CellType::Diamond:
+    return "diamond";
   case CellType::Hexahedron:
     return "hexahedron";
+  case CellType::Prism:
+    return "prism";
+  case CellType::Pyramid:
+    return "pyramid";
+  case CellType::Tetrahedron:
+    return "tetrahedron";
   default:
     return "unknown cell type";
   }

src/mesh/MeshBuilderBase.cpp

@@ -96,51 +96,95 @@ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescripto
       }
     } else if constexpr (Dimension == 3) {
       switch (descriptor.cell_type_vector[j]) {
-      case CellType::Tetrahedron: {
-        cell_nb_faces[j] = 4;
+      case CellType::Hexahedron: {
         // face 0
-        Face f0({cell_nodes[1], cell_nodes[2], cell_nodes[3]}, node_number_vector);
+        Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
         face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
 
         // face 1
-        Face f1({cell_nodes[0], cell_nodes[3], cell_nodes[2]}, node_number_vector);
+        Face f1({cell_nodes[4], cell_nodes[5], cell_nodes[6], cell_nodes[7]}, node_number_vector);
         face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
 
         // face 2
-        Face f2({cell_nodes[0], cell_nodes[1], cell_nodes[3]}, node_number_vector);
+        Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
         face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
 
         // face 3
-        Face f3({cell_nodes[0], cell_nodes[2], cell_nodes[1]}, node_number_vector);
+        Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
         face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+
+        // face 4
+        Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
+        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+
+        // face 5
+        Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
+        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+
+        cell_nb_faces[j] = 6;
         break;
       }
-      case CellType::Hexahedron: {
+      case CellType::Tetrahedron: {
+        cell_nb_faces[j] = 4;
         // face 0
-        Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
+        Face f0({cell_nodes[1], cell_nodes[2], cell_nodes[3]}, node_number_vector);
         face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
 
         // face 1
-        Face f1({cell_nodes[4], cell_nodes[5], cell_nodes[6], cell_nodes[7]}, node_number_vector);
+        Face f1({cell_nodes[0], cell_nodes[3], cell_nodes[2]}, node_number_vector);
         face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
 
         // face 2
-        Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
+        Face f2({cell_nodes[0], cell_nodes[1], cell_nodes[3]}, node_number_vector);
         face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
 
         // face 3
-        Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
+        Face f3({cell_nodes[0], cell_nodes[2], cell_nodes[1]}, node_number_vector);
         face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+        break;
+      }
+      case CellType::Pyramid: {
+        cell_nb_faces[j] = cell_nodes.size();
+        std::vector<unsigned int> base_nodes;
+        std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
 
-        // face 4
-        Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
-        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+        // base face
+        {
+          Face base_face(base_nodes, node_number_vector);
+          face_cells_map[base_face].emplace_back(std::make_tuple(j, 0, base_face.reversed()));
+        }
+        // side faces
+        const auto pyramid_vertex = cell_nodes[cell_nodes.size() - 1];
+        for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
+          Face side_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], pyramid_vertex},
+                         node_number_vector);
+          face_cells_map[side_face].emplace_back(std::make_tuple(j, i_node + 1, side_face.reversed()));
+        }
+        break;
+      }
+      case CellType::Diamond: {
+        cell_nb_faces[j] = 2 * (cell_nodes.size() - 2);
+        std::vector<unsigned int> base_nodes;
+        std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
 
-        // face 5
-        Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
-        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+        {   // top faces
+          const auto top_vertex = cell_nodes[cell_nodes.size() - 1];
+          for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
+            Face top_face({base_nodes[i_node], base_nodes[(i_node + 1) % base_nodes.size()], top_vertex},
+                          node_number_vector);
+            face_cells_map[top_face].emplace_back(std::make_tuple(j, i_node, top_face.reversed()));
+          }
+        }
 
-        cell_nb_faces[j] = 6;
+        {   // bottom faces
+          const auto bottom_vertex = cell_nodes[0];
+          for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
+            Face bottom_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], bottom_vertex},
+                             node_number_vector);
+            face_cells_map[bottom_face].emplace_back(
+              std::make_tuple(j, i_node + base_nodes.size(), bottom_face.reversed()));
+          }
+        }
         break;
       }
       default: {

src/mesh/MeshBuilderBase.hpp

@@ -44,7 +44,6 @@ class MeshBuilderBase::ConnectivityFace<2>
 {
  public:
   friend struct Hash;
-  friend class CartesianMeshBuilder;
 
   struct Hash
   {
@@ -122,10 +121,9 @@ class MeshBuilderBase::ConnectivityFace<2>
 template <>
 class MeshBuilderBase::ConnectivityFace<3>
 {
- private:
-  friend class GmshReader;
-  friend class CartesianMeshBuilder;
+ public:
   friend struct Hash;
+
   struct Hash
   {
     size_t