diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 391735aaa6e28f87b4910e0d803a123bd9fb28c0..974067b18557dbe64f6736ad718cedaa15a2d676 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -7,9 +7,13 @@ void Connectivity<3>::_computeFaceCellConnectivities()
   Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
 
   typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
+
+  const auto& cell_to_node_matrix
+      = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+
   std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
   for (unsigned int j=0; j<this->numberOfCells(); ++j) {
-    const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
+    const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
 
     switch (cell_nodes.length) {
       case 4: { // tetrahedron
@@ -182,11 +186,14 @@ void Connectivity<3>::_computeFaceCellConnectivities()
 template<>
 void Connectivity<2>::_computeFaceCellConnectivities()
 {
+  const auto& cell_to_node_matrix
+      = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+
   // In 2D faces are simply define
   typedef std::pair<unsigned int, unsigned short> CellFaceId;
   std::map<Face, std::vector<CellFaceId>> face_cells_map;
   for (unsigned int j=0; j<this->numberOfCells(); ++j) {
-    const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
+    const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
     for (unsigned short r=0; r<cell_nodes.length; ++r) {
       unsigned int node0_id = cell_nodes(r);
       unsigned int node1_id = cell_nodes((r+1)%cell_nodes.length);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index d0448615bee9e11b2f0464851977e8f2d257cfa6..931482a5f8b3963535f3ec912649487bb15e2d06 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -230,12 +230,10 @@ class Connectivity final
  public:
   static constexpr size_t dimension = Dimension;
 
- private:
-  ConnectivityMatrix m_cell_to_node_matrix;
  public:
   ConnectivityMatrix cellToNodeMatrix() const
   {
-    return m_cell_to_node_matrix;
+    return m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
   }
 
   NodeValuePerCell<unsigned short> m_cell_to_node_local_cell;
@@ -323,7 +321,9 @@ private:
   KOKKOS_INLINE_FUNCTION
   size_t numberOfCells() const
   {
-    return m_cell_to_node_matrix.numRows();
+    const auto& cell_to_node_matrix
+        = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+    return cell_to_node_matrix.numRows();
   }
 
   const Kokkos::View<const double*> invCellNbNodes() const
@@ -347,32 +347,34 @@ private:
 
   Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
   {
-    m_cell_to_node_matrix = cell_by_node_vector;
+    auto& cell_to_node_matrix
+        = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+    cell_to_node_matrix = cell_by_node_vector;
 
     Assert(this->numberOfCells()>0);
 
     {
       Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
       Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
-          const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
+          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
           inv_cell_nb_nodes[j] = 1./cell_nodes.length;
         });
       m_inv_cell_nb_nodes = inv_cell_nb_nodes;
     }
 
-    m_connectivity_computer.computeInverseConnectivityMatrix(m_cell_to_node_matrix,
+    m_connectivity_computer.computeInverseConnectivityMatrix(cell_to_node_matrix,
                                                              m_node_to_cell_matrix);
 
     m_node_to_cell_local_node = CellValuePerNode<unsigned short>(*this);
 
-    m_connectivity_computer.computeLocalChildItemNumberInItem(m_cell_to_node_matrix,
+    m_connectivity_computer.computeLocalChildItemNumberInItem(cell_to_node_matrix,
                                                               m_node_to_cell_matrix,
                                                               m_node_to_cell_local_node);
 
     m_cell_to_node_local_cell = NodeValuePerCell<unsigned short>(*this);
 
     m_connectivity_computer.computeLocalChildItemNumberInItem(m_node_to_cell_matrix,
-                                                              m_cell_to_node_matrix,
+                                                              cell_to_node_matrix,
                                                               m_cell_to_node_local_cell);
     if constexpr (Dimension>1) {
       this->_computeFaceCellConnectivities();
@@ -403,7 +405,10 @@ inline const ConnectivityMatrix&
 Connectivity<3>::itemToItemMatrix<TypeOfItem::cell,
                                   TypeOfItem::node>() const
 {
-  return m_cell_to_node_matrix;
+  const auto& cell_to_node_matrix
+      = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+
+  return cell_to_node_matrix;
 }
 
 template <>
@@ -450,7 +455,10 @@ inline const ConnectivityMatrix&
 Connectivity<2>::itemToItemMatrix<TypeOfItem::cell,
                                   TypeOfItem::node>() const
 {
-  return m_cell_to_node_matrix;
+  const auto& cell_to_node_matrix
+      = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+
+  return cell_to_node_matrix;
 }
 
 template <>
@@ -488,7 +496,10 @@ inline const ConnectivityMatrix&
 Connectivity<1>::itemToItemMatrix<TypeOfItem::cell,
                                   TypeOfItem::node>() const
 {
-  return m_cell_to_node_matrix;
+  const auto& cell_to_node_matrix
+      = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+
+  return cell_to_node_matrix;
 }
 
 template <>