diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 9c5cf3aec1e11d16f652357db53f68db54f6ed79..af1fca45e7dec5898a12731430358691728c7a2a 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -12,6 +12,7 @@
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/CartesianMeshBuilder.hpp>
 #include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
 #include <mesh/DualMeshManager.hpp>
 #include <mesh/GmshReader.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
@@ -143,6 +144,37 @@ MeshModule::MeshModule()
 
                                        ));
 
+  this->_addBuiltinFunction("check_connectivity_ordering",
+                            std::function(
+
+                              [](const std::shared_ptr<const IMesh>& i_mesh) -> bool {
+                                switch (i_mesh->dimension()) {
+                                case 1: {
+                                  using MeshType = Mesh<Connectivity<1>>;
+
+                                  std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+                                  return checkConnectivityOrdering(p_mesh->connectivity());
+                                }
+                                case 2: {
+                                  using MeshType = Mesh<Connectivity<2>>;
+
+                                  std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+                                  return checkConnectivityOrdering(p_mesh->connectivity());
+                                }
+                                case 3: {
+                                  using MeshType = Mesh<Connectivity<3>>;
+
+                                  std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+                                  return checkConnectivityOrdering(p_mesh->connectivity());
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid dimension");
+                                }
+                                }
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("cartesianMesh",
                             std::function(
 
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 66ce42ae6f9e490a0b3e6e406934478b0f026f4a..62dbd50536ad2e725e67a25d9cddf481a62b40c6 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -7,6 +7,7 @@ add_library(
   ConnectivityBuilderBase.cpp
   ConnectivityComputer.cpp
   ConnectivityDispatcher.cpp
+  ConnectivityUtils.cpp
   DiamondDualConnectivityBuilder.cpp
   DiamondDualMeshBuilder.cpp
   Dual1DConnectivityBuilder.cpp
diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
index 79fb0a2011bceae67a20e7b877f2896a2af4e66e..9f6ea50cab88390669e3471ad3c70f25e870e390 100644
--- a/src/mesh/CartesianMeshBuilder.cpp
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -5,8 +5,6 @@
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
 
-#include <utils/Timer.hpp>
-
 template <size_t Dimension>
 NodeValue<TinyVector<Dimension>>
 CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<Dimension>&,
@@ -69,9 +67,6 @@ CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<3>& a,
                                           const TinyVector<3, uint64_t>& cell_size,
                                           const IConnectivity& connectivity) const
 {
-  Timer t;
-  t.start();
-
   const TinyVector<3, uint64_t> node_size = [&] {
     TinyVector node_size{cell_size};
     for (size_t i = 0; i < 3; ++i) {
@@ -100,8 +95,6 @@ CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<3>& a,
       }
     });
 
-  std::cout << "CartesianMeshBuilder::_getNodeCoordinates t = " << t << '\n';
-
   return xr;
 }
 
@@ -111,8 +104,6 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<Dimension>& a,
                                           const TinyVector<Dimension>& b,
                                           const TinyVector<Dimension, uint64_t>& cell_size)
 {
-  Timer t;
-  t.start();
   static_assert(Dimension <= 3, "unexpected dimension");
 
   LogicalConnectivityBuilder logical_connectivity_builder{cell_size};
@@ -126,8 +117,6 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<Dimension>& a,
   NodeValue<TinyVector<Dimension>> xr = _getNodeCoordinates(a, b, cell_size, connectivity);
 
   m_mesh = std::make_shared<Mesh<ConnectivityType>>(p_connectivity, xr);
-
-  std::cout << "_buildCartesianMesh t = " << t << '\n';
 }
 
 template <size_t Dimension>
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 9dd3bf5463b094ba556895ba1fd84c2f20d74bd8..4cb11d443516aea5c02c329baa1eb5ee9b968e0b 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -13,11 +13,11 @@ template <size_t Dimension>
 void
 Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
 {
-  Assert(descriptor.cell_to_node_vector.size() == descriptor.cell_type_vector.size());
+  Assert(descriptor.cell_to_node_matrix.numberOfRows() == descriptor.cell_type_vector.size());
   Assert(descriptor.cell_number_vector.size() == descriptor.cell_type_vector.size());
   if constexpr (Dimension > 1) {
-    Assert(descriptor.cell_to_face_vector.size() == descriptor.cell_type_vector.size());
-    Assert(descriptor.face_to_node_vector.size() == descriptor.face_number_vector.size());
+    Assert(descriptor.cell_to_face_matrix.numberOfRows() == descriptor.cell_type_vector.size());
+    Assert(descriptor.face_to_node_matrix.numberOfRows() == descriptor.face_number_vector.size());
     Assert(descriptor.face_owner_vector.size() == descriptor.face_number_vector.size());
   }
 
@@ -38,7 +38,7 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
   }
 
   auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
-  cell_to_node_matrix       = descriptor.cell_to_node_vector;
+  cell_to_node_matrix       = descriptor.cell_to_node_matrix;
 
   {
     WeakCellValue<CellType> cell_type(*this);
@@ -112,20 +112,11 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
     }
 
   } else {
-    m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_vector;
+    m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_matrix;
 
-    m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_vector;
+    m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_matrix;
 
-    {
-      FaceValuePerCell<bool> cell_face_is_reversed(*this);
-      for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
-        const auto& face_cells_vector = descriptor.cell_face_is_reversed_vector[j];
-        for (unsigned short lj = 0; lj < face_cells_vector.size(); ++lj) {
-          cell_face_is_reversed(j, lj) = face_cells_vector[lj];
-        }
-      }
-      m_cell_face_is_reversed = cell_face_is_reversed;
-    }
+    m_cell_face_is_reversed = FaceValuePerCell<bool>(*this, descriptor.cell_face_is_reversed);
 
     m_face_number = WeakFaceValue<int>(*this, descriptor.face_number_vector);
 
@@ -162,22 +153,13 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
       }
 
     } else {
-      m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_vector;
+      m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_matrix;
 
-      m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.face_to_edge_vector;
+      m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.face_to_edge_matrix;
 
-      m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cell_to_edge_vector;
+      m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cell_to_edge_matrix;
 
-      {
-        EdgeValuePerFace<bool> face_edge_is_reversed(*this);
-        for (FaceId l = 0; l < descriptor.face_edge_is_reversed_vector.size(); ++l) {
-          const auto& edge_faces_vector = descriptor.face_edge_is_reversed_vector[l];
-          for (unsigned short el = 0; el < edge_faces_vector.size(); ++el) {
-            face_edge_is_reversed(l, el) = edge_faces_vector[el];
-          }
-        }
-        m_face_edge_is_reversed = face_edge_is_reversed;
-      }
+      m_face_edge_is_reversed = EdgeValuePerFace<bool>(*this, descriptor.face_edge_is_reversed);
 
       m_edge_number = WeakEdgeValue<int>(*this, descriptor.edge_number_vector);
       m_edge_owner  = WeakEdgeValue<int>(*this, descriptor.edge_owner_vector);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index b0279507ee41e6f6d3433bda2b5dd7beaf406265..ee7b07fe043a44262c23afc053e748bcac448b50 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -122,7 +122,7 @@ class Connectivity final : public IConnectivity
     const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
     if (not connectivity_matrix.isBuilt()) {
       const_cast<ConnectivityMatrix&>(connectivity_matrix) =
-        m_connectivity_computer.computeConnectivityMatrix(*this, item_type_0, item_type_1);
+        m_connectivity_computer.computeInverseConnectivityMatrix(*this, item_type_0, item_type_1);
     }
     return connectivity_matrix;
   }
diff --git a/src/mesh/ConnectivityBuilderBase.cpp b/src/mesh/ConnectivityBuilderBase.cpp
index 1657507ec030a95b6ba5238d091c4426e1764f9f..a9b99e5218a8d779552dbc208c5cd9c39955005b 100644
--- a/src/mesh/ConnectivityBuilderBase.cpp
+++ b/src/mesh/ConnectivityBuilderBase.cpp
@@ -19,8 +19,8 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
 
   size_t total_number_of_faces = 0;
 
-  for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
-    const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+  for (CellId j = 0; j < descriptor.cell_to_node_matrix.numberOfRows(); ++j) {
+    const auto& cell_nodes = descriptor.cell_to_node_matrix[j];
 
     if constexpr (Dimension == 2) {
       total_number_of_faces += cell_nodes.size();
@@ -61,7 +61,7 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
     } else {
       Assert(Dimension == 3);
       size_t count_number_of_face_by_node = 0;
-      for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
+      for (CellId j = 0; j < descriptor.cell_to_node_matrix.numberOfRows(); ++j) {
         switch (descriptor.cell_type_vector[j]) {
         case CellType::Hexahedron: {
           count_number_of_face_by_node += 6 * 4;
@@ -76,12 +76,12 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
           break;
         }
         case CellType::Pyramid: {
-          const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+          const auto& cell_nodes = descriptor.cell_to_node_matrix[j];
           count_number_of_face_by_node += 1 * cell_nodes.size() + cell_nodes.size() * 3;
           break;
         }
         case CellType::Diamond: {
-          const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+          const auto& cell_nodes = descriptor.cell_to_node_matrix[j];
           count_number_of_face_by_node += cell_nodes.size() * 3 * 2;
           break;
         }
@@ -96,65 +96,65 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
     }
   }();
 
-  Array<unsigned int> faces_node_list(total_number_of_face_by_node);
-  Array<unsigned int> face_ending(total_number_of_faces + 1);
-  size_t i_face  = 0;
-  face_ending[0] = 0;
+  Array<unsigned int> dup_faces_to_node_list(total_number_of_face_by_node);
+  Array<unsigned int> dup_face_to_node_row(total_number_of_faces + 1);
+  size_t i_face           = 0;
+  dup_face_to_node_row[0] = 0;
 
-  Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
+  Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_matrix.numberOfRows());
   {
     size_t i_face_node = 0;
-    for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
-      const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+    for (CellId j = 0; j < descriptor.cell_to_node_matrix.numberOfRows(); ++j) {
+      const auto& cell_nodes = descriptor.cell_to_node_matrix[j];
 
       if constexpr (Dimension == 2) {
         switch (descriptor.cell_type_vector[j]) {
         case CellType::Triangle: {
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-
-          face_ending[i_face + 1] = face_ending[i_face] + 2;
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
           i_face++;
-          face_ending[i_face + 1] = face_ending[i_face] + 2;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
           i_face++;
-          face_ending[i_face + 1] = face_ending[i_face] + 2;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
           i_face++;
           cell_nb_faces[j] = 3;
           break;
         }
         case CellType::Quadrangle: {
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[3];
-          faces_node_list[i_face_node++] = cell_nodes[3];
-          faces_node_list[i_face_node++] = cell_nodes[0];
-
-          face_ending[i_face + 1] = face_ending[i_face] + 2;
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
           i_face++;
-          face_ending[i_face + 1] = face_ending[i_face] + 2;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
           i_face++;
-          face_ending[i_face + 1] = face_ending[i_face] + 2;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
           i_face++;
-          face_ending[i_face + 1] = face_ending[i_face] + 2;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
           i_face++;
           cell_nb_faces[j] = 4;
           break;
         }
         case CellType::Polygon: {
           for (size_t i = 0; i < cell_nodes.size(); ++i) {
-            faces_node_list[i_face_node] = cell_nodes[i];
+            dup_faces_to_node_list[i_face_node] = cell_nodes[i];
             i_face_node++;
-            faces_node_list[i_face_node] = cell_nodes[(i + 1) % cell_nodes.size()];
+            dup_faces_to_node_list[i_face_node] = cell_nodes[(i + 1) % cell_nodes.size()];
             i_face_node++;
 
-            face_ending[i_face + 1] = face_ending[i_face] + 2;
+            dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
             i_face++;
           }
 
@@ -170,126 +170,126 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
       } else if constexpr (Dimension == 3) {
         switch (descriptor.cell_type_vector[j]) {
         case CellType::Hexahedron: {
-          faces_node_list[i_face_node++] = cell_nodes[3];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[4];
-          faces_node_list[i_face_node++] = cell_nodes[5];
-          faces_node_list[i_face_node++] = cell_nodes[6];
-          faces_node_list[i_face_node++] = cell_nodes[7];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[6];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[7];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[4];
-          faces_node_list[i_face_node++] = cell_nodes[7];
-          faces_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[7];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[6];
-          faces_node_list[i_face_node++] = cell_nodes[5];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[6];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[5];
-          faces_node_list[i_face_node++] = cell_nodes[4];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[3];
-          faces_node_list[i_face_node++] = cell_nodes[7];
-          faces_node_list[i_face_node++] = cell_nodes[6];
-          faces_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[7];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[6];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
           cell_nb_faces[j] = 6;
           break;
         }
         case CellType::Tetrahedron: {
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 3;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[3];
-          faces_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 3;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 3;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 3;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
           i_face++;
 
           cell_nb_faces[j] = 4;
           break;
         }
         case CellType::Prism: {
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 3;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[3];
-          faces_node_list[i_face_node++] = cell_nodes[4];
-          faces_node_list[i_face_node++] = cell_nodes[5];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 3;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[5];
-          faces_node_list[i_face_node++] = cell_nodes[4];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[1];
-          faces_node_list[i_face_node++] = cell_nodes[4];
-          faces_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
-          faces_node_list[i_face_node++] = cell_nodes[2];
-          faces_node_list[i_face_node++] = cell_nodes[0];
-          faces_node_list[i_face_node++] = cell_nodes[3];
-          faces_node_list[i_face_node++] = cell_nodes[5];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+          dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
 
-          face_ending[i_face + 1] = face_ending[i_face] + 4;
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
           i_face++;
 
           cell_nb_faces[j] = 5;
@@ -303,48 +303,47 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
           }
 
           for (size_t i = 0; i < base_nodes.size(); ++i) {
-            faces_node_list[i_face_node++] = base_nodes[i];
+            dup_faces_to_node_list[i_face_node++] = base_nodes[i];
           }
 
-          face_ending[i_face + 1] = face_ending[i_face] + base_nodes.size();
+          dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + base_nodes.size();
           i_face++;
 
           // 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) {
-            faces_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
-            faces_node_list[i_face_node++] = base_nodes[i_node];
-            faces_node_list[i_face_node++] = pyramid_vertex;
+            dup_faces_to_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
+            dup_faces_to_node_list[i_face_node++] = base_nodes[i_node];
+            dup_faces_to_node_list[i_face_node++] = pyramid_vertex;
 
-            face_ending[i_face + 1] = face_ending[i_face] + 3;
+            dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
             i_face++;
           }
           break;
         }
         case CellType::Diamond: {
-          std::vector<unsigned int> base_nodes;
-          std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
+          auto base_nodes = [&](size_t i) { return cell_nodes[i + 1]; };
 
           {   // 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) {
-              faces_node_list[i_face_node++] = base_nodes[i_node];
-              faces_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
-              faces_node_list[i_face_node++] = top_vertex;
+            for (size_t i_node = 0; i_node < cell_nodes.size() - 2; ++i_node) {
+              dup_faces_to_node_list[i_face_node++] = base_nodes(i_node);
+              dup_faces_to_node_list[i_face_node++] = base_nodes((i_node + 1) % (cell_nodes.size() - 2));
+              dup_faces_to_node_list[i_face_node++] = top_vertex;
 
-              face_ending[i_face + 1] = face_ending[i_face] + 3;
+              dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
               i_face++;
             }
           }
 
           {   // bottom faces
             const auto bottom_vertex = cell_nodes[0];
-            for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
-              faces_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
-              faces_node_list[i_face_node++] = base_nodes[i_node];
-              faces_node_list[i_face_node++] = bottom_vertex;
+            for (size_t i_node = 0; i_node < cell_nodes.size() - 2; ++i_node) {
+              dup_faces_to_node_list[i_face_node++] = base_nodes((i_node + 1) % (cell_nodes.size() - 2));
+              dup_faces_to_node_list[i_face_node++] = base_nodes(i_node);
+              dup_faces_to_node_list[i_face_node++] = bottom_vertex;
 
-              face_ending[i_face + 1] = face_ending[i_face] + 3;
+              dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
               i_face++;
             }
           }
@@ -361,26 +360,27 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
     }
   }
 
-  Array<bool> face_reversed(total_number_of_faces);
+  Array<bool> cell_face_is_reversed(total_number_of_faces);
 
   if constexpr (Dimension == 2) {
     for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
-      if (node_number_vector[faces_node_list[2 * i_face]] > node_number_vector[faces_node_list[2 * i_face + 1]]) {
-        std::swap(faces_node_list[2 * i_face], faces_node_list[2 * i_face + 1]);
-        face_reversed[i_face] = true;
+      if (node_number_vector[dup_faces_to_node_list[2 * i_face]] >
+          node_number_vector[dup_faces_to_node_list[2 * i_face + 1]]) {
+        std::swap(dup_faces_to_node_list[2 * i_face], dup_faces_to_node_list[2 * i_face + 1]);
+        cell_face_is_reversed[i_face] = true;
       } else {
-        face_reversed[i_face] = false;
+        cell_face_is_reversed[i_face] = false;
       }
     }
   } else if constexpr (Dimension == 3) {
     std::vector<int> buffer;
 
     for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
-      const size_t face_node_number      = face_ending[i_face + 1] - face_ending[i_face];
+      const size_t face_node_number      = dup_face_to_node_row[i_face + 1] - dup_face_to_node_row[i_face];
       size_t i_face_node_smallest_number = 0;
       for (size_t i_face_node = 1; i_face_node < face_node_number; ++i_face_node) {
-        if (node_number_vector[faces_node_list[face_ending[i_face] + i_face_node]] <
-            node_number_vector[faces_node_list[face_ending[i_face] + i_face_node_smallest_number]]) {
+        if (node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_face_node]] <
+            node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_face_node_smallest_number]]) {
           i_face_node_smallest_number = i_face_node;
         }
       }
@@ -388,27 +388,28 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
       if (i_face_node_smallest_number != 0) {
         buffer.resize(face_node_number);
         for (size_t i_node = i_face_node_smallest_number; i_node < face_node_number; ++i_node) {
-          buffer[i_node - i_face_node_smallest_number] = faces_node_list[face_ending[i_face] + i_node];
+          buffer[i_node - i_face_node_smallest_number] = dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node];
         }
         for (size_t i_node = 0; i_node < i_face_node_smallest_number; ++i_node) {
           buffer[i_node + face_node_number - i_face_node_smallest_number] =
-            faces_node_list[face_ending[i_face] + i_node];
+            dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node];
         }
 
         for (size_t i_node = 0; i_node < face_node_number; ++i_node) {
-          faces_node_list[face_ending[i_face] + i_node] = buffer[i_node];
+          dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node] = buffer[i_node];
         }
       }
 
-      if (node_number_vector[faces_node_list[face_ending[i_face] + 1]] >
-          node_number_vector[faces_node_list[face_ending[i_face + 1] - 1]]) {
+      if (node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face] + 1]] >
+          node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face + 1] - 1]]) {
         for (size_t i_node = 1; i_node <= (face_node_number + 1) / 2 - 1; ++i_node) {
-          std::swap(faces_node_list[face_ending[i_face] + i_node], faces_node_list[face_ending[i_face + 1] - i_node]);
+          std::swap(dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node],
+                    dup_faces_to_node_list[dup_face_to_node_row[i_face + 1] - i_node]);
         }
 
-        face_reversed[i_face] = true;
+        cell_face_is_reversed[i_face] = true;
       } else {
-        face_reversed[i_face] = false;
+        cell_face_is_reversed[i_face] = false;
       }
     }
   }
@@ -416,8 +417,9 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
   Array<unsigned int> node_to_duplicate_face_row(node_number_vector.size() + 1);
   node_to_duplicate_face_row.fill(0);
   for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
-    for (size_t i_node_face = face_ending[i_face]; i_node_face < face_ending[i_face + 1]; ++i_node_face) {
-      node_to_duplicate_face_row[faces_node_list[i_node_face] + 1]++;
+    for (size_t i_node_face = dup_face_to_node_row[i_face]; i_node_face < dup_face_to_node_row[i_face + 1];
+         ++i_node_face) {
+      node_to_duplicate_face_row[dup_faces_to_node_list[i_node_face] + 1]++;
     }
   }
 
@@ -428,14 +430,15 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
   Array<unsigned int> node_duplicate_face_list(node_to_duplicate_face_row[node_to_duplicate_face_row.size() - 1]);
 
   {
-    Array<unsigned int> node_face_row_idx(node_number_vector.size());
-    node_face_row_idx.fill(0);
+    Array<unsigned int> node_duplicate_face_row_idx(node_number_vector.size());
+    node_duplicate_face_row_idx.fill(0);
 
     for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
-      for (size_t i_node_face = face_ending[i_face]; i_node_face < face_ending[i_face + 1]; ++i_node_face) {
-        const size_t node_id = faces_node_list[i_node_face];
-        node_duplicate_face_list[node_to_duplicate_face_row[node_id] + node_face_row_idx[node_id]] = i_face;
-        node_face_row_idx[node_id]++;
+      for (size_t i_node_face = dup_face_to_node_row[i_face]; i_node_face < dup_face_to_node_row[i_face + 1];
+           ++i_node_face) {
+        const size_t node_id = dup_faces_to_node_list[i_node_face];
+        node_duplicate_face_list[node_to_duplicate_face_row[node_id] + node_duplicate_face_row_idx[node_id]] = i_face;
+        node_duplicate_face_row_idx[node_id]++;
       }
     }
   }
@@ -445,13 +448,14 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
     total_number_of_faces, PUGS_LAMBDA(size_t i_face) { dup_face_to_face[i_face] = i_face; });
 
   auto is_same_face = [=](const size_t i_face, const size_t j_face) {
-    if ((face_ending[i_face + 1] - face_ending[i_face]) != (face_ending[j_face + 1] - face_ending[j_face])) {
+    if ((dup_face_to_node_row[i_face + 1] - dup_face_to_node_row[i_face]) !=
+        (dup_face_to_node_row[j_face + 1] - dup_face_to_node_row[j_face])) {
       return false;
     } else {
-      auto i_face_node = face_ending[i_face];
-      auto j_face_node = face_ending[j_face];
-      while (i_face_node < face_ending[i_face + 1]) {
-        if (faces_node_list[i_face_node] != faces_node_list[j_face_node]) {
+      auto i_face_node = dup_face_to_node_row[i_face];
+      auto j_face_node = dup_face_to_node_row[j_face];
+      while (i_face_node < dup_face_to_node_row[i_face + 1]) {
+        if (dup_faces_to_node_list[i_face_node] != dup_faces_to_node_list[j_face_node]) {
           return false;
         }
         i_face_node++;
@@ -496,21 +500,54 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
     }
   }
 
-  Array<unsigned int> cell_face_begin(cell_nb_faces.size() + 1);
-  cell_face_begin[0] = 0;
-  for (size_t i = 0; i < cell_nb_faces.size(); ++i) {
-    cell_face_begin[i + 1] = cell_face_begin[i] + cell_nb_faces[i];
+  Array<unsigned int> node_to_face_row(node_to_duplicate_face_row.size());
+  {
+    unsigned int nb_faces = 0;
+    for (size_t i_node = 0; i_node < node_to_duplicate_face_row.size() - 1; ++i_node) {
+      node_to_face_row[i_node] = nb_faces;
+      for (size_t i_face = node_to_duplicate_face_row[i_node]; i_face < node_to_duplicate_face_row[i_node + 1];
+           ++i_face) {
+        if (dup_face_to_face[node_duplicate_face_list[i_face]] == node_duplicate_face_list[i_face]) {
+          ++nb_faces;
+        }
+      }
+    }
+    node_to_face_row[node_to_duplicate_face_row.size() - 1] = nb_faces;
+  }
+
+  Array<unsigned int> node_to_face_list(node_to_face_row[node_to_duplicate_face_row.size() - 1]);
+  {
+    unsigned int i_node_to_face = 0;
+    for (size_t i_node = 0; i_node < node_to_duplicate_face_row.size() - 1; ++i_node) {
+      for (size_t i_face = node_to_duplicate_face_row[i_node]; i_face < node_to_duplicate_face_row[i_node + 1];
+           ++i_face) {
+        if (dup_face_to_face[node_duplicate_face_list[i_face]] == node_duplicate_face_list[i_face]) {
+          node_to_face_list[i_node_to_face++] = dup_face_to_face_id[node_duplicate_face_list[i_face]];
+        }
+      }
+    }
   }
 
-  descriptor.cell_to_face_vector.resize(descriptor.cell_to_node_vector.size());
-  for (CellId j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
-    descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
+  descriptor.node_to_face_matrix = ConnectivityMatrix(node_to_face_row, node_to_face_list);
 
-    for (size_t i_face = 0; i_face < cell_nb_faces[j]; ++i_face) {
-      descriptor.cell_to_face_vector[j][i_face] = dup_face_to_face_id[cell_face_begin[j] + i_face];
+  Array<unsigned int> cell_to_face_row(cell_nb_faces.size() + 1);
+  cell_to_face_row[0] = 0;
+  for (size_t i = 0; i < cell_nb_faces.size(); ++i) {
+    cell_to_face_row[i + 1] = cell_to_face_row[i] + cell_nb_faces[i];
+  }
+
+  Array<unsigned int> cell_to_face_list(cell_to_face_row[cell_to_face_row.size() - 1]);
+  {
+    size_t i_cell_face = 0;
+    for (CellId cell_id = 0; cell_id < cell_nb_faces.size(); ++cell_id) {
+      for (size_t i_face = 0; i_face < cell_nb_faces[cell_id]; ++i_face) {
+        cell_to_face_list[i_cell_face++] = dup_face_to_face_id[cell_to_face_row[cell_id] + i_face];
+      }
     }
   }
 
+  descriptor.cell_to_face_matrix = ConnectivityMatrix(cell_to_face_row, cell_to_face_list);
+
   {
     // Face numbers may change if numbers are provided in the file
     descriptor.face_number_vector = [&] {
@@ -522,33 +559,40 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD
     }();
   }
 
+  Array<unsigned int> face_ending(total_number_of_faces - nb_dup_faces + 1);
   {
-    descriptor.face_to_node_vector.resize(total_number_of_faces - nb_dup_faces);
-    int l = 0;
-
-    for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
-      if (dup_face_to_face[i_face] == i_face) {
-        descriptor.face_to_node_vector[l].resize(face_ending[i_face + 1] - face_ending[i_face]);
-
-        for (size_t i_face_node = 0; i_face_node < face_ending[i_face + 1] - face_ending[i_face]; ++i_face_node) {
-          descriptor.face_to_node_vector[l][i_face_node] = faces_node_list[face_ending[i_face] + i_face_node];
-        }
-        l++;
+    size_t i_face  = 0;
+    face_ending[0] = 0;
+    for (size_t i_dup_face = 0; i_dup_face < dup_face_to_face.size(); ++i_dup_face) {
+      if (dup_face_to_face[i_dup_face] == i_dup_face) {
+        face_ending[i_face + 1] =
+          dup_face_to_node_row[i_dup_face + 1] - dup_face_to_node_row[i_dup_face] + face_ending[i_face];
+        ++i_face;
       }
     }
   }
 
+  Array<unsigned int> faces_node_list(face_ending[face_ending.size() - 1]);
   {
-    descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
-    size_t l = 0;
-    for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
-      descriptor.cell_face_is_reversed_vector[j].resize(cell_nb_faces[j]);
-      for (size_t i_face = 0; i_face < cell_nb_faces[j]; ++i_face) {
-        descriptor.cell_face_is_reversed_vector[j][i_face] = face_reversed[l];
-        ++l;
+    size_t i_face  = 0;
+    face_ending[0] = 0;
+    for (size_t i_dup_face = 0; i_dup_face < dup_face_to_face.size(); ++i_dup_face) {
+      if (dup_face_to_face[i_dup_face] == i_dup_face) {
+        size_t dup_face_node_begin = dup_face_to_node_row[i_dup_face];
+        size_t face_node_begin     = face_ending[i_face];
+
+        for (size_t i_face_node = 0;
+             i_face_node < dup_face_to_node_row[i_dup_face + 1] - dup_face_to_node_row[i_dup_face]; ++i_face_node) {
+          faces_node_list[face_node_begin + i_face_node] = dup_faces_to_node_list[dup_face_node_begin + i_face_node];
+        }
+        ++i_face;
       }
     }
   }
+
+  descriptor.face_to_node_matrix = ConnectivityMatrix(face_ending, faces_node_list);
+
+  descriptor.cell_face_is_reversed = cell_face_is_reversed;
 }
 
 template <size_t Dimension>
@@ -557,45 +601,31 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
 {
   static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
 
-  Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
-  {
-    for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-      const auto& face_nodes = descriptor.face_to_node_vector[l];
+  Array<const unsigned int> face_to_edge_row = descriptor.face_to_node_matrix.rowsMap();
 
-      face_nb_edges[l] = face_nodes.size();
-    }
-  }
-
-  Array<unsigned int> face_edge_ending(face_nb_edges.size() + 1);
-  face_edge_ending[0] = 0;
-  for (size_t i_face = 0; i_face < face_nb_edges.size(); ++i_face) {
-    face_edge_ending[i_face + 1] = face_edge_ending[i_face] + face_nb_edges[i_face];
-  }
-
-  const size_t total_number_of_face_edges = face_edge_ending[face_edge_ending.size() - 1];
+  const size_t total_number_of_face_edges = face_to_edge_row[face_to_edge_row.size() - 1];
 
   const size_t total_number_of_node_by_face_edges = 2 * total_number_of_face_edges;
 
-  Array<unsigned int> face_edges_node_list(total_number_of_node_by_face_edges);
+  Array<unsigned int> face_edge_to_node_list(total_number_of_node_by_face_edges);
   {
     size_t i_edge_node = 0;
-    for (size_t i_face = 0; i_face < face_nb_edges.size(); ++i_face) {
-      for (size_t i_node = 0; i_node < descriptor.face_to_node_vector[i_face].size() - 1; ++i_node) {
-        face_edges_node_list[i_edge_node++] = descriptor.face_to_node_vector[i_face][i_node];
-        face_edges_node_list[i_edge_node++] = descriptor.face_to_node_vector[i_face][i_node + 1];
-      }
-      face_edges_node_list[i_edge_node++] =
-        descriptor.face_to_node_vector[i_face][descriptor.face_to_node_vector[i_face].size() - 1];
-      face_edges_node_list[i_edge_node++] = descriptor.face_to_node_vector[i_face][0];
+    for (size_t i_face = 0; i_face < face_to_edge_row.size() - 1; ++i_face) {
+      for (size_t i_node = 0; i_node < descriptor.face_to_node_matrix[i_face].size() - 1; ++i_node) {
+        face_edge_to_node_list[i_edge_node++] = descriptor.face_to_node_matrix[i_face][i_node];
+        face_edge_to_node_list[i_edge_node++] = descriptor.face_to_node_matrix[i_face][i_node + 1];
+      }
+      face_edge_to_node_list[i_edge_node++] =
+        descriptor.face_to_node_matrix[i_face][descriptor.face_to_node_matrix[i_face].size() - 1];
+      face_edge_to_node_list[i_edge_node++] = descriptor.face_to_node_matrix[i_face][0];
     }
   }
 
   Array<bool> face_edge_is_reversed(total_number_of_face_edges);
-
   for (size_t i_edge = 0; i_edge < total_number_of_face_edges; ++i_edge) {
-    if (descriptor.node_number_vector[face_edges_node_list[2 * i_edge]] >
-        descriptor.node_number_vector[face_edges_node_list[2 * i_edge + 1]]) {
-      std::swap(face_edges_node_list[2 * i_edge], face_edges_node_list[2 * i_edge + 1]);
+    if (descriptor.node_number_vector[face_edge_to_node_list[2 * i_edge]] >
+        descriptor.node_number_vector[face_edge_to_node_list[2 * i_edge + 1]]) {
+      std::swap(face_edge_to_node_list[2 * i_edge], face_edge_to_node_list[2 * i_edge + 1]);
       face_edge_is_reversed[i_edge] = true;
     } else {
       face_edge_is_reversed[i_edge] = false;
@@ -607,7 +637,7 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
   node_to_face_edge_row.fill(0);
   for (size_t i_edge = 0; i_edge < total_number_of_edges; ++i_edge) {
     for (size_t i_edge_node = 0; i_edge_node < 2; ++i_edge_node) {
-      node_to_face_edge_row[face_edges_node_list[2 * i_edge + i_edge_node] + 1]++;
+      node_to_face_edge_row[face_edge_to_node_list[2 * i_edge + i_edge_node] + 1]++;
     }
   }
   for (size_t i_node = 1; i_node < node_to_face_edge_row.size(); ++i_node) {
@@ -622,7 +652,7 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
 
     for (size_t i_edge = 0; i_edge < total_number_of_edges; ++i_edge) {
       for (size_t i_edge_node = 0; i_edge_node < 2; ++i_edge_node) {
-        const size_t node_id = face_edges_node_list[2 * i_edge + i_edge_node];
+        const size_t node_id = face_edge_to_node_list[2 * i_edge + i_edge_node];
 
         node_to_face_edge_list[node_to_face_edge_row[node_id] + node_to_face_edge_row_idx[node_id]] = i_edge;
         node_to_face_edge_row_idx[node_id]++;
@@ -637,8 +667,8 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
   auto is_same_face = [=](const size_t i_edge, const size_t j_edge) {
     auto i_edge_first_node = 2 * i_edge;
     auto j_edge_first_node = 2 * j_edge;
-    return ((face_edges_node_list[i_edge_first_node] == face_edges_node_list[j_edge_first_node]) and
-            (face_edges_node_list[i_edge_first_node + 1] == face_edges_node_list[j_edge_first_node + 1]));
+    return ((face_edge_to_node_list[i_edge_first_node] == face_edge_to_node_list[j_edge_first_node]) and
+            (face_edge_to_node_list[i_edge_first_node + 1] == face_edge_to_node_list[j_edge_first_node + 1]));
   };
 
   const auto& node_number_vector = descriptor.node_number_vector;
@@ -646,14 +676,16 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
   for (size_t i_node = 0; i_node < node_number_vector.size(); ++i_node) {
     for (size_t i_node_face = node_to_face_edge_row[i_node]; i_node_face < node_to_face_edge_row[i_node + 1] - 1;
          ++i_node_face) {
+      const unsigned int i_edge    = node_to_face_edge_list[i_node_face];
+      const unsigned int i_edge_id = face_edge_to_edge[i_edge];
       for (size_t j_node_face = i_node_face + 1; j_node_face < node_to_face_edge_row[i_node + 1]; ++j_node_face) {
-        if (face_edge_to_edge[node_to_face_edge_list[i_node_face]] ==
-            face_edge_to_edge[node_to_face_edge_list[j_node_face]])
-          continue;
-        if (is_same_face(node_to_face_edge_list[i_node_face], node_to_face_edge_list[j_node_face])) {
-          face_edge_to_edge[node_to_face_edge_list[j_node_face]] =
-            face_edge_to_edge[node_to_face_edge_list[i_node_face]];
-          nb_duplicate_edges++;
+        const unsigned int j_edge = node_to_face_edge_list[j_node_face];
+        unsigned int& j_edge_id   = face_edge_to_edge[j_edge];
+        if (i_edge_id != j_edge_id) {
+          if (is_same_face(i_edge, j_edge)) {
+            j_edge_id = i_edge_id;
+            nb_duplicate_edges++;
+          }
         }
       }
     }
@@ -677,13 +709,13 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
     }
   }
 
-  Array<NodeId> edge_to_node_list(2 * (total_number_of_edges - nb_duplicate_edges));
+  Array<unsigned int> edge_to_node_list(2 * (total_number_of_edges - nb_duplicate_edges));
   {
     for (size_t i_dup_edge = 0; i_dup_edge < total_number_of_edges; ++i_dup_edge) {
       if (face_edge_to_edge[i_dup_edge] == i_dup_edge) {
         const EdgeId edge_id               = dup_edge_to_edge_id[i_dup_edge];
-        edge_to_node_list[2 * edge_id]     = face_edges_node_list[2 * i_dup_edge];
-        edge_to_node_list[2 * edge_id + 1] = face_edges_node_list[2 * i_dup_edge + 1];
+        edge_to_node_list[2 * edge_id]     = face_edge_to_node_list[2 * i_dup_edge];
+        edge_to_node_list[2 * edge_id + 1] = face_edge_to_node_list[2 * i_dup_edge + 1];
       }
     }
   }
@@ -699,46 +731,19 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
     }();
   }
 
-  descriptor.edge_to_node_vector.resize(total_number_of_edges - nb_duplicate_edges);
-  {
-    for (EdgeId edge_id = 0; edge_id < total_number_of_edges - nb_duplicate_edges; ++edge_id) {
-      descriptor.edge_to_node_vector[edge_id].resize(2);
+  Array<unsigned int> edge_to_node_row(total_number_of_edges - nb_duplicate_edges + 1);
+  parallel_for(
+    edge_to_node_row.size(), PUGS_LAMBDA(const size_t i_edge) { edge_to_node_row[i_edge] = 2 * i_edge; });
 
-      for (size_t i_edge_node = 0; i_edge_node < 2; ++i_edge_node) {
-        descriptor.edge_to_node_vector[edge_id][i_edge_node] = edge_to_node_list[2 * edge_id + i_edge_node];
-      }
-    }
-  }
+  descriptor.edge_to_node_matrix = ConnectivityMatrix(edge_to_node_row, edge_to_node_list);
 
   // Use real edge ids
   for (size_t i_edge = 0; i_edge < face_edge_to_edge.size(); ++i_edge) {
     face_edge_to_edge[i_edge] = dup_edge_to_edge_id[face_edge_to_edge[i_edge]];
   }
 
-  descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
-  {
-    size_t i_edge = 0;
-    for (FaceId face_id = 0; face_id < descriptor.face_to_edge_vector.size(); ++face_id) {
-      descriptor.face_to_edge_vector[face_id].resize(descriptor.face_to_node_vector[face_id].size());
-
-      for (size_t i_face_edge = 0; i_face_edge < descriptor.face_to_edge_vector[face_id].size(); ++i_face_edge) {
-        descriptor.face_to_edge_vector[face_id][i_face_edge] = face_edge_to_edge[i_edge++];
-      }
-    }
-  }
-
-  descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
-  {
-    size_t i_edge = 0;
-    for (FaceId face_id = 0; face_id < descriptor.face_edge_is_reversed_vector.size(); ++face_id) {
-      descriptor.face_edge_is_reversed_vector[face_id].resize(descriptor.face_to_node_vector[face_id].size());
-
-      for (size_t i_face_edge = 0; i_face_edge < descriptor.face_edge_is_reversed_vector[face_id].size();
-           ++i_face_edge) {
-        descriptor.face_edge_is_reversed_vector[face_id][i_face_edge] = face_edge_is_reversed[i_edge++];
-      }
-    }
-  }
+  descriptor.face_to_edge_matrix   = ConnectivityMatrix(face_to_edge_row, face_edge_to_edge);
+  descriptor.face_edge_is_reversed = face_edge_is_reversed;
 
   Array<size_t> node_to_duplicated_edge_id_list(node_to_face_edge_list.size());
   for (size_t i_node_edge = 0; i_node_edge < node_to_duplicated_edge_id_list.size(); ++i_node_edge) {
@@ -772,7 +777,7 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
     node_to_edge_row[node_id + 1] = node_to_edge_row[node_id] + node_nb_edges[node_id];
   }
 
-  Array<EdgeId> node_to_edge_list(node_to_edge_row[node_to_edge_row.size() - 1]);
+  Array<unsigned int> node_to_edge_list(node_to_edge_row[node_to_edge_row.size() - 1]);
   {
     size_t l = 0;
     for (size_t i_edge_id = 0; i_edge_id < node_to_duplicated_edge_id_list.size(); ++i_edge_id) {
@@ -782,6 +787,8 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
     }
   }
 
+  descriptor.node_to_edge_matrix = ConnectivityMatrix(node_to_edge_row, node_to_edge_list);
+
   auto find_edge = [=](const NodeId& node0_id, const NodeId& node1_id) -> EdgeId {
     auto [first_node_id, second_node_id] = [&] {
       if (descriptor.node_number_vector[node0_id] < descriptor.node_number_vector[node1_id]) {
@@ -802,108 +809,124 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
   };
 
   {
-    descriptor.cell_to_edge_vector.reserve(descriptor.cell_to_node_vector.size());
-    for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
-      const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+    Array<unsigned int> cell_to_edge_row(descriptor.cell_to_node_matrix.numberOfRows() + 1);
+    {
+      cell_to_edge_row[0] = 0;
 
-      switch (descriptor.cell_type_vector[j]) {
-      case CellType::Tetrahedron: {
-        constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(6);
-        for (int i_edge = 0; i_edge < 6; ++i_edge) {
-          const auto& e        = local_edge[i_edge];
-          const EdgeId edge_id = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
-          cell_edge_vector.push_back(edge_id);
-        }
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
-      }
-      case CellType::Hexahedron: {
-        constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
-                                           {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(12);
-        for (int i_edge = 0; i_edge < 12; ++i_edge) {
-          const auto& e        = local_edge[i_edge];
-          const EdgeId edge_id = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
-          cell_edge_vector.push_back(edge_id);
-        }
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
-      }
-      case CellType::Prism: {
-        constexpr int local_edge[9][2] = {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(9);
-        for (int i_edge = 0; i_edge < 9; ++i_edge) {
-          const auto& e        = local_edge[i_edge];
-          const EdgeId edge_id = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
-          cell_edge_vector.push_back(edge_id);
-        }
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
+      for (CellId cell_id = 0; cell_id < descriptor.cell_to_node_matrix.numberOfRows(); ++cell_id) {
+        switch (descriptor.cell_type_vector[cell_id]) {
+        case CellType::Tetrahedron: {
+          cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 6;
+          break;
+        }
+        case CellType::Hexahedron: {
+          cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 12;
+          break;
+        }
+        case CellType::Prism: {
+          cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 9;
+          break;
+        }
+        case CellType::Pyramid: {
+          cell_to_edge_row[cell_id + 1] =
+            cell_to_edge_row[cell_id] + 2 * (descriptor.cell_to_node_matrix[cell_id].size() - 1);
+          break;
+        }
+        case CellType::Diamond: {
+          cell_to_edge_row[cell_id + 1] =
+            cell_to_edge_row[cell_id] + 3 * (descriptor.cell_to_node_matrix[cell_id].size() - 2);
+          break;
+        }
+        default: {
+          std::stringstream error_msg;
+          error_msg << name(descriptor.cell_type_vector[cell_id]) << ": unexpected cell type in dimension 3";
+          throw NotImplementedError(error_msg.str());
+        }
+        }
       }
-      case CellType::Pyramid: {
-        const size_t number_of_edges = 2 * cell_nodes.size();
-        std::vector<unsigned int> base_nodes;
-        std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
+    }
 
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(number_of_edges);
+    Array<unsigned int> cell_to_edge_list(cell_to_edge_row[cell_to_edge_row.size() - 1]);
+    {
+      size_t i_cell_edge = 0;
+      for (CellId cell_id = 0; cell_id < descriptor.cell_to_node_matrix.numberOfRows(); ++cell_id) {
+        const auto& cell_nodes = descriptor.cell_to_node_matrix[cell_id];
 
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          const EdgeId edge_id = find_edge(base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]);
-          cell_edge_vector.push_back(edge_id);
+        switch (descriptor.cell_type_vector[cell_id]) {
+        case CellType::Tetrahedron: {
+          constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
+          for (int i_edge = 0; i_edge < 6; ++i_edge) {
+            const auto& e                    = local_edge[i_edge];
+            cell_to_edge_list[i_cell_edge++] = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
+          }
+          break;
         }
+        case CellType::Hexahedron: {
+          constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
+                                             {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
 
-        const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          const EdgeId edge_id = find_edge(base_nodes[i_edge], top_vertex);
-          cell_edge_vector.push_back(edge_id);
+          for (int i_edge = 0; i_edge < 12; ++i_edge) {
+            const auto& e                    = local_edge[i_edge];
+            cell_to_edge_list[i_cell_edge++] = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
+          }
+          break;
         }
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
-      }
-      case CellType::Diamond: {
-        const size_t number_of_edges = 3 * cell_nodes.size();
-        std::vector<unsigned int> base_nodes;
-        std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
-
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(number_of_edges);
-
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          const EdgeId edge_id = find_edge(base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]);
-          cell_edge_vector.push_back(edge_id);
+        case CellType::Prism: {
+          constexpr int local_edge[9][2] = {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
+          for (int i_edge = 0; i_edge < 9; ++i_edge) {
+            const auto& e                    = local_edge[i_edge];
+            cell_to_edge_list[i_cell_edge++] = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
+          }
+          break;
         }
+        case CellType::Pyramid: {
+          auto base_nodes = [&](size_t i) { return cell_nodes[i]; };
+
+          for (size_t i_edge = 0; i_edge < cell_nodes.size() - 1; ++i_edge) {
+            cell_to_edge_list[i_cell_edge++] =
+              find_edge(base_nodes(i_edge), base_nodes((i_edge + 1) % (cell_nodes.size() - 1)));
+          }
 
-        {
           const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
-          for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-            const EdgeId edge_id = find_edge(base_nodes[i_edge], top_vertex);
-            cell_edge_vector.push_back(edge_id);
+          for (size_t i_edge = 0; i_edge < cell_nodes.size() - 1; ++i_edge) {
+            cell_to_edge_list[i_cell_edge++] = find_edge(base_nodes(i_edge), top_vertex);
           }
+          break;
         }
+        case CellType::Diamond: {
+          auto base_nodes = [&](size_t i) { return cell_nodes[i + 1]; };
 
-        {
-          const unsigned int bottom_vertex = cell_nodes[0];
-          for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-            const EdgeId edge_id = find_edge(base_nodes[i_edge], bottom_vertex);
-            cell_edge_vector.push_back(edge_id);
+          for (size_t i_edge = 0; i_edge < cell_nodes.size() - 2; ++i_edge) {
+            cell_to_edge_list[i_cell_edge++] =
+              find_edge(base_nodes(i_edge), base_nodes((i_edge + 1) % (cell_nodes.size() - 2)));
           }
-        }
 
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
-      }
-      default: {
-        std::stringstream error_msg;
-        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
-        throw NotImplementedError(error_msg.str());
-      }
+          {
+            const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+            for (size_t i_edge = 0; i_edge < cell_nodes.size() - 2; ++i_edge) {
+              cell_to_edge_list[i_cell_edge++] = find_edge(base_nodes(i_edge), top_vertex);
+            }
+          }
+
+          {
+            const unsigned int bottom_vertex = cell_nodes[0];
+            for (size_t i_edge = 0; i_edge < cell_nodes.size() - 2; ++i_edge) {
+              cell_to_edge_list[i_cell_edge++] = find_edge(base_nodes(i_edge), bottom_vertex);
+            }
+          }
+
+          break;
+        }
+        default: {
+          std::stringstream error_msg;
+          error_msg << name(descriptor.cell_type_vector[cell_id]) << ": unexpected cell type in dimension 3";
+          throw NotImplementedError(error_msg.str());
+        }
+        }
       }
     }
+
+    descriptor.cell_to_edge_matrix = ConnectivityMatrix(cell_to_edge_row, cell_to_edge_list);
   }
 }
 
diff --git a/src/mesh/ConnectivityBuilderBase.hpp b/src/mesh/ConnectivityBuilderBase.hpp
index 9ae5df81d4792f978747877fb9c904faa1f6de07..a621e91f46ffe7e6ecb5ea6d496934a03571ba6e 100644
--- a/src/mesh/ConnectivityBuilderBase.hpp
+++ b/src/mesh/ConnectivityBuilderBase.hpp
@@ -15,9 +15,6 @@ class ConnectivityDescriptor;
 class ConnectivityBuilderBase
 {
  protected:
-  template <size_t Dimension>
-  class ConnectivityFace;
-
   std::shared_ptr<const IConnectivity> m_connectivity;
 
   template <size_t Dimension>
@@ -37,190 +34,4 @@ class ConnectivityBuilderBase
   ~ConnectivityBuilderBase() = default;
 };
 
-template <>
-class ConnectivityBuilderBase::ConnectivityFace<2>
-{
- public:
-  friend struct Hash;
-
-  struct Hash
-  {
-    size_t
-    operator()(const ConnectivityFace& f) const
-    {
-      size_t hash = 0;
-      hash ^= std::hash<unsigned int>()(f.m_node0_id);
-      hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
-      return hash;
-    }
-  };
-
- private:
-  const Array<int>& m_node_number_vector;
-
-  unsigned int m_node0_id;
-  unsigned int m_node1_id;
-
-  bool m_reversed;
-
- public:
-  std::vector<unsigned int>
-  nodeIdList() const
-  {
-    return {m_node0_id, m_node1_id};
-  }
-
-  bool
-  reversed() const
-  {
-    return m_reversed;
-  }
-
-  PUGS_INLINE
-  bool
-  operator==(const ConnectivityFace& f) const
-  {
-    return ((m_node0_id == f.m_node0_id) and (m_node1_id == f.m_node1_id));
-  }
-
-  PUGS_INLINE
-  bool
-  operator<(const ConnectivityFace& f) const
-  {
-    return ((m_node_number_vector[m_node0_id] < m_node_number_vector[f.m_node0_id]) or
-            ((m_node_number_vector[m_node0_id] == m_node_number_vector[f.m_node0_id]) and
-             (m_node_number_vector[m_node1_id] < m_node_number_vector[f.m_node1_id])));
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const std::vector<unsigned int>& node_id_list, const Array<int>& node_number_vector)
-    : m_node_number_vector(node_number_vector)
-  {
-    Assert(node_id_list.size() == 2);
-
-    if (m_node_number_vector[node_id_list[0]] < m_node_number_vector[node_id_list[1]]) {
-      m_node0_id = node_id_list[0];
-      m_node1_id = node_id_list[1];
-      m_reversed = false;
-    } else {
-      m_node0_id = node_id_list[1];
-      m_node1_id = node_id_list[0];
-      m_reversed = true;
-    }
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const ConnectivityFace&) = default;
-
-  PUGS_INLINE
-  ~ConnectivityFace() = default;
-};
-
-template <>
-class ConnectivityBuilderBase::ConnectivityFace<3>
-{
- public:
-  friend struct Hash;
-
-  struct Hash
-  {
-    size_t
-    operator()(const ConnectivityFace& f) const
-    {
-      size_t hash = 0;
-      for (size_t i = 0; i < f.m_node_id_list.size(); ++i) {
-        hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
-      }
-      return hash;
-    }
-  };
-
- private:
-  bool m_reversed;
-  std::vector<NodeId::base_type> m_node_id_list;
-  const Array<int>& m_node_number_vector;
-
-  PUGS_INLINE
-  std::vector<unsigned int>
-  _sort(const std::vector<unsigned int>& node_list)
-  {
-    const auto min_id = std::min_element(node_list.begin(), node_list.end());
-    const int shift   = std::distance(node_list.begin(), min_id);
-
-    std::vector<unsigned int> rotated_node_list(node_list.size());
-    if (node_list[(shift + 1) % node_list.size()] > node_list[(shift + node_list.size() - 1) % node_list.size()]) {
-      for (size_t i = 0; i < node_list.size(); ++i) {
-        rotated_node_list[i] = node_list[(shift + node_list.size() - i) % node_list.size()];
-        m_reversed           = true;
-      }
-    } else {
-      for (size_t i = 0; i < node_list.size(); ++i) {
-        rotated_node_list[i] = node_list[(shift + i) % node_list.size()];
-      }
-    }
-
-    return rotated_node_list;
-  }
-
- public:
-  PUGS_INLINE
-  const bool&
-  reversed() const
-  {
-    return m_reversed;
-  }
-
-  PUGS_INLINE
-  const std::vector<unsigned int>&
-  nodeIdList() const
-  {
-    return m_node_id_list;
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list, const Array<int>& node_number_vector)
-    : m_reversed(false), m_node_id_list(_sort(given_node_id_list)), m_node_number_vector(node_number_vector)
-  {
-    ;
-  }
-
- public:
-  bool
-  operator==(const ConnectivityFace& f) const
-  {
-    if (m_node_id_list.size() == f.nodeIdList().size()) {
-      for (size_t j = 0; j < m_node_id_list.size(); ++j) {
-        if (m_node_id_list[j] != f.nodeIdList()[j]) {
-          return false;
-        }
-      }
-      return true;
-    }
-    return false;
-  }
-
-  PUGS_INLINE
-  bool
-  operator<(const ConnectivityFace& f) const
-  {
-    const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
-    for (size_t i = 0; i < min_nb_nodes; ++i) {
-      if (m_node_id_list[i] < f.m_node_id_list[i])
-        return true;
-      if (m_node_id_list[i] != f.m_node_id_list[i])
-        return false;
-    }
-    return m_node_id_list.size() < f.m_node_id_list.size();
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const ConnectivityFace&) = default;
-
-  PUGS_INLINE
-  ConnectivityFace() = delete;
-
-  PUGS_INLINE
-  ~ConnectivityFace() = default;
-};
-
 #endif   // CONNECTIVITY_BUILDER_BASE_HPP
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index aa8ca779cfb965ffb3e85e4feb1b4b9feabb79ca..1d5514627a5dac75d41ba5e4d08fa182d37e552f 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -7,62 +7,90 @@
 
 template <typename ConnectivityType>
 PUGS_INLINE ConnectivityMatrix
-ConnectivityComputer::computeConnectivityMatrix(const ConnectivityType& connectivity,
-                                                ItemType item_type,
-                                                ItemType child_item_type) const
+ConnectivityComputer::computeInverseConnectivityMatrix(const ConnectivityType& connectivity,
+                                                       ItemType item_type,
+                                                       ItemType child_item_type) const
 {
-  ConnectivityMatrix item_to_child_item_matrix;
-  if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
-    const ConnectivityMatrix& child_to_item_matrix = connectivity.getMatrix(child_item_type, item_type);
+  if (item_type < child_item_type) {
+    ConnectivityMatrix item_to_child_item_matrix;
+    if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
+      const ConnectivityMatrix& child_to_item_matrix = connectivity.getMatrix(child_item_type, item_type);
+
+      switch (child_item_type) {
+      case ItemType::cell: {
+        item_to_child_item_matrix =
+          this->_computeInverseConnectivity<ConnectivityType, ItemType::cell>(connectivity, child_to_item_matrix);
+        break;
+      }
+      case ItemType::face: {
+        item_to_child_item_matrix =
+          this->_computeInverseConnectivity<ConnectivityType, ItemType::face>(connectivity, child_to_item_matrix);
+        break;
+      }
+      case ItemType::edge: {
+        item_to_child_item_matrix =
+          this->_computeInverseConnectivity<ConnectivityType, ItemType::edge>(connectivity, child_to_item_matrix);
+        break;
+      }
+        // LCOV_EXCL_START
+      default: {
+        throw UnexpectedError("node cannot be child item when computing inverse connectivity");
+      }
+        // LCOV_EXCL_STOP
+      }
+    } else {
+      std::stringstream error_msg;
+      error_msg << "unable to compute connectivity " << itemName(item_type) << " -> " << itemName(child_item_type);
+      throw UnexpectedError(error_msg.str());
+    }
 
-    item_to_child_item_matrix = this->_computeInverse(child_to_item_matrix);
+    return item_to_child_item_matrix;
   } else {
     std::stringstream error_msg;
-    error_msg << "unable to compute connectivity " << itemName(item_type) << " -> " << itemName(child_item_type);
+    error_msg << "cannot deduce " << itemName(item_type) << " -> " << itemName(child_item_type) << " connectivity";
     throw UnexpectedError(error_msg.str());
   }
-
-  return item_to_child_item_matrix;
 }
 
-template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(const Connectivity1D&,
-                                                                            ItemType,
-                                                                            ItemType) const;
+template ConnectivityMatrix ConnectivityComputer::computeInverseConnectivityMatrix(const Connectivity1D&,
+                                                                                   ItemType,
+                                                                                   ItemType) const;
 
-template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(const Connectivity2D&,
-                                                                            ItemType,
-                                                                            ItemType) const;
+template ConnectivityMatrix ConnectivityComputer::computeInverseConnectivityMatrix(const Connectivity2D&,
+                                                                                   ItemType,
+                                                                                   ItemType) const;
 
-template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(const Connectivity3D&,
-                                                                            ItemType,
-                                                                            ItemType) const;
+template ConnectivityMatrix ConnectivityComputer::computeInverseConnectivityMatrix(const Connectivity3D&,
+                                                                                   ItemType,
+                                                                                   ItemType) const;
 
+template <typename ConnectivityType, ItemType child_item_type>
 ConnectivityMatrix
-ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
+ConnectivityComputer::_computeInverseConnectivity(const ConnectivityType& connectivity,
+                                                  const ConnectivityMatrix& child_to_item_matrix) const
 {
-  const size_t& number_of_rows = item_to_child_matrix.numberOfRows();
+  const size_t number_of_transposed_columns = child_to_item_matrix.numberOfRows();
 
-  if ((item_to_child_matrix.values().size() > 0)) {
-    const size_t& number_of_columns = max(item_to_child_matrix.values());
+  if ((child_to_item_matrix.values().size() > 0)) {
+    const size_t number_of_transposed_rows = max(child_to_item_matrix.values()) + 1;
 
-    Array<uint32_t> transposed_next_free_column_index(number_of_columns + 1);
+    Array<uint32_t> transposed_next_free_column_index(number_of_transposed_rows);
     transposed_next_free_column_index.fill(0);
 
     Array<uint32_t> transposed_rows_map(transposed_next_free_column_index.size() + 1);
     transposed_rows_map.fill(0);
-    for (size_t i = 0; i < number_of_rows; ++i) {
-      for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
-        transposed_rows_map[item_to_child_matrix.values()[j] + 1]++;
+    for (size_t i = 0; i < number_of_transposed_columns; ++i) {
+      for (size_t j = child_to_item_matrix.rowsMap()[i]; j < child_to_item_matrix.rowsMap()[i + 1]; ++j) {
+        transposed_rows_map[child_to_item_matrix.values()[j] + 1]++;
       }
     }
     for (size_t i = 1; i < transposed_rows_map.size(); ++i) {
       transposed_rows_map[i] += transposed_rows_map[i - 1];
     }
     Array<uint32_t> transposed_column_indices(transposed_rows_map[transposed_rows_map.size() - 1]);
-
-    for (size_t i = 0; i < number_of_rows; ++i) {
-      for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
-        size_t i_column_index = item_to_child_matrix.values()[j];
+    for (size_t i = 0; i < number_of_transposed_columns; ++i) {
+      for (size_t j = child_to_item_matrix.rowsMap()[i]; j < child_to_item_matrix.rowsMap()[i + 1]; ++j) {
+        size_t i_column_index = child_to_item_matrix.values()[j];
         uint32_t& shift       = transposed_next_free_column_index[i_column_index];
 
         transposed_column_indices[transposed_rows_map[i_column_index] + shift] = i;
@@ -71,6 +99,18 @@ ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_ma
       }
     }
 
+    auto target_item_number = connectivity.template number<child_item_type>();
+    // Finally one sorts target item_ids for parallel reproducibility
+    for (size_t i = 0; i < number_of_transposed_rows; ++i) {
+      auto row_begining = &(transposed_column_indices[transposed_rows_map[i]]);
+      auto row_size     = (transposed_rows_map[i + 1] - transposed_rows_map[i]);
+      std::sort(row_begining, row_begining + row_size,
+                [&target_item_number](const ItemIdT<child_item_type> item0_id,
+                                      const ItemIdT<child_item_type> item1_id) {
+                  return target_item_number[item0_id] < target_item_number[item1_id];
+                });
+    }
+
     return ConnectivityMatrix{transposed_rows_map, transposed_column_indices};
   } else {
     // empty connectivity
diff --git a/src/mesh/ConnectivityComputer.hpp b/src/mesh/ConnectivityComputer.hpp
index b033a25d5ccc04db84c194e79e492a081216a49e..7093aeef17707c37452b5e6874f6ed50a285dbb3 100644
--- a/src/mesh/ConnectivityComputer.hpp
+++ b/src/mesh/ConnectivityComputer.hpp
@@ -7,13 +7,15 @@
 class ConnectivityComputer
 {
  private:
-  ConnectivityMatrix _computeInverse(const ConnectivityMatrix& item_to_child_matrix) const;
+  template <typename ConnectivityType, ItemType child_item_type>
+  ConnectivityMatrix _computeInverseConnectivity(const ConnectivityType& connectivity,
+                                                 const ConnectivityMatrix& item_to_child_matrix) const;
 
  public:
   template <typename ConnectivityType>
-  ConnectivityMatrix computeConnectivityMatrix(const ConnectivityType& connectivity,
-                                               ItemType item_type,
-                                               ItemType child_item_type) const;
+  ConnectivityMatrix computeInverseConnectivityMatrix(const ConnectivityType& connectivity,
+                                                      ItemType item_type,
+                                                      ItemType child_item_type) const;
 
   template <typename ItemOfItem, typename ConnectivityType>
   void computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const;
diff --git a/src/mesh/ConnectivityDescriptor.hpp b/src/mesh/ConnectivityDescriptor.hpp
index 6e8733895aa6c942af199a96f65780d5533e87ca..536345996ffa9a3b9fdf84e0f047074af6e1c798 100644
--- a/src/mesh/ConnectivityDescriptor.hpp
+++ b/src/mesh/ConnectivityDescriptor.hpp
@@ -2,6 +2,7 @@
 #define CONNECTIVITY_DESCRIPTOR_HPP
 
 #include <mesh/CellType.hpp>
+#include <mesh/ConnectivityMatrix.hpp>
 #include <mesh/ItemOfItemType.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/PugsTraits.hpp>
@@ -17,38 +18,41 @@ class ConnectivityDescriptor
   std::vector<RefNodeList> m_ref_node_list_vector;
 
  public:
-  std::vector<std::vector<unsigned int>> cell_to_node_vector;
-  std::vector<std::vector<unsigned int>> cell_to_face_vector;
-  std::vector<std::vector<unsigned int>> cell_to_edge_vector;
+  ConnectivityMatrix cell_to_face_matrix = ConnectivityMatrix{true};
+  ConnectivityMatrix cell_to_edge_matrix = ConnectivityMatrix{true};
+  ConnectivityMatrix cell_to_node_matrix = ConnectivityMatrix{true};
 
-  std::vector<std::vector<unsigned int>> face_to_node_vector;
-  std::vector<std::vector<unsigned int>> face_to_edge_vector;
+  ConnectivityMatrix face_to_edge_matrix = ConnectivityMatrix{true};
+  ConnectivityMatrix face_to_node_matrix = ConnectivityMatrix{true};
 
-  std::vector<std::vector<unsigned int>> edge_to_node_vector;
+  ConnectivityMatrix edge_to_node_matrix = ConnectivityMatrix{true};
+
+  ConnectivityMatrix node_to_face_matrix = ConnectivityMatrix{true};
+  ConnectivityMatrix node_to_edge_matrix = ConnectivityMatrix{true};
 
   template <typename ItemOfItemT>
   auto&
   itemOfItemVector()
   {
     if constexpr (std::is_same_v<ItemOfItemT, NodeOfCell>) {
-      return cell_to_node_vector;
+      return cell_to_node_matrix;
     } else if constexpr (std::is_same_v<ItemOfItemT, FaceOfCell>) {
-      return cell_to_face_vector;
+      return cell_to_face_matrix;
     } else if constexpr (std::is_same_v<ItemOfItemT, EdgeOfCell>) {
-      return cell_to_edge_vector;
+      return cell_to_edge_matrix;
     } else if constexpr (std::is_same_v<ItemOfItemT, EdgeOfFace>) {
-      return face_to_edge_vector;
+      return face_to_edge_matrix;
     } else if constexpr (std::is_same_v<ItemOfItemT, NodeOfFace>) {
-      return face_to_node_vector;
+      return face_to_node_matrix;
     } else if constexpr (std::is_same_v<ItemOfItemT, NodeOfEdge>) {
-      return edge_to_node_vector;
+      return edge_to_node_matrix;
     } else {
       static_assert(is_false_v<ItemOfItemT>, "Unexpected item of item type");
     }
   }
 
-  std::vector<std::vector<bool>> cell_face_is_reversed_vector;
-  std::vector<std::vector<bool>> face_edge_is_reversed_vector;
+  Array<bool> cell_face_is_reversed;
+  Array<bool> face_edge_is_reversed;
 
   Array<CellType> cell_type_vector;
 
@@ -114,7 +118,7 @@ class ConnectivityDescriptor
   }
 
   ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
-  ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
+  ConnectivityDescriptor& operator=(ConnectivityDescriptor&&)      = delete;
 
   ConnectivityDescriptor()                              = default;
   ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index 94ece0650d26333e10037d19368bb84d4850abc0..552c4ab959a18c0d8d10b281ca83ee258bdd87ce 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -145,19 +145,19 @@ template <int Dimension>
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
 ConnectivityDispatcher<Dimension>::_gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
-                                               Array<std::remove_const_t<DataType>>& gathered_vector)
+                                               Array<std::remove_const_t<DataType>>& gathered_array)
 {
   std::vector<Array<const DataType>> recv_item_data_by_proc = this->exchange(data_to_gather);
 
   const auto& recv_id_correspondance_by_proc = this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
   Assert(recv_id_correspondance_by_proc.size() == parallel::size());
 
-  gathered_vector = Array<std::remove_const_t<DataType>>(this->_dispatchedInfo<item_type>().m_number_to_id_map.size());
+  gathered_array = Array<std::remove_const_t<DataType>>(this->_dispatchedInfo<item_type>().m_number_to_id_map.size());
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
     Assert(recv_id_correspondance_by_proc[i_rank].size() == recv_item_data_by_proc[i_rank].size());
     for (size_t r = 0; r < recv_id_correspondance_by_proc[i_rank].size(); ++r) {
-      const auto& item_id      = recv_id_correspondance_by_proc[i_rank][r];
-      gathered_vector[item_id] = recv_item_data_by_proc[i_rank][r];
+      const auto& item_id     = recv_id_correspondance_by_proc[i_rank][r];
+      gathered_array[item_id] = recv_item_data_by_proc[i_rank][r];
     }
   }
 }
@@ -167,7 +167,7 @@ template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
 void
 ConnectivityDispatcher<Dimension>::_gatherFrom(
   const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
-  std::vector<std::vector<std::remove_const_t<DataType>>>& gathered_vector)
+  Array<std::remove_const_t<DataType>>& gathered_array)
 {
   using MutableDataType = std::remove_const_t<DataType>;
 
@@ -201,17 +201,23 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(
 
   parallel::exchange(data_to_send_by_proc, recv_data_to_gather_by_proc);
 
-  const auto& item_list_to_recv_size_by_proc = this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
+  const size_t recv_array_size = [&] {
+    size_t size = 0;
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      size += recv_data_to_gather_by_proc[i_rank].size();
+    }
+    return size;
+  }();
 
-  for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-    int l = 0;
-    for (size_t i = 0; i < item_list_to_recv_size_by_proc[i_rank]; ++i) {
-      std::vector<MutableDataType> data_vector(number_of_sub_item_per_item_to_recv_by_proc[i_rank][i]);
-      for (size_t k = 0; k < data_vector.size(); ++k) {
-        data_vector[k] = recv_data_to_gather_by_proc[i_rank][l++];
+  gathered_array = Array<std::remove_const_t<DataType>>(recv_array_size);
+  {
+    size_t l = 0;
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      for (size_t j = 0; j < recv_data_to_gather_by_proc[i_rank].size(); ++j) {
+        gathered_array[l++] = recv_data_to_gather_by_proc[i_rank][j];
       }
-      gathered_vector.emplace_back(data_vector);
     }
+    Assert(gathered_array.size() == l);
   }
 }
 
@@ -342,6 +348,8 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
   const auto& recv_item_of_item_numbers_by_proc =
     this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
 
+  std::vector<std::vector<unsigned int>> item_to_subitem_legacy;
+  size_t number_of_node_by_cell = 0;
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
     int l = 0;
     for (size_t i = 0; i < item_list_to_recv_size_by_proc[i_rank]; ++i) {
@@ -351,9 +359,30 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
         Assert(searched_sub_item_id != sub_item_number_id_map.end());
         sub_item_vector.push_back(searched_sub_item_id->second);
       }
-      m_new_descriptor.itemOfItemVector<ItemOfItemT>().emplace_back(sub_item_vector);
+      number_of_node_by_cell += sub_item_vector.size();
+
+      item_to_subitem_legacy.emplace_back(sub_item_vector);
     }
   }
+  Array<unsigned int> item_to_subitem_row_map(item_to_subitem_legacy.size() + 1);
+  Array<unsigned int> item_to_subitem_list(number_of_node_by_cell);
+
+  item_to_subitem_row_map.fill(10000000);
+  item_to_subitem_list.fill(10000000);
+
+  item_to_subitem_row_map[0] = 0;
+  for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
+    item_to_subitem_row_map[i + 1] = item_to_subitem_row_map[i] + item_to_subitem_legacy[i].size();
+  }
+  size_t l = 0;
+  for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
+    const auto& subitem_list = item_to_subitem_legacy[i];
+    for (size_t j = 0; j < subitem_list.size(); ++j, ++l) {
+      item_to_subitem_list[l] = subitem_list[j];
+    }
+  }
+
+  m_new_descriptor.itemOfItemVector<ItemOfItemT>() = ConnectivityMatrix(item_to_subitem_row_map, item_to_subitem_list);
 }
 
 template <int Dimension>
@@ -594,7 +623,7 @@ ConnectivityDispatcher<Dimension>::_dispatchEdges()
     this->_buildSubItemNumbersToRecvByProc<EdgeOfFace>();
     this->_buildItemToSubItemDescriptor<EdgeOfFace>();
 
-    this->_gatherFrom(m_connectivity.faceEdgeIsReversed(), m_new_descriptor.face_edge_is_reversed_vector);
+    this->_gatherFrom(m_connectivity.faceEdgeIsReversed(), m_new_descriptor.face_edge_is_reversed);
 
     this->_gatherFrom(this->_dispatchedInfo<ItemType::edge>().m_new_owner, m_new_descriptor.edge_owner_vector);
 
@@ -620,7 +649,7 @@ ConnectivityDispatcher<Dimension>::_dispatchFaces()
 
     this->_buildItemToSubItemDescriptor<FaceOfCell>();
 
-    this->_gatherFrom(m_connectivity.cellFaceIsReversed(), m_new_descriptor.cell_face_is_reversed_vector);
+    this->_gatherFrom(m_connectivity.cellFaceIsReversed(), m_new_descriptor.cell_face_is_reversed);
 
     this->_gatherFrom(this->_dispatchedInfo<ItemType::face>().m_new_owner, m_new_descriptor.face_owner_vector);
 
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index 58b9e8b5100308453b9508646a164ed1b369fe43..9ff5913b60799148f99bb6fdf4cff2cedd56973c 100644
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -180,11 +180,11 @@ class ConnectivityDispatcher
 
   template <typename DataType, ItemType item_type, typename ConnectivityPtr>
   void _gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
-                   Array<std::remove_const_t<DataType>>& gathered_vector);
+                   Array<std::remove_const_t<DataType>>& gathered_array);
 
   template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
   void _gatherFrom(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
-                   std::vector<std::vector<std::remove_const_t<DataType>>>& gathered_vector);
+                   Array<std::remove_const_t<DataType>>& gathered_array);
 
   template <typename SubItemOfItemT>
   void _buildNumberOfSubItemPerItemToRecvByProc();
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 6b777288cd6b8e3793a9addff6629c661d45d82d..eadc85347cb756932616e655d7b625a2a15505b7 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -70,34 +70,22 @@ class ConnectivityMatrix
 #endif   // NDEBUG
   }
 
-  PUGS_INLINE
-  ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer) noexcept : m_is_built{true}
-  {
-    m_row_map = [&] {
-      Array<uint32_t> row_map(initializer.size() + 1);
-      row_map[0] = 0;
-      for (size_t i = 0; i < initializer.size(); ++i) {
-        row_map[i + 1] = row_map[i] + initializer[i].size();
-      }
-      return row_map;
-    }();
-
-    m_column_indices = [&] {
-      Array<uint32_t> column_indices(m_row_map[m_row_map.size() - 1]);
-      size_t index = 0;
-      for (const auto& row : initializer) {
-        for (const auto& col_index : row) {
-          column_indices[index++] = col_index;
-        }
-      }
-      return column_indices;
-    }();
-  }
-
   ConnectivityMatrix& operator=(const ConnectivityMatrix&) = default;
-  ConnectivityMatrix& operator=(ConnectivityMatrix&&) = default;
+  ConnectivityMatrix& operator=(ConnectivityMatrix&&)      = default;
 
-  ConnectivityMatrix()                          = default;
+  ConnectivityMatrix(bool is_built = false) : m_is_built{is_built}
+  {
+    // this is useful to build
+    if (is_built) {
+      m_row_map = [&] {
+        Array<uint32_t> row_map(1);
+        row_map[0] = 0;
+        return row_map;
+      }();
+
+      m_column_indices = Array<uint32_t>(0);
+    }
+  }
   ConnectivityMatrix(const ConnectivityMatrix&) = default;
   ConnectivityMatrix(ConnectivityMatrix&&)      = default;
   ~ConnectivityMatrix()                         = default;
diff --git a/src/mesh/ConnectivityUtils.cpp b/src/mesh/ConnectivityUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..af0220ac448b76fb2e0101c278719cc07da989bc
--- /dev/null
+++ b/src/mesh/ConnectivityUtils.cpp
@@ -0,0 +1,88 @@
+#include <mesh/ConnectivityUtils.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension, ItemType SourceItemType, ItemType TargetItemType>
+bool
+checkItemToHigherDimensionItem(const Connectivity<Dimension>& connectivity)
+{
+  static_assert(SourceItemType < TargetItemType);
+  bool is_valid = true;
+
+  const auto& item_to_item_matrix = connectivity.template getItemToItemMatrix<SourceItemType, TargetItemType>();
+  const auto& item_number         = connectivity.template number<TargetItemType>();
+
+  for (ItemIdT<SourceItemType> source_item_id = 0; source_item_id < connectivity.template numberOf<SourceItemType>();
+       ++source_item_id) {
+    auto target_item_list = item_to_item_matrix[source_item_id];
+    for (size_t i = 0; i < target_item_list.size() - 1; ++i) {
+      is_valid &= (item_number[target_item_list[i + 1]] > item_number[target_item_list[i]]);
+    }
+  }
+
+  return is_valid;
+}
+
+template <size_t Dimension>
+bool checkConnectivityOrdering(const Connectivity<Dimension>&);
+
+template <>
+bool
+checkConnectivityOrdering(const Connectivity<1>& connectivity)
+{
+  bool is_valid = true;
+  is_valid &= checkItemToHigherDimensionItem<1, ItemType::node, ItemType::cell>(connectivity);
+
+  return is_valid;
+}
+template <>
+
+bool
+checkConnectivityOrdering(const Connectivity<2>& connectivity)
+{
+  bool is_valid = true;
+  is_valid &= checkItemToHigherDimensionItem<2, ItemType::node, ItemType::face>(connectivity);
+  is_valid &= checkItemToHigherDimensionItem<2, ItemType::node, ItemType::cell>(connectivity);
+
+  is_valid &= checkItemToHigherDimensionItem<2, ItemType::face, ItemType::cell>(connectivity);
+
+  return is_valid;
+}
+
+bool
+checkConnectivityOrdering(const Connectivity<3>& connectivity)
+{
+  bool is_valid = true;
+  is_valid &= checkItemToHigherDimensionItem<3, ItemType::node, ItemType::edge>(connectivity);
+  is_valid &= checkItemToHigherDimensionItem<3, ItemType::node, ItemType::face>(connectivity);
+  is_valid &= checkItemToHigherDimensionItem<3, ItemType::node, ItemType::cell>(connectivity);
+
+  is_valid &= checkItemToHigherDimensionItem<3, ItemType::edge, ItemType::face>(connectivity);
+  is_valid &= checkItemToHigherDimensionItem<3, ItemType::edge, ItemType::cell>(connectivity);
+
+  is_valid &= checkItemToHigherDimensionItem<3, ItemType::face, ItemType::cell>(connectivity);
+
+  return parallel::allReduceAnd(is_valid);
+}
+
+bool
+checkConnectivityOrdering(const IConnectivity& connecticity)
+{
+  switch (connecticity.dimension()) {
+  case 1: {
+    return checkConnectivityOrdering(dynamic_cast<const Connectivity<1>&>(connecticity));
+  }
+  case 2: {
+    return checkConnectivityOrdering(dynamic_cast<const Connectivity<2>&>(connecticity));
+  }
+  case 3: {
+    return checkConnectivityOrdering(dynamic_cast<const Connectivity<3>&>(connecticity));
+  }
+    // LCOV_EXCL_START
+  default: {
+    throw UnexpectedError("invalid dimension");
+  }
+    // LCOV_EXCL_STOP
+  }
+}
diff --git a/src/mesh/ConnectivityUtils.hpp b/src/mesh/ConnectivityUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fa359a7e1c532e3b83c91d3e42d1bd8e3d7a45d8
--- /dev/null
+++ b/src/mesh/ConnectivityUtils.hpp
@@ -0,0 +1,8 @@
+#ifndef CONNECTIVITY_UTILS_HPP
+#define CONNECTIVITY_UTILS_HPP
+
+#include "mesh/IConnectivity.hpp"
+
+bool checkConnectivityOrdering(const IConnectivity&);
+
+#endif   // CONNECTIVITY_UTILS_HPP
diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 80fbc2661694bb2ed10f06726867f935951941c3..831246a7e3b3fc07ea55f1212d1a1b2044f55e4e 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -11,6 +11,7 @@
 #include <utils/Messenger.hpp>
 #include <utils/Stringify.hpp>
 
+#include <optional>
 #include <vector>
 
 template <size_t Dimension>
@@ -106,78 +107,87 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
     }
   });
 
-  diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
-
-  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) {
-    const size_t& i_diamond_cell      = face_id;
-    const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
-    const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
-    if (primal_face_cell_list.size() == 1) {
-      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 1);
-
-      const CellId cell_id      = primal_face_cell_list[0];
-      const auto i_face_in_cell = primal_face_local_number_in_their_cells(face_id, 0);
+  Array<const unsigned int> cell_to_node_row = [&] {
+    Array<unsigned int> tmp_cell_to_node_row(primal_number_of_faces + 1);
+    tmp_cell_to_node_row[0] = 0;
+    for (FaceId face_id = 0; face_id < primal_number_of_faces; ++face_id) {
+      tmp_cell_to_node_row[face_id + 1] = tmp_cell_to_node_row[face_id] + primal_face_to_node_matrix[face_id].size() +
+                                          primal_face_to_cell_matrix[face_id].size();
+    }
 
-      for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
-        diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node] = primal_face_node_list[i_node];
-      }
-      diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size()] =
-        primal_number_of_nodes + cell_id;
+    return tmp_cell_to_node_row;
+  }();
+
+  Array<const unsigned int> cell_to_node_list = [&] {
+    Array<unsigned int> tmp_cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+    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) {
+      const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
+      const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
+      const size_t first_node           = cell_to_node_row[face_id];
+      if (primal_face_cell_list.size() == 1) {
+        const CellId cell_id      = primal_face_cell_list[0];
+        const auto i_face_in_cell = primal_face_local_number_in_their_cells(face_id, 0);
 
-      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]);
+        for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+          tmp_cell_to_node_list[first_node + i_node] = primal_face_node_list[i_node];
         }
-      } 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
-                        .cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() - 1 - i_node]);
+        tmp_cell_to_node_list[first_node + primal_face_node_list.size()] = primal_number_of_nodes + cell_id;
+
+        if constexpr (Dimension == 2) {
+          if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
+            std::swap(tmp_cell_to_node_list[first_node], tmp_cell_to_node_list[first_node + 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(tmp_cell_to_node_list[first_node + i_node],
+                        tmp_cell_to_node_list[first_node + primal_face_node_list.size() - 1 - i_node]);
+            }
           }
-        }
-      }
-    } else {
-      Assert(primal_face_cell_list.size() == 2);
-      diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 2);
-
-      const CellId cell0_id      = primal_face_cell_list[0];
-      const CellId cell1_id      = primal_face_cell_list[1];
-      const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(face_id, 0);
-
-      if constexpr (Dimension == 2) {
-        Assert(primal_face_node_list.size() == 2);
-        diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
-        diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[0];
-        diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell1_id;
-        diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1];
-
-        if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
-          std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1],
-                    diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]);
         }
       } else {
-        diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
-        for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
-          diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node + 1] = primal_face_node_list[i_node];
-        }
-        diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] =
-          primal_number_of_nodes + cell1_id;
+        Assert(primal_face_cell_list.size() == 2);
+
+        const CellId cell0_id      = primal_face_cell_list[0];
+        const CellId cell1_id      = primal_face_cell_list[1];
+        const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(face_id, 0);
+
+        if constexpr (Dimension == 2) {
+          Assert(primal_face_node_list.size() == 2);
+
+          tmp_cell_to_node_list[first_node + 0] = primal_number_of_nodes + cell0_id;
+          tmp_cell_to_node_list[first_node + 1] = primal_face_node_list[0];
+          tmp_cell_to_node_list[first_node + 2] = primal_number_of_nodes + cell1_id;
+          tmp_cell_to_node_list[first_node + 3] = primal_face_node_list[1];
+
+          if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
+            std::swap(tmp_cell_to_node_list[first_node + 1], tmp_cell_to_node_list[first_node + 3]);
+          }
+        } else {
+          tmp_cell_to_node_list[first_node + 0] = primal_number_of_nodes + cell0_id;
+          for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+            tmp_cell_to_node_list[first_node + i_node + 1] = primal_face_node_list[i_node];
+          }
+          tmp_cell_to_node_list[first_node + primal_face_node_list.size() + 1] = primal_number_of_nodes + cell1_id;
 
-        if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
-          std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
-                    diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1]);
+          if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
+            std::swap(tmp_cell_to_node_list[first_node],
+                      tmp_cell_to_node_list[first_node + primal_face_node_list.size() + 1]);
+          }
         }
       }
-    }
-  });
+    });
+
+    return tmp_cell_to_node_list;
+  }();
+
+  diamond_descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row, cell_to_node_list);
 }
 
 template <size_t Dimension>
@@ -240,18 +250,28 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
   {
     const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
 
-    using Face = ConnectivityFace<Dimension>;
-
-    const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
-      std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
-      for (FaceId l = 0; l < diamond_descriptor.face_to_node_vector.size(); ++l) {
-        const auto& node_vector = diamond_descriptor.face_to_node_vector[l];
-
-        face_to_id_map[Face(node_vector, diamond_descriptor.node_number_vector)] = l;
+    const auto find_face = [&](std::vector<uint32_t> node_list) -> std::optional<FaceId> {
+      // The node list of already sorted correctly
+      const auto& face_id_vector = diamond_descriptor.node_to_face_matrix[node_list[0]];
+
+      for (size_t i_face = 0; i_face < face_id_vector.size(); ++i_face) {
+        const FaceId face_id          = face_id_vector[i_face];
+        const auto& face_node_id_list = diamond_descriptor.face_to_node_matrix[face_id];
+        if (face_node_id_list.size() == node_list.size()) {
+          bool is_same = true;
+          for (size_t i_node = 1; i_node < face_node_id_list.size(); ++i_node) {
+            is_same &= (diamond_descriptor.face_to_node_matrix[face_id][i_node] == node_list[i_node]);
+          }
+          if (is_same) {
+            return face_id;
+          }
+        }
       }
-      return face_to_id_map;
-    }();
 
+      return {};
+    };
+
+    std::vector<unsigned int> face_node_list;
     for (size_t i_face_list = 0; i_face_list < primal_connectivity.template numberOfRefItemList<ItemType::face>();
          ++i_face_list) {
       const auto& primal_ref_face_list = primal_connectivity.template refItemList<ItemType::face>(i_face_list);
@@ -265,16 +285,15 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
 
           const auto& primal_face_node_list = primal_face_to_node_matrix[primal_face_id];
 
-          const auto i_diamond_face = [&]() {
-            std::vector<unsigned int> node_list(primal_face_node_list.size());
-            for (size_t i = 0; i < primal_face_node_list.size(); ++i) {
-              node_list[i] = primal_face_node_list[i];
-            }
-            return face_to_id_map.find(Face(node_list, diamond_descriptor.node_number_vector));
-          }();
+          face_node_list.clear();
 
-          if (i_diamond_face != face_to_id_map.end()) {
-            diamond_face_list.push_back(i_diamond_face->second);
+          for (size_t i = 0; i < primal_face_node_list.size(); ++i) {
+            face_node_list.push_back(primal_face_node_list[i]);
+          }
+
+          auto face_id = find_face(face_node_list);
+          if (face_id.has_value()) {
+            diamond_face_list.push_back(face_id.value());
           }
         }
         return diamond_face_list;
@@ -293,19 +312,9 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
 
   if constexpr (Dimension == 3) {
     const auto& primal_edge_to_node_matrix = primal_connectivity.edgeToNodeMatrix();
-    using Edge                             = ConnectivityFace<2>;
-
-    const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map = [&] {
-      std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
-      for (EdgeId l = 0; l < diamond_descriptor.edge_to_node_vector.size(); ++l) {
-        const auto& node_vector = diamond_descriptor.edge_to_node_vector[l];
-        edge_to_id_map[Edge(node_vector, diamond_descriptor.node_number_vector)] = l;
-      }
-      return edge_to_id_map;
-    }();
 
     {
-      const size_t number_of_edges          = diamond_descriptor.edge_to_node_vector.size();
+      const size_t number_of_edges          = diamond_descriptor.edge_to_node_matrix.numberOfRows();
       diamond_descriptor.edge_number_vector = [&] {
         Array<int> edge_number_vector(number_of_edges);
         for (size_t i_edge = 0; i_edge < number_of_edges; ++i_edge) {
@@ -321,6 +330,22 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
       // LCOV_EXCL_STOP
     }
 
+    const auto find_edge = [&](uint32_t node0, uint32_t node1) -> std::optional<EdgeId> {
+      if (diamond_descriptor.node_number_vector[node0] > diamond_descriptor.node_number_vector[node1]) {
+        std::swap(node0, node1);
+      }
+      const auto& edge_id_vector = diamond_descriptor.node_to_edge_matrix[node0];
+
+      for (size_t i_edge = 0; i_edge < edge_id_vector.size(); ++i_edge) {
+        const EdgeId edge_id = edge_id_vector[i_edge];
+        if (diamond_descriptor.edge_to_node_matrix[edge_id][1] == node1) {
+          return edge_id;
+        }
+      }
+
+      return {};
+    };
+
     for (size_t i_edge_list = 0; i_edge_list < primal_connectivity.template numberOfRefItemList<ItemType::edge>();
          ++i_edge_list) {
       const auto& primal_ref_edge_list = primal_connectivity.template refItemList<ItemType::edge>(i_edge_list);
@@ -334,16 +359,10 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
 
           const auto& primal_edge_node_list = primal_edge_to_node_matrix[primal_edge_id];
 
-          const auto i_diamond_edge = [&]() {
-            std::vector<unsigned int> node_list(primal_edge_node_list.size());
-            for (size_t i = 0; i < primal_edge_node_list.size(); ++i) {
-              node_list[i] = primal_edge_node_list[i];
-            }
-            return edge_to_id_map.find(Edge(node_list, diamond_descriptor.node_number_vector));
-          }();
+          const auto diamond_edge_id = find_edge(primal_edge_node_list[0], primal_edge_node_list[1]);
 
-          if (i_diamond_edge != edge_to_id_map.end()) {
-            diamond_edge_list.push_back(i_diamond_edge->second);
+          if (diamond_edge_id.has_value()) {
+            diamond_edge_list.push_back(diamond_edge_id.value());
           }
         }
         return diamond_edge_list;
@@ -388,8 +407,8 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
     diamond_descriptor.face_owner_vector = Array<int>(diamond_descriptor.face_number_vector.size());
     diamond_descriptor.face_owner_vector.fill(parallel::size());
 
-    for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
-      const auto& cell_face_list = diamond_descriptor.cell_to_face_vector[i_cell];
+    for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_matrix.numberOfRows(); ++i_cell) {
+      const auto& cell_face_list = diamond_descriptor.cell_to_face_matrix[i_cell];
       for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
         const size_t face_id = cell_face_list[i_face];
         diamond_descriptor.face_owner_vector[face_id] =
@@ -402,8 +421,8 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
     diamond_descriptor.edge_owner_vector = Array<int>(diamond_descriptor.edge_number_vector.size());
     diamond_descriptor.edge_owner_vector.fill(parallel::size());
 
-    for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
-      const auto& cell_edge_list = diamond_descriptor.cell_to_edge_vector[i_cell];
+    for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_matrix.numberOfRows(); ++i_cell) {
+      const auto& cell_edge_list = diamond_descriptor.cell_to_edge_matrix[i_cell];
       for (size_t i_edge = 0; i_edge < cell_edge_list.size(); ++i_edge) {
         const size_t edge_id = cell_edge_list[i_edge];
         diamond_descriptor.edge_owner_vector[edge_id] =
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 21e413ff4e94519b79fbc26e87d742e67c0436c2..64711edc72055effd9e5608bc01b71ea184ee1fc 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -43,8 +43,6 @@ DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const IMesh& i_mesh)
 
 DiamondDualMeshBuilder::DiamondDualMeshBuilder(const IMesh& i_mesh)
 {
-  std::cout << "building DiamondDualMesh\n";
-
   switch (i_mesh.dimension()) {
   case 2: {
     this->_buildDualDiamondMeshFrom<2>(i_mesh);
diff --git a/src/mesh/Dual1DConnectivityBuilder.cpp b/src/mesh/Dual1DConnectivityBuilder.cpp
index 4a4137b4998c852b0176207d1179d65c24aacfcf..d5ca12372dee5ecbbebb7c03e3824129a9334578 100644
--- a/src/mesh/Dual1DConnectivityBuilder.cpp
+++ b/src/mesh/Dual1DConnectivityBuilder.cpp
@@ -60,30 +60,36 @@ Dual1DConnectivityBuilder::_buildConnectivityDescriptor(const Connectivity<1>& p
   dual_descriptor.cell_type_vector = Array<CellType>(dual_number_of_cells);
   dual_descriptor.cell_type_vector.fill(CellType::Line);
 
-  dual_descriptor.cell_to_node_vector.resize(dual_number_of_cells);
-
   const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells();
 
+  Array<unsigned int> cell_to_node_row(dual_number_of_cells + 1);
+  parallel_for(
+    cell_to_node_row.size(), PUGS_LAMBDA(const CellId cell_id) { cell_to_node_row[cell_id] = 2 * cell_id; });
+
+  Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
   {
     size_t next_kept_node_id = 0;
+    size_t i_cell_node       = 0;
     for (NodeId i_node = 0; i_node < primal_connectivity.numberOfNodes(); ++i_node) {
-      const size_t& i_dual_cell         = i_node;
       const auto& primal_node_cell_list = primal_node_to_cell_matrix[i_node];
-      dual_descriptor.cell_to_node_vector[i_dual_cell].resize(2);
       if (primal_node_cell_list.size() == 1) {
         const auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0);
 
-        dual_descriptor.cell_to_node_vector[i_dual_cell][i_node_in_cell] = next_kept_node_id++;
-        dual_descriptor.cell_to_node_vector[i_dual_cell][1 - i_node_in_cell] =
-          number_of_kept_nodes + primal_node_cell_list[0];
+        cell_to_node_list[i_cell_node + i_node_in_cell]     = next_kept_node_id++;
+        cell_to_node_list[i_cell_node + 1 - i_node_in_cell] = number_of_kept_nodes + primal_node_cell_list[0];
+
+        i_cell_node += 2;
+
       } else {
         const auto i0 = primal_node_local_number_in_their_cells(i_node, 0);
 
-        dual_descriptor.cell_to_node_vector[i_dual_cell][0] = number_of_kept_nodes + primal_node_cell_list[1 - i0];
-        dual_descriptor.cell_to_node_vector[i_dual_cell][1] = number_of_kept_nodes + primal_node_cell_list[i0];
+        cell_to_node_list[i_cell_node++] = number_of_kept_nodes + primal_node_cell_list[1 - i0];
+        cell_to_node_list[i_cell_node++] = number_of_kept_nodes + primal_node_cell_list[i0];
       }
     }
   }
+
+  dual_descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row, cell_to_node_list);
 }
 
 void
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index f2030623a72e2262ce5dbb391caac0c304f846da..95bb289d3c53e3a12207b1b9be6eb7764eec749e 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -38,14 +38,22 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
   descriptor.node_number_vector = convert_to_array(gmsh_data.__verticesNumbers);
   descriptor.cell_type_vector   = Array<CellType>(nb_cells);
   descriptor.cell_number_vector = Array<int>(nb_cells);
-  descriptor.cell_to_node_vector.resize(nb_cells);
 
-  for (size_t j = 0; j < nb_cells; ++j) {
-    descriptor.cell_to_node_vector[j].resize(2);
-    for (int r = 0; r < 2; ++r) {
-      descriptor.cell_to_node_vector[j][r] = gmsh_data.__edges[j][r];
+  Array<unsigned int> cell_to_node_row(nb_cells + 1);
+  parallel_for(
+    cell_to_node_row.size(), PUGS_LAMBDA(const CellId cell_id) { cell_to_node_row[cell_id] = 2 * cell_id; });
+
+  Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+  for (CellId cell_id = 0; cell_id < nb_cells; ++cell_id) {
+    for (int i_node = 0; i_node < 2; ++i_node) {
+      cell_to_node_list[2 * cell_id + i_node] = gmsh_data.__edges[cell_id][i_node];
     }
-    descriptor.cell_type_vector[j]   = CellType::Line;
+  }
+  descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row, cell_to_node_list);
+
+  descriptor.cell_type_vector.fill(CellType::Line);
+
+  for (size_t j = 0; j < nb_cells; ++j) {
     descriptor.cell_number_vector[j] = gmsh_data.__edges_number[j];
   }
 
@@ -61,7 +69,7 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (size_t j = 0; j < nb_cells; ++j) {
     for (int r = 0; r < 2; ++r) {
-      node_nb_cell[descriptor.cell_to_node_vector[j][r]] += 1;
+      node_nb_cell[descriptor.cell_to_node_matrix[j][r]] += 1;
     }
   }
 
@@ -133,27 +141,48 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
   descriptor.node_number_vector = convert_to_array(gmsh_data.__verticesNumbers);
   descriptor.cell_type_vector   = Array<CellType>(nb_cells);
   descriptor.cell_number_vector = Array<int>(nb_cells);
-  descriptor.cell_to_node_vector.resize(nb_cells);
 
-  const size_t nb_triangles = gmsh_data.__triangles.size();
-  for (size_t j = 0; j < nb_triangles; ++j) {
-    descriptor.cell_to_node_vector[j].resize(3);
-    for (int r = 0; r < 3; ++r) {
-      descriptor.cell_to_node_vector[j][r] = gmsh_data.__triangles[j][r];
+  Array<unsigned int> cell_to_node_row(gmsh_data.__triangles.size() + gmsh_data.__quadrangles.size() + 1);
+  {
+    cell_to_node_row[0] = 0;
+    size_t i_cell       = 0;
+    for (size_t i_triangle = 0; i_triangle < gmsh_data.__triangles.size(); ++i_triangle, ++i_cell) {
+      cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 3;
+    }
+    for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.__quadrangles.size(); ++i_quadrangle, ++i_cell) {
+      cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 4;
     }
-    descriptor.cell_type_vector[j]   = CellType::Triangle;
-    descriptor.cell_number_vector[j] = gmsh_data.__triangles_number[j];
+  }
+
+  Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+  {
+    size_t i_cell_node = 0;
+    for (size_t i_triangle = 0; i_triangle < gmsh_data.__triangles.size(); ++i_triangle) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][2];
+    }
+    for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.__quadrangles.size(); ++i_quadrangle) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][3];
+    }
+  }
+  descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row, cell_to_node_list);
+
+  const size_t nb_triangles = gmsh_data.__triangles.size();
+  for (size_t i_triangle = 0; i_triangle < nb_triangles; ++i_triangle) {
+    descriptor.cell_type_vector[i_triangle]   = CellType::Triangle;
+    descriptor.cell_number_vector[i_triangle] = gmsh_data.__triangles_number[i_triangle];
   }
 
   const size_t nb_quadrangles = gmsh_data.__quadrangles.size();
-  for (size_t j = 0; j < nb_quadrangles; ++j) {
-    const size_t jq = j + nb_triangles;
-    descriptor.cell_to_node_vector[jq].resize(4);
-    for (int r = 0; r < 4; ++r) {
-      descriptor.cell_to_node_vector[jq][r] = gmsh_data.__quadrangles[j][r];
-    }
-    descriptor.cell_type_vector[jq]   = CellType::Quadrangle;
-    descriptor.cell_number_vector[jq] = gmsh_data.__quadrangles_number[j];
+  for (size_t i_quadrangle = 0; i_quadrangle < nb_quadrangles; ++i_quadrangle) {
+    const size_t i_cell = i_quadrangle + nb_triangles;
+
+    descriptor.cell_type_vector[i_cell]   = CellType::Quadrangle;
+    descriptor.cell_number_vector[i_cell] = gmsh_data.__quadrangles_number[i_quadrangle];
   }
 
   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
@@ -189,17 +218,6 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
 
-  using Face                                                                 = ConnectivityFace<2>;
-  const auto& node_number_vector                                             = descriptor.node_number_vector;
-  const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
-    std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
-    for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-      const auto& node_vector                               = descriptor.face_to_node_vector[l];
-      face_to_id_map[Face(node_vector, node_number_vector)] = l;
-    }
-    return face_to_id_map;
-  }();
-
   std::unordered_map<int, FaceId> face_number_id_map = [&] {
     std::unordered_map<int, FaceId> face_number_id_map;
     for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
@@ -209,17 +227,28 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     return face_number_id_map;
   }();
 
+  const auto find_face = [&](uint32_t node0, uint32_t node1) {
+    if (descriptor.node_number_vector[node0] > descriptor.node_number_vector[node1]) {
+      std::swap(node0, node1);
+    }
+    const auto& face_id_vector = descriptor.node_to_face_matrix[node0];
+
+    for (size_t i_face = 0; i_face < face_id_vector.size(); ++i_face) {
+      const FaceId face_id = face_id_vector[i_face];
+      if (descriptor.face_to_node_matrix[face_id][1] == node1) {
+        return face_id;
+      }
+    }
+
+    std::stringstream error_msg;
+    error_msg << "face (" << node0 << ',' << node1 << ") not found";
+    throw NormalError(error_msg.str());
+  };
+
   std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
   for (unsigned int e = 0; e < gmsh_data.__edges.size(); ++e) {
-    const unsigned int edge_id = [&] {
-      auto i = face_to_id_map.find(Face({gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]}, node_number_vector));
-      if (i == face_to_id_map.end()) {
-        std::stringstream error_msg;
-        error_msg << "face " << gmsh_data.__edges[e][0] << " not found";
-        throw NormalError(error_msg.str());
-      }
-      return i->second;
-    }();
+    const unsigned int edge_id = find_face(gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]);
+
     const unsigned int& ref = gmsh_data.__edges_ref[e];
     ref_faces_map[ref].push_back(edge_id);
 
@@ -244,9 +273,9 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
   Array<size_t> face_nb_cell(descriptor.face_number_vector.size());
   face_nb_cell.fill(0);
 
-  for (size_t j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
-    for (size_t l = 0; l < descriptor.cell_to_face_vector[j].size(); ++l) {
-      face_nb_cell[descriptor.cell_to_face_vector[j][l]] += 1;
+  for (size_t j = 0; j < descriptor.cell_to_face_matrix.numberOfRows(); ++j) {
+    for (size_t l = 0; l < descriptor.cell_to_face_matrix[j].size(); ++l) {
+      face_nb_cell[descriptor.cell_to_face_matrix[j][l]] += 1;
     }
   }
 
@@ -281,8 +310,8 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (size_t i_face = 0; i_face < face_nb_cell.size(); ++i_face) {
     if (face_nb_cell[i_face] == 1) {
-      for (size_t node_id : descriptor.face_to_node_vector[i_face]) {
-        is_boundary_node[node_id] = true;
+      for (size_t i_node = 0; i_node < descriptor.face_to_node_matrix[i_face].size(); ++i_node) {
+        is_boundary_node[descriptor.face_to_node_matrix[i_face][i_node]] = true;
       }
     }
   }
@@ -339,51 +368,94 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
   descriptor.node_number_vector = convert_to_array(gmsh_data.__verticesNumbers);
   descriptor.cell_type_vector   = Array<CellType>(nb_cells);
   descriptor.cell_number_vector = Array<int>(nb_cells);
-  descriptor.cell_to_node_vector.resize(nb_cells);
 
   const size_t nb_tetrahedra = gmsh_data.__tetrahedra.size();
-  for (size_t j = 0; j < nb_tetrahedra; ++j) {
-    descriptor.cell_to_node_vector[j].resize(4);
-    for (int r = 0; r < 4; ++r) {
-      descriptor.cell_to_node_vector[j][r] = gmsh_data.__tetrahedra[j][r];
+  const size_t nb_hexahedra  = gmsh_data.__hexahedra.size();
+  const size_t nb_prisms     = gmsh_data.__prisms.size();
+  const size_t nb_pyramids   = gmsh_data.__pyramids.size();
+
+  Array<unsigned int> cell_to_node_row(nb_cells + 1);
+  {
+    cell_to_node_row[0] = 0;
+    size_t i_cell       = 0;
+    for (size_t i_tetrahedron = 0; i_tetrahedron < nb_tetrahedra; ++i_tetrahedron, ++i_cell) {
+      cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 4;
+    }
+    for (size_t i_hexahedron = 0; i_hexahedron < nb_hexahedra; ++i_hexahedron, ++i_cell) {
+      cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 8;
+    }
+    for (size_t i_prism = 0; i_prism < nb_prisms; ++i_prism, ++i_cell) {
+      cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 6;
+    }
+    for (size_t i_pyramid = 0; i_pyramid < nb_pyramids; ++i_pyramid, ++i_cell) {
+      cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 5;
     }
+  }
+
+  Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+  {
+    size_t i_cell_node = 0;
+    for (size_t i_tetrahedron = 0; i_tetrahedron < nb_tetrahedra; ++i_tetrahedron) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][3];
+    }
+    for (size_t i_hexahedron = 0; i_hexahedron < nb_hexahedra; ++i_hexahedron) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][3];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][4];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][5];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][6];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][7];
+    }
+    for (size_t i_prism = 0; i_prism < nb_prisms; ++i_prism) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][3];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][4];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][5];
+    }
+    for (size_t i_pyramid = 0; i_pyramid < nb_pyramids; ++i_pyramid) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][3];
+      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][4];
+    }
+  }
+
+  for (size_t j = 0; j < nb_tetrahedra; ++j) {
     descriptor.cell_type_vector[j]   = CellType::Tetrahedron;
     descriptor.cell_number_vector[j] = gmsh_data.__tetrahedra_number[j];
   }
 
-  const size_t nb_hexahedra = gmsh_data.__hexahedra.size();
   for (size_t j = 0; j < nb_hexahedra; ++j) {
     const size_t jh = nb_tetrahedra + j;
-    descriptor.cell_to_node_vector[jh].resize(8);
-    for (int r = 0; r < 8; ++r) {
-      descriptor.cell_to_node_vector[jh][r] = gmsh_data.__hexahedra[j][r];
-    }
+
     descriptor.cell_type_vector[jh]   = CellType::Hexahedron;
     descriptor.cell_number_vector[jh] = gmsh_data.__hexahedra_number[j];
   }
 
-  const size_t nb_prisms = gmsh_data.__prisms.size();
   for (size_t j = 0; j < nb_prisms; ++j) {
     const size_t jp = nb_tetrahedra + nb_hexahedra + j;
-    descriptor.cell_to_node_vector[jp].resize(6);
-    for (int r = 0; r < 6; ++r) {
-      descriptor.cell_to_node_vector[jp][r] = gmsh_data.__prisms[j][r];
-    }
+
     descriptor.cell_type_vector[jp]   = CellType::Prism;
     descriptor.cell_number_vector[jp] = gmsh_data.__prisms_number[j];
   }
 
-  const size_t nb_pyramids = gmsh_data.__pyramids.size();
   for (size_t j = 0; j < nb_pyramids; ++j) {
     const size_t jh = nb_tetrahedra + nb_hexahedra + nb_prisms + j;
-    descriptor.cell_to_node_vector[jh].resize(5);
-    for (int r = 0; r < 5; ++r) {
-      descriptor.cell_to_node_vector[jh][r] = gmsh_data.__pyramids[j][r];
-    }
+
     descriptor.cell_type_vector[jh]   = CellType::Pyramid;
     descriptor.cell_number_vector[jh] = gmsh_data.__pyramids_number[j];
   }
 
+  descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row, cell_to_node_list);
+
   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
   for (unsigned int j = 0; j < gmsh_data.__tetrahedra_ref.size(); ++j) {
     const unsigned int elem_number = j;
@@ -434,23 +506,13 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
   Array<size_t> face_nb_cell(descriptor.face_number_vector.size());
   face_nb_cell.fill(0);
 
-  for (size_t j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
-    for (size_t l = 0; l < descriptor.cell_to_face_vector[j].size(); ++l) {
-      face_nb_cell[descriptor.cell_to_face_vector[j][l]] += 1;
+  for (size_t j = 0; j < descriptor.cell_to_face_matrix.numberOfRows(); ++j) {
+    for (size_t l = 0; l < descriptor.cell_to_face_matrix[j].size(); ++l) {
+      face_nb_cell[descriptor.cell_to_face_matrix[j][l]] += 1;
     }
   }
 
   {
-    using Face                                                                 = ConnectivityFace<3>;
-    const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
-      std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
-      for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-        const auto& node_vector                               = descriptor.face_to_node_vector[l];
-        face_to_id_map[Face(node_vector, node_number_vector)] = l;
-      }
-      return face_to_id_map;
-    }();
-
     std::unordered_map<int, FaceId> face_number_id_map = [&] {
       std::unordered_map<int, FaceId> face_number_id_map;
       for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
@@ -460,17 +522,63 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       return face_number_id_map;
     }();
 
+    const auto find_face = [&](std::vector<uint32_t> node_list) {
+      size_t i_node_smallest_number = 0;
+      for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+        if (node_number_vector[node_list[i_node]] < node_number_vector[node_list[i_node_smallest_number]]) {
+          i_node_smallest_number = i_node;
+        }
+      }
+
+      if (i_node_smallest_number != 0) {
+        std::vector<uint64_t> buffer(node_list.size());
+        for (size_t i_node = i_node_smallest_number; i_node < buffer.size(); ++i_node) {
+          buffer[i_node - i_node_smallest_number] = node_list[i_node];
+        }
+        for (size_t i_node = 0; i_node < i_node_smallest_number; ++i_node) {
+          buffer[i_node + node_list.size() - i_node_smallest_number] = node_list[i_node];
+        }
+
+        for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+          node_list[i_node] = buffer[i_node];
+        }
+      }
+
+      if (node_number_vector[node_list[1]] > node_number_vector[node_list[node_list.size() - 1]]) {
+        for (size_t i_node = 1; i_node <= (node_list.size() + 1) / 2 - 1; ++i_node) {
+          std::swap(node_list[i_node], node_list[node_list.size() - i_node]);
+        }
+      }
+
+      const auto& face_id_vector = descriptor.node_to_face_matrix[node_list[0]];
+
+      for (size_t i_face = 0; i_face < face_id_vector.size(); ++i_face) {
+        const FaceId face_id          = face_id_vector[i_face];
+        const auto& face_node_id_list = descriptor.face_to_node_matrix[face_id];
+        if (face_node_id_list.size() == node_list.size()) {
+          bool is_same = true;
+          for (size_t i_node = 1; i_node < face_node_id_list.size(); ++i_node) {
+            is_same &= (descriptor.face_to_node_matrix[face_id][i_node] == node_list[i_node]);
+          }
+          if (is_same) {
+            return face_id;
+          }
+        }
+      }
+
+      std::stringstream error_msg;
+      error_msg << "face (" << node_list[0];
+      for (size_t i = 1; i < node_list.size(); ++i) {
+        error_msg << ',' << node_list[i];
+      }
+      error_msg << ") not found";
+      throw NormalError(error_msg.str());
+    };
+
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
     for (unsigned int f = 0; f < gmsh_data.__triangles.size(); ++f) {
-      const unsigned int face_id = [&] {
-        auto i = face_to_id_map.find(
-          Face({gmsh_data.__triangles[f][0], gmsh_data.__triangles[f][1], gmsh_data.__triangles[f][2]},
-               node_number_vector));
-        if (i == face_to_id_map.end()) {
-          throw NormalError("face not found");
-        }
-        return i->second;
-      }();
+      const unsigned int face_id =
+        find_face({gmsh_data.__triangles[f][0], gmsh_data.__triangles[f][1], gmsh_data.__triangles[f][2]});
 
       const unsigned int& ref = gmsh_data.__triangles_ref[f];
       ref_faces_map[ref].push_back(face_id);
@@ -495,15 +603,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     }
 
     for (unsigned int f = 0; f < gmsh_data.__quadrangles.size(); ++f) {
-      const unsigned int face_id = [&] {
-        auto i = face_to_id_map.find(Face({gmsh_data.__quadrangles[f][0], gmsh_data.__quadrangles[f][1],
-                                           gmsh_data.__quadrangles[f][2], gmsh_data.__quadrangles[f][3]},
-                                          node_number_vector));
-        if (i == face_to_id_map.end()) {
-          throw NormalError("face not found");
-        }
-        return i->second;
-      }();
+      const unsigned int face_id = find_face({gmsh_data.__quadrangles[f][0], gmsh_data.__quadrangles[f][1],
+                                              gmsh_data.__quadrangles[f][2], gmsh_data.__quadrangles[f][3]});
 
       const unsigned int& ref = gmsh_data.__quadrangles_ref[f];
       ref_faces_map[ref].push_back(face_id);
@@ -561,22 +662,30 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
     for (size_t i_face = 0; i_face < face_nb_cell.size(); ++i_face) {
       if (face_nb_cell[i_face] == 1) {
-        for (size_t edge_id : descriptor.face_to_edge_vector[i_face]) {
-          is_boundary_edge[edge_id] = true;
+        auto face_edges = descriptor.face_to_edge_matrix[i_face];
+        for (size_t i_edge = 0; i_edge < face_edges.size(); ++i_edge) {
+          is_boundary_edge[face_edges[i_edge]] = true;
         }
       }
     }
 
-    using Edge                                                                 = ConnectivityFace<2>;
-    const auto& node_number_vector                                             = descriptor.node_number_vector;
-    const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map = [&] {
-      std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
-      for (EdgeId l = 0; l < descriptor.edge_to_node_vector.size(); ++l) {
-        const auto& node_vector                               = descriptor.edge_to_node_vector[l];
-        edge_to_id_map[Edge(node_vector, node_number_vector)] = l;
+    const auto find_edge = [&](uint32_t node0, uint32_t node1) {
+      if (descriptor.node_number_vector[node0] > descriptor.node_number_vector[node1]) {
+        std::swap(node0, node1);
       }
-      return edge_to_id_map;
-    }();
+      const auto& edge_id_vector = descriptor.node_to_edge_matrix[node0];
+
+      for (size_t i_edge = 0; i_edge < edge_id_vector.size(); ++i_edge) {
+        const EdgeId edge_id = edge_id_vector[i_edge];
+        if (descriptor.edge_to_node_matrix[edge_id][1] == node1) {
+          return edge_id;
+        }
+      }
+
+      std::stringstream error_msg;
+      error_msg << "edge (" << node0 << ',' << node1 << ") not found";
+      throw NormalError(error_msg.str());
+    };
 
     std::unordered_map<int, EdgeId> edge_number_id_map = [&] {
       std::unordered_map<int, EdgeId> edge_number_id_map;
@@ -589,16 +698,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
     std::map<unsigned int, std::vector<unsigned int>> ref_edges_map;
     for (unsigned int e = 0; e < gmsh_data.__edges.size(); ++e) {
-      const unsigned int edge_id = [&] {
-        auto i = edge_to_id_map.find(Edge({gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]}, node_number_vector));
-        if (i == edge_to_id_map.end()) {
-          std::stringstream error_msg;
-          error_msg << "edge " << gmsh_data.__edges[e][0] << " not found";
-          throw NormalError(error_msg.str());
-        }
-        return i->second;
-      }();
-      const unsigned int& ref = gmsh_data.__edges_ref[e];
+      const unsigned int edge_id = find_edge(gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]);
+      const unsigned int& ref    = gmsh_data.__edges_ref[e];
       ref_edges_map[ref].push_back(edge_id);
 
       if (descriptor.edge_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
@@ -650,7 +751,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (size_t i_face = 0; i_face < face_nb_cell.size(); ++i_face) {
     if (face_nb_cell[i_face] == 1) {
-      for (size_t node_id : descriptor.face_to_node_vector[i_face]) {
+      for (size_t i_node = 0; i_node < descriptor.face_to_node_matrix[i_face].size(); ++i_node) {
+        const NodeId node_id      = descriptor.face_to_node_matrix[i_face][i_node];
         is_boundary_node[node_id] = true;
       }
     }
@@ -958,255 +1060,6 @@ GmshReader::__readPeriodic2_2()
   }
 }
 
-// std::shared_ptr<IConnectivity>
-// GmshReader::_buildConnectivity3D(const size_t nb_cells)
-// {
-//   ConnectivityDescriptor descriptor;
-
-//   descriptor.node_number_vector = m_mesh_data.__verticesNumbers;
-//   descriptor.cell_type_vector.resize(nb_cells);
-//   descriptor.cell_number_vector.resize(nb_cells);
-//   descriptor.cell_to_node_vector.resize(nb_cells);
-
-//   const size_t nb_tetrahedra = m_mesh_data.__tetrahedra.size();
-//   for (size_t j = 0; j < nb_tetrahedra; ++j) {
-//     descriptor.cell_to_node_vector[j].resize(4);
-//     for (int r = 0; r < 4; ++r) {
-//       descriptor.cell_to_node_vector[j][r] = m_mesh_data.__tetrahedra[j][r];
-//     }
-//     descriptor.cell_type_vector[j]   = CellType::Tetrahedron;
-//     descriptor.cell_number_vector[j] = m_mesh_data.__tetrahedra_number[j];
-//   }
-//   const size_t nb_hexahedra = m_mesh_data.__hexahedra.size();
-//   for (size_t j = 0; j < nb_hexahedra; ++j) {
-//     const size_t jh = nb_tetrahedra + j;
-//     descriptor.cell_to_node_vector[jh].resize(8);
-//     for (int r = 0; r < 8; ++r) {
-//       descriptor.cell_to_node_vector[jh][r] = m_mesh_data.__hexahedra[j][r];
-//     }
-//     descriptor.cell_type_vector[jh]   = CellType::Hexahedron;
-//     descriptor.cell_number_vector[jh] = m_mesh_data.__hexahedra_number[j];
-//   }
-
-//   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
-//   for (unsigned int r = 0; r < m_mesh_data.__tetrahedra_ref.size(); ++r) {
-//     const unsigned int elem_number = m_mesh_data.__tetrahedra_ref[r];
-//     const unsigned int& ref        = m_mesh_data.__tetrahedra_ref[r];
-//     ref_cells_map[ref].push_back(elem_number);
-//   }
-
-//   for (unsigned int j = 0; j < m_mesh_data.__hexahedra_ref.size(); ++j) {
-//     const size_t elem_number = nb_tetrahedra + j;
-//     const unsigned int& ref  = m_mesh_data.__hexahedra_ref[j];
-//     ref_cells_map[ref].push_back(elem_number);
-//   }
-
-//   for (const auto& ref_cell_list : ref_cells_map) {
-//     Array<CellId> cell_list(ref_cell_list.second.size());
-//     for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
-//       cell_list[j] = ref_cell_list.second[j];
-//     }
-//     const PhysicalRefId& physical_ref_id = m_mesh_data.m_physical_ref_map.at(ref_cell_list.first);
-//     descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
-//   }
-
-//   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
-
-//   const auto& node_number_vector = descriptor.node_number_vector;
-
-//   {
-//     using Face                                                                 = ConnectivityFace<3>;
-//     const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
-//       std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
-//       for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-//         const auto& node_vector                               = descriptor.face_to_node_vector[l];
-//         face_to_id_map[Face(node_vector, node_number_vector)] = l;
-//       }
-//       return face_to_id_map;
-//     }();
-
-//     std::unordered_map<int, FaceId> face_number_id_map = [&] {
-//       std::unordered_map<int, FaceId> face_number_id_map;
-//       for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
-//         face_number_id_map[descriptor.face_number_vector[l]] = l;
-//       }
-//       Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
-//       return face_number_id_map;
-//     }();
-
-//     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
-//     for (unsigned int f = 0; f < m_mesh_data.__triangles.size(); ++f) {
-//       const unsigned int face_id = [&] {
-//         auto i = face_to_id_map.find(
-//           Face({m_mesh_data.__triangles[f][0], m_mesh_data.__triangles[f][1], m_mesh_data.__triangles[f][2]},
-//                node_number_vector));
-//         if (i == face_to_id_map.end()) {
-//           throw NormalError("face not found");
-//         }
-//         return i->second;
-//       }();
-
-//       const unsigned int& ref = m_mesh_data.__triangles_ref[f];
-//       ref_faces_map[ref].push_back(face_id);
-
-//       if (descriptor.face_number_vector[face_id] != m_mesh_data.__quadrangles_number[f]) {
-//         if (auto i_face = face_number_id_map.find(m_mesh_data.__quadrangles_number[f]);
-//             i_face != face_number_id_map.end()) {
-//           const int other_face_id = i_face->second;
-//           std::swap(descriptor.face_number_vector[face_id], descriptor.face_number_vector[other_face_id]);
-
-//           face_number_id_map.erase(descriptor.face_number_vector[face_id]);
-//           face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
-
-//           face_number_id_map[descriptor.face_number_vector[face_id]]       = face_id;
-//           face_number_id_map[descriptor.face_number_vector[other_face_id]] = other_face_id;
-//         } else {
-//           face_number_id_map.erase(descriptor.face_number_vector[face_id]);
-//           descriptor.face_number_vector[face_id]                     = m_mesh_data.__quadrangles_number[f];
-//           face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
-//         }
-//       }
-//     }
-
-//     for (unsigned int f = 0; f < m_mesh_data.__quadrangles.size(); ++f) {
-//       const unsigned int face_id = [&] {
-//         auto i = face_to_id_map.find(Face({m_mesh_data.__quadrangles[f][0], m_mesh_data.__quadrangles[f][1],
-//                                            m_mesh_data.__quadrangles[f][2], m_mesh_data.__quadrangles[f][3]},
-//                                           node_number_vector));
-//         if (i == face_to_id_map.end()) {
-//           throw NormalError("face not found");
-//         }
-//         return i->second;
-//       }();
-
-//       const unsigned int& ref = m_mesh_data.__quadrangles_ref[f];
-//       ref_faces_map[ref].push_back(face_id);
-
-//       if (descriptor.face_number_vector[face_id] != m_mesh_data.__quadrangles_number[f]) {
-//         if (auto i_face = face_number_id_map.find(m_mesh_data.__quadrangles_number[f]);
-//             i_face != face_number_id_map.end()) {
-//           const int other_face_id = i_face->second;
-//           std::swap(descriptor.face_number_vector[face_id], descriptor.face_number_vector[other_face_id]);
-
-//           face_number_id_map.erase(descriptor.face_number_vector[face_id]);
-//           face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
-
-//           face_number_id_map[descriptor.face_number_vector[face_id]]       = face_id;
-//           face_number_id_map[descriptor.face_number_vector[other_face_id]] = other_face_id;
-//         } else {
-//           face_number_id_map.erase(descriptor.face_number_vector[face_id]);
-//           descriptor.face_number_vector[face_id]                     = m_mesh_data.__quadrangles_number[f];
-//           face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
-//         }
-//       }
-//     }
-
-//     for (const auto& ref_face_list : ref_faces_map) {
-//       Array<FaceId> face_list(ref_face_list.second.size());
-//       for (size_t j = 0; j < ref_face_list.second.size(); ++j) {
-//         face_list[j] = ref_face_list.second[j];
-//       }
-//       const PhysicalRefId& physical_ref_id = m_mesh_data.m_physical_ref_map.at(ref_face_list.first);
-//       descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
-//     }
-//   }
-
-//   ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(descriptor);
-
-//   {
-//     using Edge                                                                 = ConnectivityFace<2>;
-//     const auto& node_number_vector                                             = descriptor.node_number_vector;
-//     const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map = [&] {
-//       std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
-//       for (EdgeId l = 0; l < descriptor.edge_to_node_vector.size(); ++l) {
-//         const auto& node_vector                               = descriptor.edge_to_node_vector[l];
-//         edge_to_id_map[Edge(node_vector, node_number_vector)] = l;
-//       }
-//       return edge_to_id_map;
-//     }();
-
-//     std::unordered_map<int, EdgeId> edge_number_id_map = [&] {
-//       std::unordered_map<int, EdgeId> edge_number_id_map;
-//       for (size_t l = 0; l < descriptor.edge_number_vector.size(); ++l) {
-//         edge_number_id_map[descriptor.edge_number_vector[l]] = l;
-//       }
-//       Assert(edge_number_id_map.size() == descriptor.edge_number_vector.size());
-//       return edge_number_id_map;
-//     }();
-
-//     std::map<unsigned int, std::vector<unsigned int>> ref_edges_map;
-//     for (unsigned int e = 0; e < m_mesh_data.__edges.size(); ++e) {
-//       const unsigned int edge_id = [&] {
-//         auto i = edge_to_id_map.find(Edge({m_mesh_data.__edges[e][0], m_mesh_data.__edges[e][1]},
-//         node_number_vector)); if (i == edge_to_id_map.end()) {
-//           std::stringstream error_msg;
-//           error_msg << "edge " << m_mesh_data.__edges[e][0] << " not found";
-//           throw NormalError(error_msg.str());
-//         }
-//         return i->second;
-//       }();
-//       const unsigned int& ref = m_mesh_data.__edges_ref[e];
-//       ref_edges_map[ref].push_back(edge_id);
-
-//       if (descriptor.edge_number_vector[edge_id] != m_mesh_data.__edges_number[e]) {
-//         if (auto i_edge = edge_number_id_map.find(m_mesh_data.__edges_number[e]); i_edge != edge_number_id_map.end())
-//         {
-//           const int other_edge_id = i_edge->second;
-//           std::swap(descriptor.edge_number_vector[edge_id], descriptor.edge_number_vector[other_edge_id]);
-
-//           edge_number_id_map.erase(descriptor.edge_number_vector[edge_id]);
-//           edge_number_id_map.erase(descriptor.edge_number_vector[other_edge_id]);
-
-//           edge_number_id_map[descriptor.edge_number_vector[edge_id]]       = edge_id;
-//           edge_number_id_map[descriptor.edge_number_vector[other_edge_id]] = other_edge_id;
-//         } else {
-//           edge_number_id_map.erase(descriptor.edge_number_vector[edge_id]);
-//           descriptor.edge_number_vector[edge_id]                     = m_mesh_data.__edges_number[e];
-//           edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_id;
-//         }
-//       }
-//     }
-
-//     for (const auto& ref_edge_list : ref_edges_map) {
-//       Array<EdgeId> edge_list(ref_edge_list.second.size());
-//       for (size_t j = 0; j < ref_edge_list.second.size(); ++j) {
-//         edge_list[j] = ref_edge_list.second[j];
-//       }
-//       const PhysicalRefId& physical_ref_id = m_mesh_data.m_physical_ref_map.at(ref_edge_list.first);
-//       descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list});
-//     }
-//   }
-
-//   std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
-//   for (unsigned int r = 0; r < m_mesh_data.__points.size(); ++r) {
-//     const unsigned int point_number = m_mesh_data.__points[r];
-//     const unsigned int& ref         = m_mesh_data.__points_ref[r];
-//     ref_points_map[ref].push_back(point_number);
-//   }
-
-//   for (const auto& ref_point_list : ref_points_map) {
-//     Array<NodeId> point_list(ref_point_list.second.size());
-//     for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
-//       point_list[j] = ref_point_list.second[j];
-//     }
-//     const PhysicalRefId& physical_ref_id = m_mesh_data.m_physical_ref_map.at(ref_point_list.first);
-//     descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list));
-//   }
-
-//   descriptor.cell_owner_vector.resize(nb_cells);
-//   std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
-
-//   descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
-//   std::fill(descriptor.face_owner_vector.begin(), descriptor.face_owner_vector.end(), parallel::rank());//
-//   descriptor.edge_owner_vector.resize(descriptor.edge_number_vector.size());
-//   std::fill(descriptor.edge_owner_vector.begin(), descriptor.edge_owner_vector.end(), parallel::rank());
-
-//   descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
-//   std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
-
-//   return Connectivity3D::build(descriptor);
-// }
-
 void
 GmshReader::__proceedData()
 {
diff --git a/src/mesh/LogicalConnectivityBuilder.cpp b/src/mesh/LogicalConnectivityBuilder.cpp
index 319e29f6da89f429b3e0ecc657d929497a0f7832..1fd5a950a304a3aec89a10ad769407c4987abd6a 100644
--- a/src/mesh/LogicalConnectivityBuilder.cpp
+++ b/src/mesh/LogicalConnectivityBuilder.cpp
@@ -6,8 +6,6 @@
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
 
-#include <utils/Timer.hpp>
-
 #include <unordered_map>
 
 template <size_t Dimension>
@@ -143,28 +141,25 @@ void
 LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
                                                    ConnectivityDescriptor& descriptor)
 {
-  using Edge                                                                 = ConnectivityFace<2>;
-  const auto& node_number_vector                                             = descriptor.node_number_vector;
-  const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map = [&] {
-    std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
-    for (EdgeId l = 0; l < descriptor.edge_to_node_vector.size(); ++l) {
-      const auto& node_vector                               = descriptor.edge_to_node_vector[l];
-      edge_to_id_map[Edge(node_vector, node_number_vector)] = l;
+  const auto& node_number_vector = descriptor.node_number_vector;
+
+  const auto find_edge = [&](uint32_t node0, uint32_t node1) {
+    if (node_number_vector[node0] > node_number_vector[node1]) {
+      std::swap(node0, node1);
     }
-    return edge_to_id_map;
-  }();
+    const auto& edge_id_vector = descriptor.node_to_edge_matrix[node0];
 
-  std::unordered_map<int, EdgeId> edge_number_id_map = [&] {
-    std::unordered_map<int, EdgeId> edge_number_id_map;
-    for (size_t l = 0; l < descriptor.edge_number_vector.size(); ++l) {
-      edge_number_id_map[descriptor.edge_number_vector[l]] = l;
+    for (size_t i_edge = 0; i_edge < edge_id_vector.size(); ++i_edge) {
+      const EdgeId edge_id = edge_id_vector[i_edge];
+      if (descriptor.edge_to_node_matrix[edge_id][1] == node1) {
+        return edge_id;
+      }
     }
-    Assert(edge_number_id_map.size() == descriptor.edge_number_vector.size());
-    return edge_number_id_map;
-  }();
+    throw UnexpectedError("Cannot find edge");
+  };
 
   const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
-  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+  const auto get_node_id = [&](const TinyVector<3, uint64_t>& node_logic_id) {
     return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
   };
 
@@ -174,13 +169,10 @@ LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>
       Array<EdgeId> boundary_edges(cell_size[2]);
       size_t l = 0;
       for (size_t k = 0; k < cell_size[2]; ++k) {
-        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+        const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+        const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i, j, k + 1});
 
-        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
-        Assert(i_edge != edge_to_id_map.end());
-
-        boundary_edges[l++] = i_edge->second;
+        boundary_edges[l++] = find_edge(node_0_id, node_1_id);
       }
       Assert(l == cell_size[2]);
       descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges, true});
@@ -198,13 +190,10 @@ LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>
       Array<EdgeId> boundary_edges(cell_size[1]);
       size_t l = 0;
       for (size_t j = 0; j < cell_size[1]; ++j) {
-        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
-
-        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
-        Assert(i_edge != edge_to_id_map.end());
+        const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+        const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k});
 
-        boundary_edges[l++] = i_edge->second;
+        boundary_edges[l++] = find_edge(node_0_id, node_1_id);
       }
       Assert(l == cell_size[1]);
       descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges, true});
@@ -222,13 +211,10 @@ LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>
       Array<EdgeId> boundary_edges(cell_size[0]);
       size_t l = 0;
       for (size_t i = 0; i < cell_size[0]; ++i) {
-        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
-
-        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
-        Assert(i_edge != edge_to_id_map.end());
+        const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+        const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k});
 
-        boundary_edges[l++] = i_edge->second;
+        boundary_edges[l++] = find_edge(node_0_id, node_1_id);
       }
       Assert(l == cell_size[0]);
       descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges, true});
@@ -265,7 +251,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
       constexpr size_t left_face = 3;
 
       const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][left_face];
+      const size_t face_id = descriptor.cell_to_face_matrix[cell_id][left_face];
 
       boundary_faces[j] = face_id;
     }
@@ -279,7 +265,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
       constexpr size_t right_face = 1;
 
       const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][right_face];
+      const size_t face_id = descriptor.cell_to_face_matrix[cell_id][right_face];
 
       boundary_faces[j] = face_id;
     }
@@ -293,7 +279,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
       constexpr size_t bottom_face = 0;
 
       const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][bottom_face];
+      const size_t face_id = descriptor.cell_to_face_matrix[cell_id][bottom_face];
 
       boundary_faces[i] = face_id;
     }
@@ -307,7 +293,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
       constexpr size_t top_face = 2;
 
       const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][top_face];
+      const size_t face_id = descriptor.cell_to_face_matrix[cell_id][top_face];
 
       boundary_faces[i] = face_id;
     }
@@ -320,28 +306,51 @@ void
 LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
                                                    ConnectivityDescriptor& descriptor)
 {
-  using Face = ConnectivityFace<3>;
+  const auto& node_number_vector = descriptor.node_number_vector;
 
-  const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
-    std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
-    for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-      const auto& node_vector                                          = descriptor.face_to_node_vector[l];
-      face_to_id_map[Face(node_vector, descriptor.node_number_vector)] = l;
+  const auto find_face = [&](std::array<uint32_t, 4> node_list) {
+    size_t i_node_smallest_number = 0;
+    for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+      if (node_number_vector[node_list[i_node]] < node_number_vector[node_list[i_node_smallest_number]]) {
+        i_node_smallest_number = i_node;
+      }
     }
-    return face_to_id_map;
-  }();
 
-  const std::unordered_map<int, FaceId> face_number_id_map = [&] {
-    std::unordered_map<int, FaceId> face_number_id_map;
-    for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
-      face_number_id_map[descriptor.face_number_vector[l]] = l;
+    if (i_node_smallest_number != 0) {
+      std::array<uint64_t, 4> buffer;
+      for (size_t i_node = i_node_smallest_number; i_node < buffer.size(); ++i_node) {
+        buffer[i_node - i_node_smallest_number] = node_list[i_node];
+      }
+      for (size_t i_node = 0; i_node < i_node_smallest_number; ++i_node) {
+        buffer[i_node + node_list.size() - i_node_smallest_number] = node_list[i_node];
+      }
+
+      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+        node_list[i_node] = buffer[i_node];
+      }
     }
-    Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
-    return face_number_id_map;
-  }();
+
+    if (node_number_vector[node_list[1]] > node_number_vector[node_list[node_list.size() - 1]]) {
+      for (size_t i_node = 1; i_node <= (node_list.size() + 1) / 2 - 1; ++i_node) {
+        std::swap(node_list[i_node], node_list[node_list.size() - i_node]);
+      }
+    }
+
+    const auto& face_id_vector = descriptor.node_to_face_matrix[node_list[0]];
+
+    for (size_t i_face = 0; i_face < face_id_vector.size(); ++i_face) {
+      const FaceId face_id = face_id_vector[i_face];
+      if ((descriptor.face_to_node_matrix[face_id][1] == node_list[1]) and
+          (descriptor.face_to_node_matrix[face_id][2] == node_list[2]) and
+          (descriptor.face_to_node_matrix[face_id][3] == node_list[3])) {
+        return face_id;
+      }
+    }
+    throw UnexpectedError("Cannot find edge");
+  };
 
   const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
-  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+  const auto get_node_id = [&](const TinyVector<3, uint64_t>& node_logic_id) {
     return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
   };
 
@@ -351,16 +360,12 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>
       size_t l = 0;
       for (size_t j = 0; j < cell_size[1]; ++j) {
         for (size_t k = 0; k < cell_size[2]; ++k) {
-          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
-          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k + 1});
-          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
-
-          auto i_face =
-            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
-          Assert(i_face != face_to_id_map.end());
+          const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+          const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k});
+          const uint32_t node_2_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k + 1});
+          const uint32_t node_3_id = get_node_id(TinyVector<3, uint64_t>{i, j, k + 1});
 
-          boundary_faces[l++] = i_face->second;
+          boundary_faces[l++] = find_face({node_0_id, node_1_id, node_2_id, node_3_id});
         }
       }
       Assert(l == cell_size[1] * cell_size[2]);
@@ -377,16 +382,12 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>
       size_t l = 0;
       for (size_t i = 0; i < cell_size[0]; ++i) {
         for (size_t k = 0; k < cell_size[2]; ++k) {
-          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
-          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k + 1});
-          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
-
-          auto i_face =
-            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
-          Assert(i_face != face_to_id_map.end());
+          const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+          const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k});
+          const uint32_t node_2_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k + 1});
+          const uint32_t node_3_id = get_node_id(TinyVector<3, uint64_t>{i, j, k + 1});
 
-          boundary_faces[l++] = i_face->second;
+          boundary_faces[l++] = find_face({node_0_id, node_1_id, node_2_id, node_3_id});
         }
       }
       Assert(l == cell_size[0] * cell_size[2]);
@@ -403,16 +404,12 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>
       size_t l = 0;
       for (size_t i = 0; i < cell_size[0]; ++i) {
         for (size_t j = 0; j < cell_size[1]; ++j) {
-          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
-          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k});
-          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+          const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+          const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k});
+          const uint32_t node_2_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j + 1, k});
+          const uint32_t node_3_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k});
 
-          auto i_face =
-            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
-          Assert(i_face != face_to_id_map.end());
-
-          boundary_faces[l++] = i_face->second;
+          boundary_faces[l++] = find_face({node_0_id, node_1_id, node_2_id, node_3_id});
         }
       }
       Assert(l == cell_size[0] * cell_size[1]);
@@ -452,15 +449,17 @@ LogicalConnectivityBuilder::_buildConnectivity(
   descriptor.cell_type_vector = Array<CellType>(number_of_cells);
   descriptor.cell_type_vector.fill(CellType::Line);
 
-  descriptor.cell_to_node_vector.resize(number_of_cells);
-  constexpr size_t nb_node_per_cell = 2;
-  for (size_t j = 0; j < number_of_cells; ++j) {
-    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
-    for (size_t r = 0; r < nb_node_per_cell; ++r) {
-      descriptor.cell_to_node_vector[j][0] = j;
-      descriptor.cell_to_node_vector[j][1] = j + 1;
-    }
+  Array<unsigned int> cell_to_node_row_map(number_of_cells + 1);
+  for (size_t i = 0; i < cell_to_node_row_map.size(); ++i) {
+    cell_to_node_row_map[i] = 2 * i;
   }
+  Array<unsigned int> cell_to_node_list(2 * number_of_cells);
+  for (size_t i = 0; i < number_of_cells; ++i) {
+    cell_to_node_list[2 * i]     = i;
+    cell_to_node_list[2 * i + 1] = i + 1;
+  }
+
+  descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row_map, cell_to_node_list);
 
   this->_buildBoundaryNodeList(cell_size, descriptor);
 
@@ -515,19 +514,22 @@ LogicalConnectivityBuilder::_buildConnectivity(
     return TinyVector<Dimension, uint64_t>{j0, j1};
   };
 
-  descriptor.cell_to_node_vector.resize(number_of_cells);
-  constexpr size_t nb_node_per_cell = 1 << Dimension;
+  Array<unsigned int> cell_to_node_row_map(number_of_cells + 1);
+  for (size_t i = 0; i < cell_to_node_row_map.size(); ++i) {
+    cell_to_node_row_map[i] = 4 * i;
+  }
+  Array<unsigned int> cell_to_node_list(4 * number_of_cells);
   for (size_t j = 0; j < number_of_cells; ++j) {
     TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
-    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
-    for (size_t r = 0; r < nb_node_per_cell; ++r) {
-      descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0});
-      descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0});
-      descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1});
-      descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1});
-    }
+
+    cell_to_node_list[4 * j]     = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0});
+    cell_to_node_list[4 * j + 1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0});
+    cell_to_node_list[4 * j + 2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1});
+    cell_to_node_list[4 * j + 3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1});
   }
 
+  descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row_map, cell_to_node_list);
+
   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
 
   this->_buildBoundaryNodeList(cell_size, descriptor);
@@ -549,9 +551,6 @@ template <>
 void
 LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& cell_size)
 {
-  Timer t;
-  t.start();
-
   constexpr size_t Dimension = 3;
 
   const TinyVector<Dimension, uint64_t> node_size = [&] {
@@ -609,57 +608,33 @@ LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& ce
     return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
   };
 
-  std::cout << "LogicalConnectivityBuilder:\n - prepare t = " << t << '\n';
-  t.stop();
-  t.start();
-
-  descriptor.cell_to_node_vector.resize(number_of_cells);
-  constexpr size_t nb_node_per_cell = 1 << Dimension;
+  Array<unsigned int> cell_to_node_row_map(number_of_cells + 1);
+  for (size_t i = 0; i < cell_to_node_row_map.size(); ++i) {
+    cell_to_node_row_map[i] = 8 * i;
+  }
+  Array<unsigned int> cell_to_node_list(8 * number_of_cells);
   for (size_t j = 0; j < number_of_cells; ++j) {
     TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
-    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
-    for (size_t r = 0; r < nb_node_per_cell; ++r) {
-      static_assert(Dimension == 3, "unexpected dimension");
-      descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 0});
-      descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 0});
-      descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 0});
-      descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 0});
-      descriptor.cell_to_node_vector[j][4] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 1});
-      descriptor.cell_to_node_vector[j][5] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 1});
-      descriptor.cell_to_node_vector[j][6] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 1});
-      descriptor.cell_to_node_vector[j][7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
-    }
+
+    cell_to_node_list[8 * j]     = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 0});
+    cell_to_node_list[8 * j + 1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 0});
+    cell_to_node_list[8 * j + 2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 0});
+    cell_to_node_list[8 * j + 3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 0});
+    cell_to_node_list[8 * j + 4] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 1});
+    cell_to_node_list[8 * j + 5] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 1});
+    cell_to_node_list[8 * j + 6] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 1});
+    cell_to_node_list[8 * j + 7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
   }
-  std::cout << " - descriptor.cell_to_node_vector t = " << t << '\n';
-  t.stop();
 
-  t.start();
-  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
-  std::cout << " - cell->face, face->node connectivities t = " << t << '\n';
-  t.stop();
+  descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row_map, cell_to_node_list);
 
-  t.start();
+  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
   ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
-  std::cout << " - face->edge, edge->node, cell->edge connectivities t = " << t << '\n';
-  t.stop();
-
-  t.start();
 
   this->_buildBoundaryNodeList(cell_size, descriptor);
-  std::cout << " - boundaryNodeList t = " << t << '\n';
-  t.stop();
-
-  t.start();
   this->_buildBoundaryEdgeList(cell_size, descriptor);
-  std::cout << " - boundaryEdgeList t = " << t << '\n';
-  t.stop();
-
-  t.start();
   this->_buildBoundaryFaceList(cell_size, descriptor);
-  std::cout << " - boundaryFaceList t = " << t << '\n';
-  t.stop();
 
-  t.start();
   descriptor.cell_owner_vector = Array<int>(number_of_cells);
   descriptor.cell_owner_vector.fill(parallel::rank());
 
@@ -671,12 +646,8 @@ LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& ce
 
   descriptor.node_owner_vector = Array<int>(descriptor.node_number_vector.size());
   descriptor.node_owner_vector.fill(parallel::rank());
-  std::cout << " - final infos t = " << t << '\n';
-  t.stop();
 
-  t.start();
   m_connectivity = Connectivity<Dimension>::build(descriptor);
-  std::cout << " - Connectivity<Dimension>::build(descriptor) t = " << t << '\n';
 }
 
 template <size_t Dimension>
diff --git a/src/mesh/MedianDualConnectivityBuilder.cpp b/src/mesh/MedianDualConnectivityBuilder.cpp
index 9c8762b606c3e79b4d86f72139a775c5c197fbf7..fa3490b07b520d93ff33a6f8adeaac467ecac3b2 100644
--- a/src/mesh/MedianDualConnectivityBuilder.cpp
+++ b/src/mesh/MedianDualConnectivityBuilder.cpp
@@ -133,7 +133,6 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
     }
   });
 
-  dual_descriptor.cell_to_node_vector.resize(dual_number_of_cells);
   const auto& primal_cell_to_face_matrix               = primal_connectivity.cellToFaceMatrix();
   const auto& primal_node_to_face_matrix               = primal_connectivity.nodeToFaceMatrix();
   const auto& primal_face_to_cell_matrix               = primal_connectivity.faceToCellMatrix();
@@ -170,16 +169,35 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
     // LCOV_EXCL_STOP
   };
 
+  Array<unsigned int> cell_to_node_row(dual_number_of_cells + 1);
+  cell_to_node_row[0] = 0;
+  {
+    for (NodeId node_id = 0; node_id < primal_number_of_nodes; ++node_id) {
+      // const size_t i_dual_cell             = node_id;
+      const auto& primal_node_to_cell_list = primal_node_to_cell_matrix[node_id];
+      const auto& primal_node_to_face_list = primal_node_to_face_matrix[node_id];
+
+      if (primal_node_to_cell_list.size() != primal_node_to_face_list.size()) {
+        // boundary cell
+        cell_to_node_row[node_id + 1] =
+          cell_to_node_row[node_id] + 1 + primal_node_to_cell_list.size() + primal_node_to_face_list.size();
+      } else {
+        // inner cell
+        cell_to_node_row[node_id + 1] =
+          cell_to_node_row[node_id] + primal_node_to_cell_list.size() + primal_node_to_face_list.size();
+      }
+    }
+  }
+
+  Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
   parallel_for(primal_number_of_nodes, [&](NodeId node_id) {
-    const size_t i_dual_cell             = node_id;
     const auto& primal_node_to_cell_list = primal_node_to_cell_matrix[node_id];
     const auto& primal_node_to_face_list = primal_node_to_face_matrix[node_id];
 
-    auto& dual_cell_node_list = dual_descriptor.cell_to_node_vector[i_dual_cell];
+    size_t i_dual_cell_node = cell_to_node_row[node_id];
 
     if (primal_node_to_cell_list.size() != primal_node_to_face_list.size()) {
       // boundary cell
-      dual_cell_node_list.reserve(1 + primal_node_to_cell_list.size() + primal_node_to_face_list.size());
 
       auto [face_id, cell_id] = [&]() -> std::pair<FaceId, CellId> {
         for (size_t i_face = 0; i_face < primal_node_to_face_list.size(); ++i_face) {
@@ -200,24 +218,24 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
         // LCOV_EXCL_STOP
       }();
 
-      dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
-      dual_cell_node_list.push_back(m_primal_cell_to_dual_node_map[cell_id].second);
+      cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+      cell_to_node_list[i_dual_cell_node++] = m_primal_cell_to_dual_node_map[cell_id].second;
 
       face_id = next_face(cell_id, face_id, node_id);
 
       while (primal_face_to_cell_matrix[face_id].size() > 1) {
-        dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
-        cell_id = next_cell(cell_id, face_id);
-        dual_cell_node_list.push_back(m_primal_cell_to_dual_node_map[cell_id].second);
+        cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+
+        cell_id                               = next_cell(cell_id, face_id);
+        cell_to_node_list[i_dual_cell_node++] = m_primal_cell_to_dual_node_map[cell_id].second;
+
         face_id = next_face(cell_id, face_id, node_id);
       }
-      dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
-      dual_cell_node_list.push_back(node_to_dual_node_correpondance[node_id]);
+      cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+      cell_to_node_list[i_dual_cell_node++] = node_to_dual_node_correpondance[node_id];
 
-      Assert(dual_cell_node_list.size() == 1 + primal_node_to_cell_list.size() + primal_node_to_face_list.size());
     } else {
       // inner cell
-      dual_cell_node_list.reserve(primal_node_to_cell_list.size() + primal_node_to_face_list.size());
 
       auto [face_id, cell_id] = [&]() -> std::pair<FaceId, CellId> {
         const FaceId face_id = primal_node_to_face_list[0];
@@ -237,14 +255,16 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
 
       const FaceId first_face_id = face_id;
       do {
-        dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
-        dual_cell_node_list.push_back(m_primal_cell_to_dual_node_map[cell_id].second);
+        cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+        cell_to_node_list[i_dual_cell_node++] = m_primal_cell_to_dual_node_map[cell_id].second;
 
         face_id = next_face(cell_id, face_id, node_id);
         cell_id = next_cell(cell_id, face_id);
       } while (face_id != first_face_id);
     }
   });
+
+  dual_descriptor.cell_to_node_matrix = ConnectivityMatrix(cell_to_node_row, cell_to_node_list);
 }
 
 template <>
@@ -299,18 +319,6 @@ MedianDualConnectivityBuilder::_buildConnectivityFrom<2>(const IConnectivity& i_
     }
   }
 
-  using Face = ConnectivityFace<2>;
-
-  const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
-    std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
-    for (FaceId l = 0; l < dual_descriptor.face_to_node_vector.size(); ++l) {
-      const auto& node_vector = dual_descriptor.face_to_node_vector[l];
-
-      face_to_id_map[Face(node_vector, dual_descriptor.node_number_vector)] = l;
-    }
-    return face_to_id_map;
-  }();
-
   for (size_t i_face_list = 0; i_face_list < primal_connectivity.template numberOfRefItemList<ItemType::face>();
        ++i_face_list) {
     const auto& primal_ref_face_list = primal_connectivity.template refItemList<ItemType::face>(i_face_list);
@@ -334,9 +342,9 @@ MedianDualConnectivityBuilder::_buildConnectivityFrom<2>(const IConnectivity& i_
 
       std::vector<FaceId> dual_face_list;
       dual_face_list.reserve(2 * primal_face_list.size());
-      for (size_t i_dual_face = 0; i_dual_face < dual_descriptor.face_to_node_vector.size(); ++i_dual_face) {
-        const NodeId dual_node_0 = dual_descriptor.face_to_node_vector[i_dual_face][0];
-        const NodeId dual_node_1 = dual_descriptor.face_to_node_vector[i_dual_face][1];
+      for (size_t i_dual_face = 0; i_dual_face < dual_descriptor.face_to_node_matrix.numberOfRows(); ++i_dual_face) {
+        const NodeId dual_node_0 = dual_descriptor.face_to_node_matrix[i_dual_face][0];
+        const NodeId dual_node_1 = dual_descriptor.face_to_node_matrix[i_dual_face][1];
 
         if ((is_dual_node_from_boundary_face[dual_node_0] and is_dual_node_from_boundary_node[dual_node_1]) or
             (is_dual_node_from_boundary_node[dual_node_0] and is_dual_node_from_boundary_face[dual_node_1])) {
@@ -389,8 +397,8 @@ MedianDualConnectivityBuilder::_buildConnectivityFrom<2>(const IConnectivity& i_
   {
     dual_descriptor.face_owner_vector = Array<int>(dual_descriptor.face_number_vector.size());
     dual_descriptor.face_owner_vector.fill(parallel::size());
-    for (size_t i_cell = 0; i_cell < dual_descriptor.cell_to_face_vector.size(); ++i_cell) {
-      const auto& cell_face_list = dual_descriptor.cell_to_face_vector[i_cell];
+    for (size_t i_cell = 0; i_cell < dual_descriptor.cell_to_face_matrix.numberOfRows(); ++i_cell) {
+      const auto& cell_face_list = dual_descriptor.cell_to_face_matrix[i_cell];
       for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
         const size_t face_id = cell_face_list[i_face];
         dual_descriptor.face_owner_vector[face_id] =
diff --git a/tests/test_Connectivity.cpp b/tests/test_Connectivity.cpp
index 60e894a83a367b6d07ef80939ac0860095dfd322..b0465d08489c644685500d59b145e12594966522 100644
--- a/tests/test_Connectivity.cpp
+++ b/tests/test_Connectivity.cpp
@@ -3,6 +3,7 @@
 
 #include <MeshDataBaseForTests.hpp>
 #include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
 #include <mesh/ItemValue.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
@@ -670,6 +671,318 @@ TEST_CASE("Connectivity", "[mesh]")
     }
   }
 
+  SECTION("item ordering")
+  {
+    SECTION("1D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          SECTION("cell -> nodes")
+          {
+            auto mesh = named_mesh.mesh();
+            auto xr   = mesh->xr();
+
+            const Connectivity<1>& connectivity = mesh->connectivity();
+
+            auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+            bool is_correct = true;
+
+            for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+              if (xr[cell_to_node_matrix[cell_id][1]][0] < xr[cell_to_node_matrix[cell_id][0]][0]) {
+                is_correct = false;
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("node -> cells")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<1>& connectivity = mesh->connectivity();
+
+            auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+            auto cell_number         = connectivity.cellNumber();
+
+            bool is_correct = true;
+            for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+              auto cell_node_list = node_to_cell_matrix[node_id];
+              for (size_t i_node = 0; i_node < cell_node_list.size() - 1; ++i_node) {
+                is_correct &= (cell_number[cell_node_list[i_node]] < cell_number[cell_node_list[i_node + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          SECTION("face -> nodes")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<2>& connectivity = mesh->connectivity();
+
+            auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+            auto node_number         = connectivity.nodeNumber();
+
+            bool is_correct = true;
+
+            for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+              auto face_node_list = face_to_node_matrix[face_id];
+              if (node_number[face_node_list[1]] < node_number[face_node_list[0]]) {
+                is_correct = false;
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("node -> faces")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<2>& connectivity = mesh->connectivity();
+
+            auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
+            auto face_number         = connectivity.faceNumber();
+
+            bool is_correct = true;
+            for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+              auto face_node_list = node_to_face_matrix[node_id];
+              for (size_t i_node = 0; i_node < face_node_list.size() - 1; ++i_node) {
+                is_correct &= (face_number[face_node_list[i_node]] < face_number[face_node_list[i_node + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("node -> cells")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<2>& connectivity = mesh->connectivity();
+
+            auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+            auto cell_number         = connectivity.cellNumber();
+
+            bool is_correct = true;
+            for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+              auto cell_node_list = node_to_cell_matrix[node_id];
+              for (size_t i_node = 0; i_node < cell_node_list.size() - 1; ++i_node) {
+                is_correct &= (cell_number[cell_node_list[i_node]] < cell_number[cell_node_list[i_node + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("face -> cells")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<2>& connectivity = mesh->connectivity();
+
+            auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+            auto cell_number         = connectivity.cellNumber();
+
+            bool is_correct = true;
+            for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+              auto cell_face_list = face_to_cell_matrix[face_id];
+              for (size_t i_face = 0; i_face < cell_face_list.size() - 1; ++i_face) {
+                is_correct &= (cell_number[cell_face_list[i_face]] < cell_number[cell_face_list[i_face + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (const auto& named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          SECTION("edge -> nodes")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+            auto node_number         = connectivity.nodeNumber();
+
+            bool is_correct = true;
+
+            for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+              auto edge_node_list = edge_to_node_matrix[edge_id];
+              if (node_number[edge_node_list[1]] < node_number[edge_node_list[0]]) {
+                is_correct = false;
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("face -> nodes")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+            auto node_number         = connectivity.nodeNumber();
+
+            bool is_correct = true;
+
+            for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+              auto face_node_list = face_to_node_matrix[face_id];
+              for (size_t i = 1; i < face_node_list.size() - 1; ++i) {
+                if (node_number[face_node_list[i]] < node_number[face_node_list[0]]) {
+                  is_correct = false;
+                }
+              }
+              for (size_t i = 2; i < face_node_list.size() - 1; ++i) {
+                if (node_number[face_node_list[i]] < node_number[face_node_list[1]]) {
+                  is_correct = false;
+                }
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("node -> edges")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            REQUIRE(checkConnectivityOrdering(connectivity));
+
+            auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
+            auto edge_number         = connectivity.edgeNumber();
+
+            bool is_correct = true;
+            for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+              auto edge_node_list = node_to_edge_matrix[node_id];
+              for (size_t i_node = 0; i_node < edge_node_list.size() - 1; ++i_node) {
+                is_correct &= (edge_number[edge_node_list[i_node]] < edge_number[edge_node_list[i_node + 1]]);
+              }
+            }
+
+            REQUIRE(is_correct);
+          }
+
+          SECTION("node -> faces")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
+            auto face_number         = connectivity.faceNumber();
+
+            bool is_correct = true;
+            for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+              auto face_node_list = node_to_face_matrix[node_id];
+              for (size_t i_node = 0; i_node < face_node_list.size() - 1; ++i_node) {
+                is_correct &= (face_number[face_node_list[i_node]] < face_number[face_node_list[i_node + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("node -> cells")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+            auto cell_number         = connectivity.cellNumber();
+
+            bool is_correct = true;
+            for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+              auto cell_node_list = node_to_cell_matrix[node_id];
+              for (size_t i_node = 0; i_node < cell_node_list.size() - 1; ++i_node) {
+                is_correct &= (cell_number[cell_node_list[i_node]] < cell_number[cell_node_list[i_node + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("edge -> faces")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            auto edge_to_face_matrix = connectivity.edgeToFaceMatrix();
+            auto face_number         = connectivity.faceNumber();
+
+            bool is_correct = true;
+            for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+              auto face_edge_list = edge_to_face_matrix[edge_id];
+              for (size_t i_edge = 0; i_edge < face_edge_list.size() - 1; ++i_edge) {
+                is_correct &= (face_number[face_edge_list[i_edge]] < face_number[face_edge_list[i_edge + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("edge -> cells")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
+            auto cell_number         = connectivity.cellNumber();
+
+            bool is_correct = true;
+            for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+              auto cell_edge_list = edge_to_cell_matrix[edge_id];
+              for (size_t i_edge = 0; i_edge < cell_edge_list.size() - 1; ++i_edge) {
+                is_correct &= (cell_number[cell_edge_list[i_edge]] < cell_number[cell_edge_list[i_edge + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+
+          SECTION("face -> cells")
+          {
+            auto mesh = named_mesh.mesh();
+
+            const Connectivity<3>& connectivity = mesh->connectivity();
+
+            auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+            auto cell_number         = connectivity.cellNumber();
+
+            bool is_correct = true;
+            for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+              auto cell_face_list = face_to_cell_matrix[face_id];
+              for (size_t i_face = 0; i_face < cell_face_list.size() - 1; ++i_face) {
+                is_correct &= (cell_number[cell_face_list[i_face]] < cell_number[cell_face_list[i_face + 1]]);
+              }
+            }
+            REQUIRE(is_correct);
+          }
+        }
+      }
+    }
+  }
+
   SECTION("ItemLocalNumbersInTheirSubItems")
   {
     auto check_item_local_numbers_in_their_subitems = [](auto item_to_subitem_matrix, auto subitem_to_item_matrix,
diff --git a/tests/test_DiamondDualConnectivityBuilder.cpp b/tests/test_DiamondDualConnectivityBuilder.cpp
index 70b8f75440061530f386654a4bbf9e3f5342d33e..45d9cd700ca7b5ae73e6c6fd22db024e7de3126c 100644
--- a/tests/test_DiamondDualConnectivityBuilder.cpp
+++ b/tests/test_DiamondDualConnectivityBuilder.cpp
@@ -5,6 +5,7 @@
 #include <mesh/DualConnectivityManager.hpp>
 
 #include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 
@@ -54,6 +55,8 @@ TEST_CASE("DiamondDualConnectivityBuilder", "[mesh]")
     REQUIRE(dual_connectivity.numberOfFaces() == 220);
     REQUIRE(dual_connectivity.numberOfCells() == 110);
 
+    REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
     SECTION("ref node list")
     {
       REQUIRE(primal_connectivity.numberOfRefItemList<ItemType::node>() == 4);
@@ -159,6 +162,8 @@ TEST_CASE("DiamondDualConnectivityBuilder", "[mesh]")
       DualConnectivityManager::instance().getDiamondDualConnectivity(primal_connectivity);
     const ConnectivityType& dual_connectivity = *p_diamond_dual_connectivity;
 
+    REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
     REQUIRE(dual_connectivity.numberOfNodes() == 331);
     REQUIRE(dual_connectivity.numberOfEdges() == 1461);
     REQUIRE(dual_connectivity.numberOfFaces() == 1651);
diff --git a/tests/test_Dual1DConnectivityBuilder.cpp b/tests/test_Dual1DConnectivityBuilder.cpp
index 2b194a92baa821d6da4bac1e6277aaea6c1ddb1b..b7a35abfc577d93f729a0befc166c16c04f47fd3 100644
--- a/tests/test_Dual1DConnectivityBuilder.cpp
+++ b/tests/test_Dual1DConnectivityBuilder.cpp
@@ -5,6 +5,7 @@
 #include <mesh/DualConnectivityManager.hpp>
 
 #include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 
@@ -48,6 +49,8 @@ TEST_CASE("Dual1DConnectivityBuilder", "[mesh]")
     DualConnectivityManager::instance().getDual1DConnectivity(primal_connectivity);
   const ConnectivityType& dual_connectivity = *p_dual_1d_connectivity;
 
+  REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
   REQUIRE(dual_connectivity.numberOfNodes() == 36);
   REQUIRE(dual_connectivity.numberOfCells() == 35);
 
diff --git a/tests/test_MedianDualConnectivityBuilder.cpp b/tests/test_MedianDualConnectivityBuilder.cpp
index a95e6ad22f8a9b2605dec5c28e7c7704d40aed6f..6d091bd037e4ac57635757b98671eff9fe6eb643 100644
--- a/tests/test_MedianDualConnectivityBuilder.cpp
+++ b/tests/test_MedianDualConnectivityBuilder.cpp
@@ -5,6 +5,7 @@
 #include <mesh/DualConnectivityManager.hpp>
 
 #include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 
@@ -50,6 +51,8 @@ TEST_CASE("MedianDualConnectivityBuilder", "[mesh]")
       DualConnectivityManager::instance().getMedianDualConnectivity(primal_connectivity);
     const ConnectivityType& dual_connectivity = *p_median_dual_connectivity;
 
+    REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
     REQUIRE(dual_connectivity.numberOfNodes() == 192);
     REQUIRE(dual_connectivity.numberOfFaces() == 244);
     REQUIRE(dual_connectivity.numberOfCells() == 53);