Fix pyramid node ordering when building the diamond mesh

src/mesh/DiamondDualConnectivityBuilder.cpp

@@ -86,21 +86,27 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
   diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);
 
   const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
+  const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
 
+  static_assert(Dimension > 1);
   parallel_for(primal_number_of_faces, [&](FaceId face_id) {
     const size_t i_cell               = face_id;
     const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
 
+    if constexpr (Dimension == 2) {
       if (primal_face_cell_list.size() == 1) {
         diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
       } else {
         Assert(primal_face_cell_list.size() == 2);
         diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
       }
-
-    if constexpr (Dimension == 3) {
+    } else if constexpr (Dimension == 3) {
       if (primal_face_cell_list.size() == 1) {
+        if (primal_face_to_node_matrix[face_id].size() == 3) {
+          diamond_descriptor.cell_type_vector[i_cell] = CellType::Tetrahedron;
+        } else {
           diamond_descriptor.cell_type_vector[i_cell] = CellType::Pyramid;
+        }
       } else {
         Assert(primal_face_cell_list.size() == 2);
         diamond_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
@@ -110,7 +116,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
 
   diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
 
-  const auto& primal_face_to_node_matrix              = primal_connectivity.faceToNodeMatrix();
   const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
   const auto& cell_face_is_reversed                   = primal_connectivity.cellFaceIsReversed();
   parallel_for(primal_number_of_faces, [&](FaceId face_id) {
@@ -129,12 +134,17 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
       diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size()] =
         primal_number_of_nodes + cell_id;
 
-      if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
       if constexpr (Dimension == 2) {
+        if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
           std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
                     diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]);
-
+        }
       } else {
+        if (not cell_face_is_reversed(cell_id, i_face_in_cell)) {
+          // In 3D the basis of the pyramid is described in the
+          // indirect way IF the face is not reversed. In other words
+          // the "topological normal" must point to the "top" of the
+          // pyramid.
           for (size_t i_node = 0; i_node < primal_face_node_list.size() / 2; ++i_node) {
             std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node],
                       diamond_descriptor