diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1c4e72c356c1337caf1fcf72d644cdf00d9b42cd..9a397a2753ce5021034007d171cf3d41c74db138 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -275,10 +275,10 @@ add_executable(
 # Libraries
 target_link_libraries(
   pastis
-  kokkos
+  PastisMesh
   PastisUtils
+  PastisLanguage
+  kokkos
   ${PARMETIS_LIBRARIES}
   ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
-  PastisLanguage
-  PastisMesh
-  PastisUtils)
+)
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index faf177c5e19c5e86e1694904a806d963655598d9..2dde6fc1ac9d72a298b4157d0388fe3045103f16 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -296,7 +296,7 @@ template <size_t N, typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<N,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
   static_assert(std::is_floating_point<T>::value, "determinent for arbitrary dimension N is defined for floating point types only");
   TinyMatrix<N,T> M = A;
 
@@ -340,7 +340,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<1,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
   return A(0,0);
 }
 
@@ -348,7 +348,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<2,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
   return A(0,0)*A(1,1)-A(1,0)*A(0,1);
 }
 
@@ -356,7 +356,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr T det(const TinyMatrix<3,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
   return
     A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
     -A(1,0)*(A(0,1)*A(2,2)-A(2,1)*A(0,2))
@@ -399,7 +399,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr TinyMatrix<1,T> inverse(const TinyMatrix<1,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
   static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
 
   TinyMatrix<1,T> A_1(1./A(0,0));
@@ -412,7 +412,7 @@ constexpr T cofactor(const TinyMatrix<N,T>& A,
                      const size_t& i,
                      const size_t& j)
 {
-  static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non-arithmetic types");
   const T sign = ((i+j)%2) ? -1 : 1;
 
   return sign * det(getMinor(A, i, j));
@@ -422,7 +422,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
   static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
 
   const T determinent = det(A);
@@ -438,7 +438,7 @@ template <typename T>
 PASTIS_INLINE
 constexpr TinyMatrix<3,T> inverse(const TinyMatrix<3,T>& A)
 {
-  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
   static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
 
   const T determinent = det(A);
diff --git a/src/main.cpp b/src/main.cpp
index 5f9949819724839153c1bf3233a5250c0f325699..6620cf4c68e3fdd41ec11b665a363d17ef20f9f2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -83,9 +83,9 @@ int main(int argc, char *argv[])
               case BoundaryConditionDescriptor::Type::symmetry: {
                 const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
                     = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-                for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefNodeList();
+                for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefItemList<ItemType::node>();
                      ++i_ref_node_list) {
-                  const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(i_ref_node_list);
+                  const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
                   const RefId& ref = ref_node_list.refId();
                   if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                     SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
@@ -200,9 +200,9 @@ int main(int argc, char *argv[])
               case BoundaryConditionDescriptor::Type::symmetry: {
                 const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
                     = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefItemList<ItemType::face>();
                      ++i_ref_face_list) {
-                  const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+                  const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
                   const RefId& ref = ref_face_list.refId();
                   if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                     SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
@@ -304,9 +304,9 @@ int main(int argc, char *argv[])
               case BoundaryConditionDescriptor::Type::symmetry: {
                 const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
                     = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+                for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefItemList<ItemType::face>();
                      ++i_ref_face_list) {
-                  const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+                  const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
                   const RefId& ref = ref_face_list.refId();
                   if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                     SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 0ea63e1418f8c879a43b59873d6ef16ad5bf526a..3a4b45766d36453e862aaeca2b5d8d855942e495 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -3,6 +3,8 @@
 
 #include <Messenger.hpp>
 
+#include <ConnectivityDescriptor.hpp>
+
 template<size_t Dimension>
 Connectivity<Dimension>::Connectivity() {}
 
@@ -10,8 +12,7 @@ template<size_t Dimension>
 void Connectivity<Dimension>::
 _buildFrom(const ConnectivityDescriptor& descriptor)
 {
-#warning All of these should be checked by ConnectivityDescriptor
-  Assert(descriptor.cell_by_node_vector.size() == descriptor.cell_type_vector.size());
+  Assert(descriptor.cell_to_node_vector.size() == 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());
@@ -21,7 +22,7 @@ _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_by_node_vector;
+  cell_to_node_matrix = descriptor.cell_to_node_vector;
 
   {
     WeakCellValue<CellType> cell_type(*this);
@@ -45,8 +46,6 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
 
   {
     WeakCellValue<int> cell_global_index(*this);
-#warning index must start accounting number of global indices of other procs
-#warning must take care of ghost cells
     int first_index = 0;
     parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
       cell_global_index[j] = first_index+j;
@@ -54,16 +53,6 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
     m_cell_global_index = cell_global_index;
   }
 
-
-  {
-    WeakCellValue<double> inv_cell_nb_nodes(*this);
-    parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& 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;
-  }
-
   {
     WeakCellValue<int> cell_owner(*this);
     cell_owner = convert_to_array(descriptor.cell_owner_vector);
@@ -94,14 +83,13 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
     m_node_is_owned = node_is_owned;
   }
 
+  m_ref_node_list_vector = descriptor.template refItemListVector<ItemType::node>();
+  m_ref_cell_list_vector = descriptor.template refItemListVector<ItemType::cell>();
+
   if constexpr (Dimension>1) {
-    auto& face_to_node_matrix
-        = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
-    face_to_node_matrix = descriptor.face_to_node_vector;
+    m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_vector;
 
-    auto& cell_to_face_matrix
-        = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
-    cell_to_face_matrix = descriptor.cell_to_face_vector;
+    m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_vector;
 
     {
       FaceValuePerCell<bool> cell_face_is_reversed(*this);
@@ -111,14 +99,15 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
           cell_face_is_reversed(j,lj) = face_cells_vector[lj];
         }
       }
-
       m_cell_face_is_reversed = cell_face_is_reversed;
     }
+
     {
       WeakFaceValue<int> face_number(*this);
       face_number = convert_to_array(descriptor.face_number_vector);
       m_face_number = face_number;
     }
+
     {
       WeakFaceValue<int> face_owner(*this);
       face_owner = convert_to_array(descriptor.face_owner_vector);
@@ -134,8 +123,49 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
       m_face_is_owned = face_is_owned;
     }
 
+    m_ref_face_list_vector = descriptor.template refItemListVector<ItemType::face>();
+
+    if constexpr (Dimension>2) {
+      m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_vector;
 
-    m_ref_face_list = descriptor.refFaceList();
+      m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.face_to_edge_vector;
+
+      m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cell_to_edge_vector;
+
+      {
+        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;
+      }
+
+      {
+        WeakEdgeValue<int> edge_number(*this);
+        edge_number = convert_to_array(descriptor.edge_number_vector);
+        m_edge_number = edge_number;
+      }
+
+      {
+        WeakEdgeValue<int> edge_owner(*this);
+        edge_owner = convert_to_array(descriptor.edge_owner_vector);
+        m_edge_owner = edge_owner;
+      }
+
+      {
+        const int rank = parallel::rank();
+        WeakEdgeValue<bool> edge_is_owned(*this);
+        parallel_for(this->numberOfEdges(), PASTIS_LAMBDA(const EdgeId& e) {
+            edge_is_owned[e] = (m_edge_owner[e] == rank);
+          });
+        m_edge_is_owned = edge_is_owned;
+      }
+
+      m_ref_edge_list_vector = descriptor.template refItemListVector<ItemType::edge>();
+    }
   }
 }
 
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 4c19ae6b53a97523f977fdb56a5d310c84227b1d..7e1ebd1198094a1f214e87e26106ea28cfb0a4f1 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -7,6 +7,8 @@
 #include <PastisOStream.hpp>
 #include <PastisUtils.hpp>
 
+#include <PastisTraits.hpp>
+
 #include <TinyVector.hpp>
 
 #include <ItemValue.hpp>
@@ -28,58 +30,14 @@
 
 #include <RefId.hpp>
 #include <ItemType.hpp>
-#include <RefNodeList.hpp>
-#include <RefFaceList.hpp>
+#include <RefItemList.hpp>
+
 
 #include <SynchronizerManager.hpp>
 
-#include <tuple>
-#include <algorithm>
 #include <set>
 
-class ConnectivityDescriptor
-{
- private:
-  std::vector<RefFaceList> m_ref_face_list;
-
- public:
-  std::vector<std::vector<unsigned int>> cell_by_node_vector;
-
-#warning should be set in 2d
-  std::vector<std::vector<unsigned int>> cell_to_face_vector;
-  std::vector<std::vector<bool>> cell_face_is_reversed_vector;
-
-  // std::vector<std::vector<unsigned int>> face_to_cell_vector;
-  std::vector<std::vector<unsigned int>> face_to_node_vector;
-
-  std::vector<CellType> cell_type_vector;
-
-  std::vector<int> cell_number_vector;
-  std::vector<int> face_number_vector;
-  std::vector<int> node_number_vector;
-
-  std::vector<int> cell_owner_vector;
-  std::vector<int> face_owner_vector;
-  std::vector<int> node_owner_vector;
-
-  const std::vector<RefFaceList>& refFaceList() const
-  {
-    return m_ref_face_list;
-  }
-
-  void addRefFaceList(const RefFaceList& ref_face_list)
-  {
-    m_ref_face_list.push_back(ref_face_list);
-  }
-
-  ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
-  ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
-
-  ConnectivityDescriptor() = default;
-  ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
-  ConnectivityDescriptor(ConnectivityDescriptor&&) = delete;
-  ~ConnectivityDescriptor() = default;
-};
+class ConnectivityDescriptor;
 
 template <size_t Dim>
 class Connectivity final
@@ -105,12 +63,10 @@ class Connectivity final
   ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
   WeakCellValue<const CellType> m_cell_type;
 
-#warning is m_cell_global_index really necessary? should it be computed on demand instead?
   WeakCellValue<const int> m_cell_global_index;
 
   WeakCellValue<const int> m_cell_number;
   WeakFaceValue<const int> m_face_number;
-#warning check that m_edge_number is filled correctly
   WeakEdgeValue<const int> m_edge_number;
   WeakNodeValue<const int> m_node_number;
 
@@ -127,6 +83,7 @@ class Connectivity final
   WeakNodeValue<const bool> m_node_is_owned;
 
   WeakFaceValuePerCell<const bool> m_cell_face_is_reversed;
+  WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed;
 
   WeakNodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
   WeakEdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;
@@ -146,8 +103,10 @@ class Connectivity final
 
   ConnectivityComputer m_connectivity_computer;
 
-  std::vector<RefFaceList> m_ref_face_list;
-  std::vector<RefNodeList> m_ref_node_list;
+  std::vector<RefCellList> m_ref_cell_list_vector;
+  std::vector<RefFaceList> m_ref_face_list_vector;
+  std::vector<RefEdgeList> m_ref_edge_list_vector;
+  std::vector<RefNodeList> m_ref_node_list_vector;
 
   WeakCellValue<const double> m_inv_cell_nb_nodes;
 
@@ -228,8 +187,7 @@ class Connectivity final
     } else if constexpr(item_type == ItemType::node) {
       return m_node_number;
     } else {
-      static_assert(item_type == ItemType::cell, "unknown ItemType");
-      return m_cell_number;
+      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
     }
   }
 
@@ -272,8 +230,7 @@ class Connectivity final
     } else if constexpr(item_type == ItemType::node) {
       return m_node_owner;
     } else {
-      static_assert(item_type == ItemType::cell, "unknown ItemType");
-      return m_cell_owner;
+      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
     }
   }
 
@@ -316,8 +273,7 @@ class Connectivity final
     } else if constexpr(item_type == ItemType::node) {
       return m_node_is_owned;
     } else {
-      static_assert(item_type == ItemType::cell, "unknown ItemType");
-      return m_cell_is_owned;
+      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
     }
   }
 
@@ -418,6 +374,13 @@ class Connectivity final
     return m_cell_face_is_reversed;
   }
 
+  PASTIS_INLINE
+  const auto& faceEdgeIsReversed() const
+  {
+    static_assert(Dimension>2, "reversed edges makes no sense in dimension 1 or 2");
+    return m_face_edge_is_reversed;
+  }
+
   PASTIS_INLINE
   const auto& cellLocalNumbersInTheirNodes() const
   {
@@ -512,29 +475,53 @@ class Connectivity final
     return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
   }
 
-  size_t numberOfRefFaceList() const
-  {
-    return m_ref_face_list.size();
-  }
-
-  const RefFaceList& refFaceList(const size_t& i) const
-  {
-    return m_ref_face_list[i];
-  }
+  template <ItemType item_type>
+  size_t numberOfRefItemList() const
+  {
+    if constexpr (item_type == ItemType::cell) {
+      return m_ref_cell_list_vector.size();
+    } else if constexpr (item_type == ItemType::face) {
+      return m_ref_face_list_vector.size();
+    } else if constexpr (item_type == ItemType::edge) {
+      return m_ref_edge_list_vector.size();
+    } else if constexpr (item_type == ItemType::node) {
+      return m_ref_node_list_vector.size();
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
 
-  void addRefNodeList(const RefNodeList& ref_node_list)
-  {
-    m_ref_node_list.push_back(ref_node_list);
   }
 
-  size_t numberOfRefNodeList() const
-  {
-    return m_ref_node_list.size();
+  template <ItemType item_type>
+  const RefItemList<item_type>& refItemList(const size_t& i) const
+  {
+    if constexpr (item_type == ItemType::cell) {
+      return m_ref_cell_list_vector[i];
+    } else if constexpr (item_type == ItemType::face) {
+      return m_ref_face_list_vector[i];
+    } else if constexpr (item_type == ItemType::edge) {
+      return m_ref_edge_list_vector[i];
+    } else if constexpr (item_type == ItemType::node) {
+      return m_ref_node_list_vector[i];
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
   }
 
-  const RefNodeList& refNodeList(const size_t& i) const
-  {
-    return m_ref_node_list[i];
+  template <ItemType item_type>
+  void addRefItemList(const RefItemList<item_type>& ref_item_list)
+  {
+    if constexpr (item_type == ItemType::cell) {
+      m_ref_cell_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::face) {
+      m_ref_face_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::edge) {
+      m_ref_edge_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::node) {
+      m_ref_node_list_vector.push_back(ref_item_list);
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
   }
 
   PASTIS_INLINE
@@ -624,7 +611,18 @@ class Connectivity final
 
   CellValue<const double> invCellNbNodes() const
   {
-#warning add calculation on demand when variables will be defined
+    if (not m_inv_cell_nb_nodes.isBuilt()) {
+      const auto& cell_to_node_matrix
+          = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
+
+      WeakCellValue<double> inv_cell_nb_nodes(*this);
+      parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+          inv_cell_nb_nodes[j] = 1./cell_nodes.length;
+        });
+      const_cast<WeakCellValue<const double>&>(m_inv_cell_nb_nodes) = inv_cell_nb_nodes;
+    }
+
     return m_inv_cell_nb_nodes;
   }
 
diff --git a/src/mesh/ConnectivityDescriptor.hpp b/src/mesh/ConnectivityDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e20b05c11ac585b4293fc233bd0919efce48e5ea
--- /dev/null
+++ b/src/mesh/ConnectivityDescriptor.hpp
@@ -0,0 +1,121 @@
+#ifndef CONNECTIVITY_DESCRIPTOR_HPP
+#define CONNECTIVITY_DESCRIPTOR_HPP
+
+#include <RefItemList.hpp>
+#include <PastisTraits.hpp>
+#include <ItemOfItemType.hpp>
+
+#include <vector>
+
+class ConnectivityDescriptor
+{
+ private:
+  std::vector<RefCellList> m_ref_cell_list_vector;
+  std::vector<RefFaceList> m_ref_face_list_vector;
+  std::vector<RefEdgeList> m_ref_edge_list_vector;
+  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;
+
+  std::vector<std::vector<unsigned int>> face_to_node_vector;
+  std::vector<std::vector<unsigned int>> face_to_edge_vector;
+
+  std::vector<std::vector<unsigned int>> edge_to_node_vector;
+
+  template <typename ItemOfItemT>
+  auto& itemOfItemVector()
+  {
+    if constexpr (std::is_same_v<ItemOfItemT,NodeOfCell>) {
+      return cell_to_node_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,FaceOfCell>) {
+      return cell_to_face_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfCell>) {
+      return cell_to_edge_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfFace>) {
+      return face_to_edge_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,NodeOfFace>) {
+      return face_to_node_vector;
+    } else if constexpr (std::is_same_v<ItemOfItemT,NodeOfEdge>) {
+      return edge_to_node_vector;
+    } else {
+      static_assert(is_false_v<ItemOfItemT>, "Unexpected item of item type");
+    }
+  }
+
+  std::vector<Array<bool>> cell_face_is_reversed_vector;
+  std::vector<Array<bool>> face_edge_is_reversed_vector;
+
+  std::vector<CellType> cell_type_vector;
+
+  std::vector<int> cell_number_vector;
+  std::vector<int> face_number_vector;
+  std::vector<int> edge_number_vector;
+  std::vector<int> node_number_vector;
+
+  template <ItemType item_type>
+  const std::vector<int>& itemNumberVector() const
+  {
+    if constexpr (item_type == ItemType::cell) {
+      return cell_number_vector;
+    } else if constexpr (item_type == ItemType::face) {
+      return face_number_vector;
+    } else if constexpr (item_type == ItemType::edge) {
+      return edge_number_vector;
+    } else if constexpr (item_type == ItemType::node) {
+      return node_number_vector;
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
+  }
+
+  std::vector<int> cell_owner_vector;
+  std::vector<int> face_owner_vector;
+  std::vector<int> edge_owner_vector;
+  std::vector<int> node_owner_vector;
+
+  template <ItemType item_type>
+  const std::vector<RefItemList<item_type>>& refItemListVector() const
+  {
+    if constexpr (item_type == ItemType::cell) {
+      return m_ref_cell_list_vector;
+    } else if constexpr (item_type == ItemType::face) {
+      return m_ref_face_list_vector;
+    } else if constexpr (item_type == ItemType::edge) {
+      return m_ref_edge_list_vector;
+    } else if constexpr (item_type == ItemType::node) {
+      return m_ref_node_list_vector;
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
+  }
+
+  template <ItemType item_type>
+  void addRefItemList(const RefItemList<item_type>& ref_item_list)
+  {
+    if constexpr (item_type == ItemType::cell) {
+      m_ref_cell_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::face) {
+      m_ref_face_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::edge) {
+      m_ref_edge_list_vector.push_back(ref_item_list);
+    } else if constexpr (item_type == ItemType::node) {
+      m_ref_node_list_vector.push_back(ref_item_list);
+    } else {
+      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
+    }
+  }
+
+  ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
+  ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
+
+  ConnectivityDescriptor() = default;
+  ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
+  ConnectivityDescriptor(ConnectivityDescriptor&&) = delete;
+  ~ConnectivityDescriptor() = default;
+};
+
+
+#endif // CONNECTIVITY_DESCRIPTOR_HPP
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index 69d97399d5b0bb5533758f94593f2703ce8bfcb9..aa26cb6f71422cfaf646f3d52ab22ff882e0f72a 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -23,7 +23,6 @@ ConnectivityDispatcher<Dimension>::_buildNewOwner()
         = m_connectivity.template getItemToItemMatrix<item_type,ItemType::cell>();
 
     const auto& cell_number = m_connectivity.cellNumber();
-
     const auto& cell_new_owner = this->_dispatchedInfo<ItemType::cell>().m_new_owner;
 
     using ItemId = ItemIdT<item_type>;
@@ -41,12 +40,24 @@ ConnectivityDispatcher<Dimension>::_buildNewOwner()
         item_new_owner[l] = cell_new_owner[Jmin];
       });
 
-#warning Add missing synchronize (fix it!)
-    // synchronize(face_new_owner);
+    synchronize(item_new_owner);
     this->_dispatchedInfo<item_type>().m_new_owner = item_new_owner;
   }
 }
 
+template <int Dimension>
+template <ItemType item_type>
+void ConnectivityDispatcher<Dimension>::
+_buildItemToExchangeLists()
+{
+  this->_buildItemListToSend<item_type>();
+  this->_buildNumberOfItemToExchange<item_type>();
+  if constexpr (item_type == ItemType::cell) {
+    this->_buildCellNumberIdMap();
+  }
+  this->_buildRecvItemIdCorrespondanceByProc<item_type>();
+}
+
 template <int Dimension>
 template <ItemType item_type>
 void
@@ -117,17 +128,22 @@ ConnectivityDispatcher<Dimension>::_buildItemListToSend()
 }
 
 template <int Dimension>
-Array<const unsigned int>
-ConnectivityDispatcher<Dimension>::_buildNbCellToSend()
+template <ItemType item_type>
+void
+ConnectivityDispatcher<Dimension>::_buildNumberOfItemToExchange()
 {
-  const auto& item_list_to_send_by_proc = this->_dispatchedInfo<ItemType::cell>().m_list_to_send_by_proc;
-  Array<unsigned int> nb_cell_to_send_by_proc(parallel::size());
+  const auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
+  Array<unsigned int> nb_item_to_send_by_proc(parallel::size());
   for (size_t i=0; i<parallel::size(); ++i) {
-    nb_cell_to_send_by_proc[i] = item_list_to_send_by_proc[i].size();
+    nb_item_to_send_by_proc[i] = item_list_to_send_by_proc[i].size();
   }
-  return nb_cell_to_send_by_proc;
+  this->_dispatchedInfo<item_type>().m_list_to_send_size_by_proc = nb_item_to_send_by_proc;
+
+  this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc
+      = parallel::allToAll(nb_item_to_send_by_proc);
 }
 
+
 template <int Dimension>
 template<typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
@@ -150,6 +166,59 @@ _gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gathe
   }
 }
 
+template <int Dimension>
+template<typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+void ConnectivityDispatcher<Dimension>::
+_gatherFrom(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
+            std::vector<Array<std::remove_const_t<DataType>>>& gathered_vector)
+{
+  using MutableDataType = std::remove_const_t<DataType>;
+
+  constexpr ItemType item_type = ItemOfItem::item_type;
+  using ItemId = ItemIdT<item_type>;
+
+  const auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
+
+  std::vector<Array<MutableDataType>> data_to_send_by_proc(parallel::size());
+
+  for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
+    std::vector<MutableDataType> data_by_item_vector;
+    for (size_t j=0; j<item_list_to_send_by_proc[i_rank].size(); ++j) {
+      const ItemId& item_id = item_list_to_send_by_proc[i_rank][j];
+      const auto& item_data = data_to_gather.itemValues(item_id);
+      for (size_t l=0; l<item_data.size(); ++l) {
+        data_by_item_vector.push_back(item_data[l]);
+      }
+    }
+    data_to_send_by_proc[i_rank] = convert_to_array(data_by_item_vector);
+  }
+
+  const auto& number_of_sub_item_per_item_to_recv_by_proc =
+      this->_dispatchedInfo<ItemOfItem>().m_number_of_sub_item_per_item_to_recv_by_proc;
+
+  std::vector<Array<MutableDataType>> recv_data_to_gather_by_proc(parallel::size());
+  for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
+    recv_data_to_gather_by_proc[i_rank]
+        = Array<MutableDataType>(sum(number_of_sub_item_per_item_to_recv_by_proc[i_rank]));
+  }
+
+  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;
+
+  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) {
+      Array<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_vector.emplace_back(data_vector);
+    }
+  }
+}
+
 template <int Dimension>
 void
 ConnectivityDispatcher<Dimension>::
@@ -170,16 +239,18 @@ _buildCellNumberIdMap()
 template <int Dimension>
 template <typename ItemOfItemT>
 void
-ConnectivityDispatcher<Dimension>::
-_buildSubItemNumberToIdMap(const std::vector<Array<const int>>& recv_cell_sub_item_number_by_proc)
+ConnectivityDispatcher<Dimension>::_buildSubItemNumberToIdMap()
 {
-  static_assert(ItemOfItemT::item_type == ItemType::cell, "Dispatcher requires to be build using cell as master entities");
+  static_assert(ItemOfItemT::item_type == ItemType::cell, "Dispatcher requires to be built using cell as master entities");
+
+  const auto& cell_sub_item_number_to_recv_by_proc
+      = this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
 
   auto& sub_item_number_id_map = this->_dispatchedInfo<ItemOfItemT::sub_item_type>().m_number_to_id_map;
   for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
     int sub_item_id=0;
-    for (size_t i=0; i<recv_cell_sub_item_number_by_proc[i_rank].size(); ++i) {
-      int sub_item_number = recv_cell_sub_item_number_by_proc[i_rank][i];
+    for (size_t i=0; i<cell_sub_item_number_to_recv_by_proc[i_rank].size(); ++i) {
+      int sub_item_number = cell_sub_item_number_to_recv_by_proc[i_rank][i];
       auto [iterator, inserted] = sub_item_number_id_map.insert(std::make_pair(sub_item_number, sub_item_id));
       if (inserted) sub_item_id++;
     }
@@ -188,8 +259,8 @@ _buildSubItemNumberToIdMap(const std::vector<Array<const int>>& recv_cell_sub_it
 
 template <int Dimension>
 template <typename SubItemOfItemT>
-std::vector<Array<const int>>
-ConnectivityDispatcher<Dimension>::_getRecvNumberOfSubItemPerItemByProc()
+void
+ConnectivityDispatcher<Dimension>::_buildNumberOfSubItemPerItemToRecvByProc()
 {
   const auto& item_to_sub_item_matrix
       = m_connectivity.template getItemToItemMatrix<SubItemOfItemT::item_type,
@@ -201,58 +272,98 @@ ConnectivityDispatcher<Dimension>::_getRecvNumberOfSubItemPerItemByProc()
   parallel_for(number_of_sub_item_per_item.size(), PASTIS_LAMBDA(const ItemId& j){
       number_of_sub_item_per_item[j] = item_to_sub_item_matrix[j].size();
     });
-  return this->exchange(number_of_sub_item_per_item);
+
+  this->_dispatchedInfo<SubItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc =
+      this->exchange(number_of_sub_item_per_item);
 }
 
 template <int Dimension>
 template <typename SubItemOfItemT>
-std::vector<Array<const int>>
-ConnectivityDispatcher<Dimension>::
-_getRecvItemSubItemNumberingByProc(const std::vector<Array<const int>>& recv_number_of_sub_item_per_item_by_proc)
+void
+ConnectivityDispatcher<Dimension>::_buildSubItemNumbersToRecvByProc()
 {
-  std::vector<Array<const int>> item_sub_item_numbering_to_send_by_proc =
+  const std::vector<Array<const int>> sub_item_numbers_to_send_by_proc =
       [&] () {
         const auto& item_to_sub_item_matrix
             = m_connectivity.template getItemToItemMatrix<SubItemOfItemT::item_type,
                                                           SubItemOfItemT::sub_item_type>();
 
-        const ItemValue<const int, SubItemOfItemT::sub_item_type>& sub_item_number =
-            m_connectivity.template number<SubItemOfItemT::sub_item_type>();
+        const auto& sub_item_number = m_connectivity.template number<SubItemOfItemT::sub_item_type>();
 
         using ItemId = ItemIdT<SubItemOfItemT::item_type>;
         using SubItemId = ItemIdT<SubItemOfItemT::sub_item_type>;
 
-        std::vector<Array<const int>> item_sub_item_numbering_to_send_by_proc(parallel::size());
+        std::vector<Array<const int>> sub_item_numbers_to_send_by_proc(parallel::size());
         for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
           const auto& item_list_to_send_by_proc
               = this->_dispatchedInfo<SubItemOfItemT::item_type>().m_list_to_send_by_proc;
 
-          std::vector<int> sub_item_numbering_by_item_vector;
+          std::vector<int> sub_item_numbers_by_item_vector;
           for (size_t j=0; j<item_list_to_send_by_proc[i_rank].size(); ++j) {
             const ItemId& item_id = item_list_to_send_by_proc[i_rank][j];
             const auto& sub_item_list = item_to_sub_item_matrix[item_id];
             for (size_t r=0; r<sub_item_list.size(); ++r) {
               const SubItemId& sub_item_id = sub_item_list[r];
-              sub_item_numbering_by_item_vector.push_back(sub_item_number[sub_item_id]);
+              sub_item_numbers_by_item_vector.push_back(sub_item_number[sub_item_id]);
             }
           }
-          item_sub_item_numbering_to_send_by_proc[i_rank] = convert_to_array(sub_item_numbering_by_item_vector);
+          sub_item_numbers_to_send_by_proc[i_rank] = convert_to_array(sub_item_numbers_by_item_vector);
         }
-        return item_sub_item_numbering_to_send_by_proc;
+        return sub_item_numbers_to_send_by_proc;
       } ();
 
-  std::vector<Array<int>> recv_item_sub_item_numbering_by_proc(parallel::size());
+  const auto& number_of_sub_item_per_item_to_recv_by_proc =
+      this->_dispatchedInfo<SubItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc;
+
+  std::vector<Array<int>> sub_item_numbers_to_recv_by_proc(parallel::size());
+  for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
+    sub_item_numbers_to_recv_by_proc[i_rank]
+        = Array<int>(sum(number_of_sub_item_per_item_to_recv_by_proc[i_rank]));
+  }
+  parallel::exchange(sub_item_numbers_to_send_by_proc, sub_item_numbers_to_recv_by_proc);
+
+  auto& const_sub_item_numbers_to_recv_by_proc =
+      this->_dispatchedInfo<SubItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
+
+  const_sub_item_numbers_to_recv_by_proc.resize(parallel::size());
   for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-    recv_item_sub_item_numbering_by_proc[i_rank]
-        = Array<int>(sum(recv_number_of_sub_item_per_item_by_proc[i_rank]));
+    const_sub_item_numbers_to_recv_by_proc[i_rank] = sub_item_numbers_to_recv_by_proc[i_rank];
   }
-  parallel::exchange(item_sub_item_numbering_to_send_by_proc, recv_item_sub_item_numbering_by_proc);
+}
+
+template <int Dimension>
+template <typename ItemOfItemT>
+void
+ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
+{
+  constexpr ItemType item_type = ItemOfItemT::item_type;
+  constexpr ItemType sub_item_type = ItemOfItemT::sub_item_type;
+
+  const auto& item_list_to_recv_size_by_proc =
+      this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
+
+  const auto& number_of_sub_item_per_item_to_recv_by_proc =
+      this->_dispatchedInfo<ItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc;
+
+  const auto& sub_item_number_id_map =
+      this->_dispatchedInfo<sub_item_type>().m_number_to_id_map;
+
+  const auto& recv_item_of_item_numbers_by_proc =
+      this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
 
-  std::vector<Array<const int>> const_recv_item_sub_item_numbering_by_proc(parallel::size());
   for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-    const_recv_item_sub_item_numbering_by_proc[i_rank] = recv_item_sub_item_numbering_by_proc[i_rank];
+    int l=0;
+    for (size_t i=0; i<item_list_to_recv_size_by_proc[i_rank]; ++i) {
+      std::vector<unsigned int> sub_item_vector;
+      for (int k=0; k<number_of_sub_item_per_item_to_recv_by_proc[i_rank][i]; ++k) {
+        const auto& searched_sub_item_id =
+            sub_item_number_id_map.find(recv_item_of_item_numbers_by_proc[i_rank][l++]);
+        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);
+    }
   }
-  return const_recv_item_sub_item_numbering_by_proc;
 }
 
 template <int Dimension>
@@ -300,274 +411,228 @@ _buildRecvItemIdCorrespondanceByProc()
 }
 
 template <int Dimension>
+template <ItemType item_type>
 void
-ConnectivityDispatcher<Dimension>::_dispatchFaces()
+ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
 {
-  if constexpr (Dimension>1) {
-    std::vector<Array<const int>> recv_number_of_face_per_cell_by_proc =
-        _getRecvNumberOfSubItemPerItemByProc<FaceOfCell>();
-
-    std::vector<Array<const int>> recv_cell_face_numbering_by_proc
-      = this->_getRecvItemSubItemNumberingByProc<FaceOfCell>(recv_number_of_face_per_cell_by_proc);
-
-    this->_buildSubItemNumberToIdMap<FaceOfCell>(recv_cell_face_numbering_by_proc);
-
-    const std::unordered_map<int, int>& face_number_id_map =
-        this->_dispatchedInfo<ItemType::face>().m_number_to_id_map;
+  using ItemId = ItemIdT<item_type>;
 
-    this->_buildItemListToSend<ItemType::face>();
+  // Getting references
+  Array<const size_t> number_of_item_ref_list_per_proc
+      = parallel::allGather(m_connectivity.template numberOfRefItemList<item_type>());
 
-#warning this patterns repeats for each item type. Should be factorized
-    Array<unsigned int> nb_face_to_send_by_proc(parallel::size());
-    for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
-      nb_face_to_send_by_proc[i_rank] = m_dispatched_face_info.m_list_to_send_by_proc[i_rank].size();
+  const size_t number_of_item_list_sender
+      = [&] () {
+          size_t number_of_item_list_sender=0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            number_of_item_list_sender
+                += (number_of_item_ref_list_per_proc[i_rank] > 0);
+          }
+          return number_of_item_list_sender;
+        }();
+
+  if (number_of_item_list_sender > 0) {
+    if (number_of_item_list_sender > 1) {
+      perr() << __FILE__ << ':' << __LINE__ << ": "
+             << rang::fgB::red
+             <<"need to check that knowing procs know the same item_ref_lists!"
+             << rang::fg::reset << '\n';
     }
-    this->_dispatchedInfo<ItemType::face>().m_list_to_recv_size_by_proc
-        = parallel::allToAll(nb_face_to_send_by_proc);
 
-    this->_buildRecvItemIdCorrespondanceByProc<ItemType::face>();
-    this->_gatherFrom(m_connectivity.template number<ItemType::face>(), m_new_descriptor.face_number_vector);
-
-    {
-      const auto& cell_list_to_recv_size_by_proc =
-          this->_dispatchedInfo<ItemType::cell>().m_list_to_recv_size_by_proc;
-
-      for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-        int l=0;
-        for (size_t i=0; i<cell_list_to_recv_size_by_proc[i_rank]; ++i) {
-          std::vector<unsigned int> face_vector;
-          for (int k=0; k<recv_number_of_face_per_cell_by_proc[i_rank][i]; ++k) {
-            const auto& searched_face_id = face_number_id_map.find(recv_cell_face_numbering_by_proc[i_rank][l++]);
-            Assert(searched_face_id != face_number_id_map.end());
-            face_vector.push_back(searched_face_id->second);
+    if (number_of_item_list_sender < parallel::size()) {
+      const size_t sender_rank
+          = [&] () {
+              size_t i_rank = 0;
+              for (; i_rank < parallel::size(); ++i_rank) {
+                if (number_of_item_ref_list_per_proc[i_rank] > 0) {
+                  break;
+                }
+              }
+              return i_rank;
+            }();
+
+      Assert(number_of_item_list_sender < parallel::size());
+
+      // sending references tags
+      Array<RefId::TagNumberType> ref_tag_list{number_of_item_ref_list_per_proc[sender_rank]};
+      if (parallel::rank() == sender_rank){
+        for (size_t i_item_ref_list=0; i_item_ref_list<m_connectivity.template numberOfRefItemList<item_type>();
+             ++i_item_ref_list) {
+          auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
+          ref_tag_list[i_item_ref_list] = item_ref_list.refId().tagNumber();
+        }
+      }
+      parallel::broadcast(ref_tag_list, sender_rank);
+
+      // sending references name size
+      Array<size_t> ref_name_size_list{number_of_item_ref_list_per_proc[sender_rank]};
+      if (parallel::rank() == sender_rank){
+        for (size_t i_item_ref_list=0; i_item_ref_list<m_connectivity.template numberOfRefItemList<item_type>();
+             ++i_item_ref_list) {
+          auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
+          ref_name_size_list[i_item_ref_list] = item_ref_list.refId().tagName().size();
+        }
+      }
+      parallel::broadcast(ref_name_size_list, sender_rank);
+
+      // sending references name size
+      Array<RefId::TagNameType::value_type> ref_name_cat{sum(ref_name_size_list)};
+      if (parallel::rank() == sender_rank){
+        size_t i_char=0;
+        for (size_t i_item_ref_list=0; i_item_ref_list<m_connectivity.template numberOfRefItemList<item_type>();
+             ++i_item_ref_list) {
+          auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
+          for (auto c : item_ref_list.refId().tagName()) {
+            ref_name_cat[i_char++] = c;
           }
-          m_new_descriptor.cell_to_face_vector.emplace_back(face_vector);
         }
       }
-    }
+      parallel::broadcast(ref_name_cat, sender_rank);
+
+      std::vector<RefId> ref_id_list
+          = [&] () {
+              std::vector<RefId> ref_id_list;
+              ref_id_list.reserve(ref_name_size_list.size());
+              size_t begining=0;
+              for (size_t i_ref=0; i_ref < ref_name_size_list.size(); ++i_ref) {
+                const size_t size = ref_name_size_list[i_ref];
+                ref_id_list.emplace_back(ref_tag_list[i_ref],
+                                         std::string{&(ref_name_cat[begining]), size});
+                begining += size;
+              }
+              return ref_id_list;
+            } ();
 
-    {
-      std::vector<Array<bool>> cell_face_is_reversed_to_send_by_proc(parallel::size());
-      {
-        const auto& cell_list_to_send_by_proc = this->_dispatchedInfo<ItemType::cell>().m_list_to_send_by_proc;
+      using block_type = int32_t;
+      constexpr size_t block_size = sizeof(block_type);
+      const size_t nb_block = ref_id_list.size()/block_size + (ref_id_list.size()%block_size != 0);
+      for (size_t i_block=0; i_block<nb_block; ++i_block) {
+        ItemValue<block_type, item_type> item_references(m_connectivity);
+        item_references.fill(0);
 
-        const auto& cell_face_is_reversed = m_connectivity.cellFaceIsReversed();
-        for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-          std::vector<bool> face_is_reversed_by_cell_vector;
-          for (size_t j=0; j<cell_list_to_send_by_proc[i_rank].size(); ++j) {
-            const CellId& cell_id = cell_list_to_send_by_proc[i_rank][j];
-            const auto& face_is_reversed = cell_face_is_reversed.itemValues(cell_id);
-            for (size_t L=0; L<face_is_reversed.size(); ++L) {
-              face_is_reversed_by_cell_vector.push_back(face_is_reversed[L]);
+        if (m_connectivity.template numberOfRefItemList<item_type>() > 0) {
+          const size_t max_i_ref = std::min(ref_id_list.size(), block_size*(i_block+1));
+          for (size_t i_ref=block_size*i_block, i=0; i_ref<max_i_ref; ++i_ref, ++i) {
+            block_type ref_bit{1<<i};
+            auto item_ref_list = m_connectivity.template refItemList<item_type>(i_ref);
+
+            const auto& item_list = item_ref_list.list();
+            for (size_t i_item=0; i_item<item_list.size(); ++i_item) {
+              const ItemId& item_id = item_list[i_item];
+              item_references[item_id] |= ref_bit;
             }
           }
-          cell_face_is_reversed_to_send_by_proc[i_rank] = convert_to_array(face_is_reversed_by_cell_vector);
         }
-      }
 
-      std::vector<Array<bool>> recv_cell_face_is_reversed_by_proc(parallel::size());
-      for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-        recv_cell_face_is_reversed_by_proc[i_rank]
-            = Array<bool>(sum(recv_number_of_face_per_cell_by_proc[i_rank]));
-      }
+        const auto& nb_item_to_send_by_proc =
+            this->_dispatchedInfo<item_type>().m_list_to_send_size_by_proc;
 
-      parallel::exchange(cell_face_is_reversed_to_send_by_proc, recv_cell_face_is_reversed_by_proc);
+        const auto& send_item_id_by_proc =
+            this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
 
-      const auto& cell_list_to_recv_size_by_proc =
-          this->_dispatchedInfo<ItemType::cell>().m_list_to_recv_size_by_proc;
+        std::vector<Array<const block_type>> send_item_refs_by_proc(parallel::size());
 
-      for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-        int l=0;
-        for (size_t i=0; i<cell_list_to_recv_size_by_proc[i_rank]; ++i) {
-          std::vector<bool> face_is_reversed_vector;
-          for (int k=0; k<recv_number_of_face_per_cell_by_proc[i_rank][i]; ++k) {
-            face_is_reversed_vector.push_back(recv_cell_face_is_reversed_by_proc[i_rank][l++]);
+        for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+          Array<block_type> send_item_refs(nb_item_to_send_by_proc[i_rank]);
+          const Array<const ItemId> send_item_id = send_item_id_by_proc[i_rank];
+          parallel_for(send_item_id.size(), PASTIS_LAMBDA(const size_t& l) {
+              const ItemId& item_id = send_item_id[l];
+              send_item_refs[l] = item_references[item_id];
+            });
+          send_item_refs_by_proc[i_rank] = send_item_refs;
+        }
+
+        std::vector<Array<block_type>> recv_item_refs_by_proc(parallel::size());
+        for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+          recv_item_refs_by_proc[i_rank] = Array<block_type>(this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc[i_rank]);
+        }
+        parallel::exchange(send_item_refs_by_proc, recv_item_refs_by_proc);
+
+        const auto& recv_item_id_correspondance_by_proc =
+            this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
+        std::vector<block_type> item_refs(m_new_descriptor.template itemNumberVector<item_type>().size());
+        for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+          for (size_t r=0; r<recv_item_refs_by_proc[i_rank].size(); ++r) {
+            const ItemId& item_id = recv_item_id_correspondance_by_proc[i_rank][r];
+            item_refs[item_id] = recv_item_refs_by_proc[i_rank][r];
           }
-          m_new_descriptor.cell_face_is_reversed_vector.emplace_back(face_is_reversed_vector);
         }
-      }
-    }
 
-    this->_gatherFrom(this->_dispatchedInfo<ItemType::face>().m_new_owner, m_new_descriptor.face_owner_vector);
+        const size_t max_i_ref = std::min(ref_id_list.size(), block_size*(i_block+1));
+        for (size_t i_ref=block_size*i_block, i=0; i_ref<max_i_ref; ++i_ref, ++i) {
+          block_type ref_bit{1<<i};
+
+          std::vector<ItemId> item_id_vector;
 
-    std::vector<Array<const int>> recv_number_of_node_per_face_by_proc =
-        _getRecvNumberOfSubItemPerItemByProc<NodeOfFace>();
-
-    std::vector<Array<const int>> recv_face_node_numbering_by_proc
-        = this->_getRecvItemSubItemNumberingByProc<NodeOfFace>(recv_number_of_node_per_face_by_proc);
-
-    {
-      const auto& node_number_id_map = this->_dispatchedInfo<ItemType::node>().m_number_to_id_map;
-      for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-        int l=0;
-        for (size_t i=0; i<recv_number_of_node_per_face_by_proc[i_rank].size(); ++i) {
-          std::vector<unsigned int> node_vector;
-          for (int k=0; k<recv_number_of_node_per_face_by_proc[i_rank][i]; ++k) {
-            const auto& searched_node_id = node_number_id_map.find(recv_face_node_numbering_by_proc[i_rank][l++]);
-            Assert(searched_node_id != node_number_id_map.end());
-            node_vector.push_back(searched_node_id->second);
+          for (uint32_t i_item=0; i_item<item_refs.size(); ++i_item) {
+            const ItemId item_id{i_item};
+            if (item_refs[item_id] & ref_bit) {
+              item_id_vector.push_back(item_id);
+            }
           }
-          m_new_descriptor.face_to_node_vector.emplace_back(node_vector);
+
+          Array<const ItemId> item_id_array = convert_to_array(item_id_vector);
+
+          m_new_descriptor.addRefItemList(RefItemList<item_type>(ref_id_list[i_ref], item_id_array));
         }
       }
     }
+  }
+}
 
-    // Getting references
-    Array<const size_t> number_of_ref_face_list_per_proc
-        = parallel::allGather(m_connectivity.numberOfRefFaceList());
+template <int Dimension>
+void
+ConnectivityDispatcher<Dimension>::_dispatchEdges()
+{
+  if constexpr (Dimension>2) {
+    this->_buildNumberOfSubItemPerItemToRecvByProc<EdgeOfCell>();
+    this->_buildSubItemNumbersToRecvByProc<EdgeOfCell>();
+    this->_buildSubItemNumberToIdMap<EdgeOfCell>();
+    this->_buildItemToExchangeLists<ItemType::edge>();
 
-    const size_t number_of_face_list_sender
-        = [&] () {
-            size_t number_of_face_list_sender=0;
-            for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-              number_of_face_list_sender
-                  += (number_of_ref_face_list_per_proc[i_rank] > 0);
-            }
-            return number_of_face_list_sender;
-          }();
-
-    if (number_of_face_list_sender > 0) {
-      if (number_of_face_list_sender > 1) {
-        perr() << __FILE__ << ':' << __LINE__ << ": "
-               << rang::fgB::red
-               <<"need to check that knowing procs know the same ref_face_lists!"
-               << rang::fg::reset << '\n';
-      }
+    this->_gatherFrom(m_connectivity.template number<ItemType::edge>(), m_new_descriptor.edge_number_vector);
 
-      if (number_of_face_list_sender < parallel::size()) {
-        const size_t sender_rank
-            = [&] () {
-                size_t i_rank = 0;
-                for (; i_rank < parallel::size(); ++i_rank) {
-                  if (number_of_ref_face_list_per_proc[i_rank] > 0) {
-                    break;
-                  }
-                }
-                return i_rank;
-              }();
-
-        Assert(number_of_face_list_sender < parallel::size());
-
-        // sending references tags
-        Array<RefId::TagNumberType> ref_tag_list{number_of_ref_face_list_per_proc[sender_rank]};
-        if (parallel::rank() == sender_rank){
-          for (size_t i_ref_face_list=0; i_ref_face_list<m_connectivity.numberOfRefFaceList();
-               ++i_ref_face_list) {
-            auto ref_face_list = m_connectivity.refFaceList(i_ref_face_list);
-            ref_tag_list[i_ref_face_list] = ref_face_list.refId().tagNumber();
-          }
-        }
-        parallel::broadcast(ref_tag_list, sender_rank);
-
-        // sending references name size
-        Array<size_t> ref_name_size_list{number_of_ref_face_list_per_proc[sender_rank]};
-        if (parallel::rank() == sender_rank){
-          for (size_t i_ref_face_list=0; i_ref_face_list<m_connectivity.numberOfRefFaceList();
-               ++i_ref_face_list) {
-            auto ref_face_list = m_connectivity.refFaceList(i_ref_face_list);
-            ref_name_size_list[i_ref_face_list] = ref_face_list.refId().tagName().size();
-          }
-        }
-        parallel::broadcast(ref_name_size_list, sender_rank);
-
-        // sending references name size
-        Array<RefId::TagNameType::value_type> ref_name_cat{sum(ref_name_size_list)};
-        if (parallel::rank() == sender_rank){
-          size_t i_char=0;
-          for (size_t i_ref_face_list=0; i_ref_face_list<m_connectivity.numberOfRefFaceList();
-               ++i_ref_face_list) {
-            auto ref_face_list = m_connectivity.refFaceList(i_ref_face_list);
-            for (auto c : ref_face_list.refId().tagName()) {
-              ref_name_cat[i_char++] = c;
-            }
-          }
-        }
-        parallel::broadcast(ref_name_cat, sender_rank);
-
-        std::vector<RefId> ref_id_list
-            = [&] () {
-                std::vector<RefId> ref_id_list;
-                ref_id_list.reserve(ref_name_size_list.size());
-                size_t begining=0;
-                for (size_t i_ref=0; i_ref < ref_name_size_list.size(); ++i_ref) {
-                  const size_t size = ref_name_size_list[i_ref];
-                  ref_id_list.emplace_back(ref_tag_list[i_ref],
-                                           std::string{&(ref_name_cat[begining]), size});
-                  begining += size;
-                }
-                return ref_id_list;
-              } ();
-
-        using block_type = int32_t;
-        constexpr size_t block_size = sizeof(block_type);
-        const size_t nb_block = ref_id_list.size()/block_size + (ref_id_list.size()%block_size != 0);
-        for (size_t i_block=0; i_block<nb_block; ++i_block) {
-          FaceValue<block_type> face_references(m_connectivity);
-          face_references.fill(0);
-
-          if (m_connectivity.numberOfRefFaceList() > 0) {
-            const size_t max_i_ref = std::min(ref_id_list.size(), block_size*(i_block+1));
-            for (size_t i_ref=block_size*i_block, i=0; i_ref<max_i_ref; ++i_ref, ++i) {
-              block_type ref_bit{1<<i};
-              auto ref_face_list = m_connectivity.refFaceList(i_ref);
-
-              const auto& face_list = ref_face_list.faceList();
-              for (size_t i_face=0; i_face<face_list.size(); ++i_face) {
-                const FaceId& face_id = face_list[i_face];
-                face_references[face_id] |= ref_bit;
-              }
-            }
-          }
+    this->_buildItemToSubItemDescriptor<EdgeOfCell>();
 
-          const auto& send_face_id_by_proc = m_dispatched_face_info.m_list_to_send_by_proc;
-          std::vector<Array<const block_type>> send_face_refs_by_proc(parallel::size());
-          for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
-            Array<block_type> send_face_refs(nb_face_to_send_by_proc[i_rank]);
-            const Array<const FaceId> send_face_id = send_face_id_by_proc[i_rank];
-            parallel_for(send_face_id.size(), PASTIS_LAMBDA(const size_t& l) {
-                const FaceId& face_id = send_face_id[l];
-                send_face_refs[l] = face_references[face_id];
-              });
-            send_face_refs_by_proc[i_rank] = send_face_refs;
-          }
+    this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfEdge>();
+    this->_buildSubItemNumbersToRecvByProc<NodeOfEdge>();
+    this->_buildItemToSubItemDescriptor<NodeOfEdge>();
 
-          std::vector<Array<block_type>> recv_face_refs_by_proc(parallel::size());
-          for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
-            recv_face_refs_by_proc[i_rank] = Array<block_type>(m_dispatched_face_info.m_list_to_recv_size_by_proc[i_rank]);
-          }
-          parallel::exchange(send_face_refs_by_proc, recv_face_refs_by_proc);
-
-          const auto& recv_face_id_correspondance_by_proc = m_dispatched_face_info.m_recv_id_correspondance_by_proc;
-          std::vector<block_type> face_refs(m_new_descriptor.face_number_vector.size());
-          for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
-            for (size_t r=0; r<recv_face_refs_by_proc[i_rank].size(); ++r) {
-              const FaceId& face_id = recv_face_id_correspondance_by_proc[i_rank][r];
-              face_refs[face_id] = recv_face_refs_by_proc[i_rank][r];
-            }
-          }
+    this->_buildNumberOfSubItemPerItemToRecvByProc<EdgeOfFace>();
+    this->_buildSubItemNumbersToRecvByProc<EdgeOfFace>();
+    this->_buildItemToSubItemDescriptor<EdgeOfFace>();
 
-          const size_t max_i_ref = std::min(ref_id_list.size(), block_size*(i_block+1));
-          for (size_t i_ref=block_size*i_block, i=0; i_ref<max_i_ref; ++i_ref, ++i) {
-            block_type ref_bit{1<<i};
+    this->_gatherFrom(m_connectivity.faceEdgeIsReversed(), m_new_descriptor.face_edge_is_reversed_vector);
 
-            std::vector<FaceId> face_id_vector;
+    this->_gatherFrom(this->_dispatchedInfo<ItemType::edge>().m_new_owner, m_new_descriptor.edge_owner_vector);
 
-            for (uint32_t i_face=0; i_face<face_refs.size(); ++i_face) {
-              const FaceId face_id{i_face};
-              if (face_refs[face_id] & ref_bit) {
-                face_id_vector.push_back(face_id);
-              }
-            }
+    this->_buildItemReferenceList<ItemType::edge>();
+  }
+}
 
-            Array<const FaceId> face_id_array = convert_to_array(face_id_vector);
+template <int Dimension>
+void
+ConnectivityDispatcher<Dimension>::_dispatchFaces()
+{
+  if constexpr (Dimension>1) {
+    this->_buildNumberOfSubItemPerItemToRecvByProc<FaceOfCell>();
+    this->_buildSubItemNumbersToRecvByProc<FaceOfCell>();
+    this->_buildSubItemNumberToIdMap<FaceOfCell>();
+    this->_buildItemToExchangeLists<ItemType::face>();
 
-            m_new_descriptor.addRefFaceList(RefFaceList(ref_id_list[i_ref], face_id_array));
-          }
+    this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfFace>();
+    this->_buildSubItemNumbersToRecvByProc<NodeOfFace>();
+    this->_buildItemToSubItemDescriptor<NodeOfFace>();
 
-          pout() << __FILE__ << ':' << __LINE__ << ": remains to build lists\n";
-        }
-      }
-    }
+    this->_gatherFrom(m_connectivity.template number<ItemType::face>(), m_new_descriptor.face_number_vector);
+
+    this->_buildItemToSubItemDescriptor<FaceOfCell>();
+
+    this->_gatherFrom(m_connectivity.cellFaceIsReversed(), m_new_descriptor.cell_face_is_reversed_vector);
+
+    this->_gatherFrom(this->_dispatchedInfo<ItemType::face>().m_new_owner, m_new_descriptor.face_owner_vector);
+
+    this->_buildItemReferenceList<ItemType::face>();
   }
 }
 
@@ -578,35 +643,20 @@ ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType
 {
   this->_buildNewOwner<ItemType::cell>();
   this->_buildNewOwner<ItemType::face>();
-  // this->_buildNewOwner<ItemType::edge>();
+  this->_buildNewOwner<ItemType::edge>();
   this->_buildNewOwner<ItemType::node>();
 
-  this->_buildItemListToSend<ItemType::cell>();
-  this->_dispatchedInfo<ItemType::cell>().m_list_to_recv_size_by_proc
-      = parallel::allToAll(this->_buildNbCellToSend());
-
-  this->_buildCellNumberIdMap();
+  this->_buildItemToExchangeLists<ItemType::cell>();
 
-  const std::vector<Array<const int>> recv_number_of_node_per_cell_by_proc
-      = this->_getRecvNumberOfSubItemPerItemByProc<NodeOfCell>();
+  this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfCell>();
 
-  const std::vector<Array<const int>> recv_cell_node_numbering_by_proc
-      = this->_getRecvItemSubItemNumberingByProc<NodeOfCell>(recv_number_of_node_per_cell_by_proc);
+  this->_buildSubItemNumbersToRecvByProc<NodeOfCell>();
 
-  this->_buildRecvItemIdCorrespondanceByProc<ItemType::cell>();
   this->_gatherFrom(m_connectivity.template number<ItemType::cell>(), m_new_descriptor.cell_number_vector);
 
-  this->_buildSubItemNumberToIdMap<NodeOfCell>(recv_cell_node_numbering_by_proc);
-  this->_buildItemListToSend<ItemType::node>();
+  this->_buildSubItemNumberToIdMap<NodeOfCell>();
 
-  Array<unsigned int> nb_node_to_send_by_proc(parallel::size());
-  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
-    nb_node_to_send_by_proc[i_rank] = m_dispatched_node_info.m_list_to_send_by_proc[i_rank].size();
-  }
-  this->_dispatchedInfo<ItemType::node>().m_list_to_recv_size_by_proc
-      = parallel::allToAll(nb_node_to_send_by_proc);
-
-  this->_buildRecvItemIdCorrespondanceByProc<ItemType::node>();
+  this->_buildItemToExchangeLists<ItemType::node>();
 
   // Fill new descriptor
   this->_gatherFrom(m_connectivity.cellType(), m_new_descriptor.cell_type_vector);
@@ -615,26 +665,16 @@ ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType
   this->_gatherFrom(m_connectivity.template number<ItemType::node>(), m_new_descriptor.node_number_vector);
   this->_gatherFrom(this->_dispatchedInfo<ItemType::node>().m_new_owner, m_new_descriptor.node_owner_vector);
 
-  { // build cells connectivity
-    const auto& cell_list_to_recv_size_by_proc =
-        this->_dispatchedInfo<ItemType::cell>().m_list_to_recv_size_by_proc;
-
-    const auto& node_number_id_map  = this->_dispatchedInfo<ItemType::node>().m_number_to_id_map;
-    for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
-      int l=0;
-      for (size_t i=0; i<cell_list_to_recv_size_by_proc[i_rank]; ++i) {
-        std::vector<unsigned int> node_vector;
-        for (int k=0; k<recv_number_of_node_per_cell_by_proc[i_rank][i]; ++k) {
-          const auto& searched_node_id = node_number_id_map.find(recv_cell_node_numbering_by_proc[i_rank][l++]);
-          Assert(searched_node_id != node_number_id_map.end());
-          node_vector.push_back(searched_node_id->second);
-        }
-        m_new_descriptor.cell_by_node_vector.emplace_back(node_vector);
-      }
-    }
-  }
+  this->_buildItemToSubItemDescriptor<NodeOfCell>();
+
+  this->_buildItemReferenceList<ItemType::cell>();
+
   this->_dispatchFaces();
 
+  this->_dispatchEdges();
+
+  this->_buildItemReferenceList<ItemType::node>();
+
   m_dispatched_connectivity = ConnectivityType::build(m_new_descriptor);
 }
 
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index 947821f08c4f7c070d3919651567bb0c60210047..e0ce2e88e2690093b6c8d4c16f7f2abc70346674 100644
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -6,6 +6,7 @@
 #include <ItemValueUtils.hpp>
 
 #include <unordered_map>
+#include <ConnectivityDescriptor.hpp>
 
 template <int Dimension>
 class ConnectivityDispatcher
@@ -23,8 +24,8 @@ class ConnectivityDispatcher
   {
     using ItemId = ItemIdT<item_type>;
     ItemValue<const int, item_type> m_new_owner;
+    Array<const unsigned int> m_list_to_send_size_by_proc;
     std::vector<Array<const ItemId>> m_list_to_send_by_proc;
-#warning is m_list_to_recv_size_by_proc really necessary?
     Array<const unsigned int> m_list_to_recv_size_by_proc;
     std::unordered_map<int, int> m_number_to_id_map;
     std::vector<Array<const ItemId>> m_recv_id_correspondance_by_proc;
@@ -65,36 +66,138 @@ class ConnectivityDispatcher
     }
   }
 
+  template <typename ItemToItem>
+  struct DispatchedItemOfItemInfo
+  {
+    std::vector<Array<const int>> m_number_of_sub_item_per_item_to_recv_by_proc;
+    std::vector<Array<const int>> m_sub_item_numbers_to_recv_by_proc;
+  };
+
+  DispatchedItemOfItemInfo<NodeOfCell> m_dispatched_node_of_cell_info;
+  DispatchedItemOfItemInfo<EdgeOfCell> m_dispatched_edge_of_cell_info;
+  DispatchedItemOfItemInfo<FaceOfCell> m_dispatched_face_of_cell_info;
+
+  DispatchedItemOfItemInfo<NodeOfEdge> m_dispatched_node_of_edge_info;
+  DispatchedItemOfItemInfo<FaceOfEdge> m_dispatched_face_of_edge_info;
+  DispatchedItemOfItemInfo<CellOfEdge> m_dispatched_cell_of_edge_info;
+
+  DispatchedItemOfItemInfo<NodeOfFace> m_dispatched_node_of_face_info;
+  DispatchedItemOfItemInfo<EdgeOfFace> m_dispatched_edge_of_face_info;
+  DispatchedItemOfItemInfo<CellOfFace> m_dispatched_cell_of_face_info;
+
+  DispatchedItemOfItemInfo<EdgeOfNode> m_dispatched_edge_of_node_info;
+  DispatchedItemOfItemInfo<FaceOfNode> m_dispatched_face_of_node_info;
+  DispatchedItemOfItemInfo<CellOfNode> m_dispatched_cell_of_node_info;
+
+  template <typename ItemOfItem>
+  PASTIS_INLINE
+  DispatchedItemOfItemInfo<ItemOfItem>& _dispatchedInfo()
+  {
+    if constexpr (std::is_same_v<NodeOfCell, ItemOfItem>) {
+      return m_dispatched_node_of_cell_info;
+    } else if constexpr (std::is_same_v<EdgeOfCell, ItemOfItem>) {
+      return m_dispatched_edge_of_cell_info;
+    } else if constexpr (std::is_same_v<FaceOfCell, ItemOfItem>) {
+      return m_dispatched_face_of_cell_info;
+    } else if constexpr (std::is_same_v<NodeOfEdge, ItemOfItem>) {
+      return m_dispatched_node_of_edge_info;
+    } else if constexpr (std::is_same_v<FaceOfEdge, ItemOfItem>) {
+      return m_dispatched_face_of_edge_info;
+    } else if constexpr (std::is_same_v<CellOfEdge, ItemOfItem>) {
+      return m_dispatched_cell_of_edge_info;
+    } else if constexpr (std::is_same_v<NodeOfFace, ItemOfItem>) {
+      return m_dispatched_node_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfFace, ItemOfItem>) {
+      return m_dispatched_edge_of_face_info;
+    } else if constexpr (std::is_same_v<CellOfFace, ItemOfItem>) {
+      return m_dispatched_cell_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfNode, ItemOfItem>) {
+      return m_dispatched_edge_of_node_info;
+    } else if constexpr (std::is_same_v<FaceOfNode, ItemOfItem>) {
+      return m_dispatched_face_of_node_info;
+    } else if constexpr (std::is_same_v<CellOfNode, ItemOfItem>) {
+      return m_dispatched_cell_of_node_info;
+    } else {
+      static_assert(is_false_v<ItemOfItem>, "Unexpected ItemOfItem type");
+    }
+  }
+
+  template <typename ItemOfItem>
+  PASTIS_INLINE
+  const DispatchedItemOfItemInfo<ItemOfItem>& _dispatchedInfo() const
+  {
+    if constexpr (std::is_same_v<NodeOfCell, ItemOfItem>) {
+      return m_dispatched_node_of_cell_info;
+    } else if constexpr (std::is_same_v<EdgeOfCell, ItemOfItem>) {
+      return m_dispatched_edge_of_cell_info;
+    } else if constexpr (std::is_same_v<FaceOfCell, ItemOfItem>) {
+      return m_dispatched_face_of_cell_info;
+    } else if constexpr (std::is_same_v<NodeOfEdge, ItemOfItem>) {
+      return m_dispatched_node_of_edge_info;
+    } else if constexpr (std::is_same_v<FaceOfEdge, ItemOfItem>) {
+      return m_dispatched_face_of_edge_info;
+    } else if constexpr (std::is_same_v<CellOfEdge, ItemOfItem>) {
+      return m_dispatched_cell_of_edge_info;
+    } else if constexpr (std::is_same_v<NodeOfFace, ItemOfItem>) {
+      return m_dispatched_node_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfFace, ItemOfItem>) {
+      return m_dispatched_edge_of_face_info;
+    } else if constexpr (std::is_same_v<CellOfFace, ItemOfItem>) {
+      return m_dispatched_cell_of_face_info;
+    } else if constexpr (std::is_same_v<EdgeOfNode, ItemOfItem>) {
+      return m_dispatched_edge_of_node_info;
+    } else if constexpr (std::is_same_v<FaceOfNode, ItemOfItem>) {
+      return m_dispatched_face_of_node_info;
+    } else if constexpr (std::is_same_v<CellOfNode, ItemOfItem>) {
+      return m_dispatched_cell_of_node_info;
+    } else {
+      static_assert(is_false_v<ItemOfItem>, "Unexpected ItemOfItem type");
+    }
+  }
+
   template <ItemType item_type>
   void _buildNewOwner();
 
   template <ItemType item_type>
   void _buildItemListToSend();
 
-  Array<const unsigned int> _buildNbCellToSend();
-
   void _buildCellNumberIdMap();
 
   template <typename ItemOfItemT>
-  void _buildSubItemNumberToIdMap(const std::vector<Array<const int>>& recv_cell_node_number_by_proc);
+  void _buildSubItemNumberToIdMap();
+
+  template <ItemType item_type>
+  void _buildItemToExchangeLists();
+
+  template <ItemType item_type>
+  void _buildNumberOfItemToExchange();
+
+  template <typename ItemOfItemT>
+  void _buildItemToSubItemDescriptor();
 
+  void _dispatchEdges();
   void _dispatchFaces();
 
   template<typename DataType, ItemType item_type, typename ConnectivityPtr>
   void _gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
                    std::vector<std::remove_const_t<DataType>>& gathered_vector);
 
+  template<typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+  void _gatherFrom(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
+                   std::vector<Array<std::remove_const_t<DataType>>>& gathered_vector);
+
   template <typename SubItemOfItemT>
-  std::vector<Array<const int>>
-  _getRecvNumberOfSubItemPerItemByProc();
+  void _buildNumberOfSubItemPerItemToRecvByProc();
 
   template <typename SubItemOfItemT>
-  std::vector<Array<const int>>
-  _getRecvItemSubItemNumberingByProc(const std::vector<Array<const int>>& recv_number_of_sub_item_per_item_by_proc);
+  void _buildSubItemNumbersToRecvByProc();
 
   template <ItemType item_type>
   void _buildRecvItemIdCorrespondanceByProc();
 
+  template <ItemType item_type>
+  void _buildItemReferenceList();
+
  public:
   std::shared_ptr<const ConnectivityType>
   dispatchedConnectivity() const
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index fde501ac1fea1516e6bcb1577be60ea0c1d1efa3..1881f99a4eb55e8f6de284ee2b91a05a412aba04 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -12,7 +12,7 @@
 #include <Mesh.hpp>
 #include <MeshData.hpp>
 
-#include <RefFaceList.hpp>
+#include <RefItemList.hpp>
 #include <Messenger.hpp>
 
 #include <ArrayUtils.hpp>
@@ -208,12 +208,6 @@ class ConnectivityFace<2>
              (m_node_number_vector[m_node1_id]<m_node_number_vector[f.m_node1_id])));
   }
 
-  PASTIS_INLINE
-  ConnectivityFace& operator=(const ConnectivityFace&) = default;
-
-  PASTIS_INLINE
-  ConnectivityFace& operator=(ConnectivityFace&&) = default;
-
   PASTIS_INLINE
   ConnectivityFace(const std::vector<unsigned int>& node_id_list,
                    const std::vector<int>& node_number_vector)
@@ -235,9 +229,6 @@ class ConnectivityFace<2>
   PASTIS_INLINE
   ConnectivityFace(const ConnectivityFace&) = default;
 
-  PASTIS_INLINE
-  ConnectivityFace(ConnectivityFace&&) = default;
-
   PASTIS_INLINE
   ~ConnectivityFace() = default;
 };
@@ -333,19 +324,9 @@ class ConnectivityFace<3>
     return m_node_id_list.size() < f.m_node_id_list.size();
   }
 
-  PASTIS_INLINE
-  ConnectivityFace& operator=(const ConnectivityFace&) = default;
-
-  PASTIS_INLINE
-  ConnectivityFace& operator=(ConnectivityFace&&) = default;
-
   PASTIS_INLINE
   ConnectivityFace(const ConnectivityFace&) = default;
 
-  PASTIS_INLINE
-  ConnectivityFace(ConnectivityFace&&) = default;
-
-
   PASTIS_INLINE
   ConnectivityFace() = delete;
 
@@ -797,10 +778,10 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
   using Face = ConnectivityFace<Dimension>;
 
   const auto& node_number_vector = descriptor.node_number_vector;
-  Array<unsigned short> cell_nb_faces(descriptor.cell_by_node_vector.size());
+  Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
   std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
-  for (CellId j=0; j<descriptor.cell_by_node_vector.size(); ++j) {
-    const auto& cell_nodes = descriptor.cell_by_node_vector[j];
+  for (CellId j=0; j<descriptor.cell_to_node_vector.size(); ++j) {
+    const auto& cell_nodes = descriptor.cell_to_node_vector[j];
 
     if constexpr (Dimension==2) {
       switch (descriptor.cell_type_vector[j]) {
@@ -929,7 +910,7 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
   }
 
   {
-    descriptor.cell_to_face_vector.resize(descriptor.cell_by_node_vector.size());
+    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]);
     }
@@ -937,7 +918,7 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
     for (const auto& face_cells_vector : face_cells_map) {
       const auto& cells_vector = face_cells_vector.second;
       for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
-        const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
+        const auto& [cell_number, cell_local_face,  reversed] = cells_vector[lj];
         descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
       }
       ++l;
@@ -945,9 +926,9 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
   }
 
   {
-    descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_by_node_vector.size());
+    descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
     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]);
+      descriptor.cell_face_is_reversed_vector[j] = Array<bool>(cell_nb_faces[j]);
     }
     for (const auto& face_cells_vector : face_cells_map) {
       const auto& cells_vector = face_cells_vector.second;
@@ -978,6 +959,134 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
 }
 
 
+template <size_t Dimension>
+void
+GmshReader::__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
+{
+  static_assert(Dimension==3, "Invalid dimension to compute face-edge connectivities");
+  using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
+  using Edge = ConnectivityFace<2>;
+
+  const auto& node_number_vector = descriptor.node_number_vector;
+  Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
+  std::map<Edge, std::vector<FaceEdgeInfo>> edge_faces_map;
+  for (FaceId l=0; l<descriptor.face_to_node_vector.size(); ++l) {
+    const auto& face_nodes = descriptor.face_to_node_vector[l];
+
+    face_nb_edges[l] = face_nodes.size();
+    for (size_t r=0; r<face_nodes.size()-1; ++r) {
+      Edge e({face_nodes[r], face_nodes[r+1]}, node_number_vector);
+      edge_faces_map[e].emplace_back(std::make_tuple(l, r, e.reversed()));
+    }
+    {
+      Edge e({face_nodes[face_nodes.size()-1], face_nodes[0]}, node_number_vector);
+      edge_faces_map[e].emplace_back(std::make_tuple(l, face_nodes.size()-1, e.reversed()));
+    }
+  }
+
+  std::unordered_map<Edge,EdgeId,typename Edge::Hash> edge_id_map;
+  {
+    descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
+    for (FaceId l=0; l<descriptor.face_to_node_vector.size(); ++l) {
+      descriptor.face_to_edge_vector[l].resize(face_nb_edges[l]);
+    }
+    EdgeId e=0;
+    for (const auto& edge_faces_vector : edge_faces_map) {
+      const auto& faces_vector = edge_faces_vector.second;
+      for (unsigned short l=0; l<faces_vector.size(); ++l) {
+        const auto& [face_number, face_local_edge,  reversed] = faces_vector[l];
+        descriptor.face_to_edge_vector[face_number][face_local_edge] = e;
+      }
+      edge_id_map[edge_faces_vector.first] = e;
+      ++e;
+    }
+  }
+
+  {
+    descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
+    for (FaceId j=0; j<descriptor.face_edge_is_reversed_vector.size(); ++j) {
+      descriptor.face_edge_is_reversed_vector[j] = Array<bool>(face_nb_edges[j]);
+    }
+    for (const auto& edge_faces_vector : edge_faces_map) {
+      const auto& faces_vector = edge_faces_vector.second;
+      for (unsigned short lj=0; lj<faces_vector.size(); ++lj) {
+        const auto& [face_number, face_local_edge, reversed] = faces_vector[lj];
+        descriptor.face_edge_is_reversed_vector[face_number][face_local_edge] = reversed;
+      }
+    }
+  }
+
+  {
+    descriptor.edge_to_node_vector.resize(edge_faces_map.size());
+    int e=0;
+    for (const auto& edge_info : edge_faces_map) {
+      const Edge& edge = edge_info.first;
+      descriptor.edge_to_node_vector[e] = edge.nodeIdList();
+      ++e;
+    }
+  }
+
+  {
+    // Edge numbers may change if numbers are provided in the file
+    descriptor.edge_number_vector.resize(edge_faces_map.size());
+    for (size_t e=0; e<edge_faces_map.size(); ++e) {
+      descriptor.edge_number_vector[e] = e;
+    }
+  }
+
+  {
+    descriptor.cell_to_node_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];
+
+      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];
+          Edge edge{{cell_nodes[e[0]],cell_nodes[e[1]]}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            std::cerr << "could not find this edge!\n";
+            std::terminate();
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+        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];
+          Edge edge{{cell_nodes[e[0]],cell_nodes[e[1]]}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            std::cerr << "could not find this edge!\n";
+            std::terminate();
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+        break;
+      }
+      default: {
+        perr() << name(descriptor.cell_type_vector[j])
+               << ": unexpected cell type in dimension 3!\n";
+        std::terminate();
+      }
+      }
+    }
+  }
+}
+
+
 void
 GmshReader::__proceedData()
 {
@@ -1082,7 +1191,6 @@ GmshReader::__proceedData()
                        ErrorHandler::normal);
   }
 
-#warning should use an unordered_map
   // A value of -1 means that the vertex is unknown
   __verticesCorrepondance.resize(maxNumber+1,-1);
 
@@ -1263,13 +1371,13 @@ GmshReader::__proceedData()
     descriptor.node_number_vector = __verticesNumbers;
     descriptor.cell_type_vector.resize(nb_cells);
     descriptor.cell_number_vector.resize(nb_cells);
-    descriptor.cell_by_node_vector.resize(nb_cells);
+    descriptor.cell_to_node_vector.resize(nb_cells);
 
     const size_t nb_tetrahedra = __tetrahedra.size();
     for (size_t j=0; j<nb_tetrahedra; ++j) {
-      descriptor.cell_by_node_vector[j].resize(4);
+      descriptor.cell_to_node_vector[j].resize(4);
       for (int r=0; r<4; ++r) {
-        descriptor.cell_by_node_vector[j][r] = __tetrahedra[j][r];
+        descriptor.cell_to_node_vector[j][r] = __tetrahedra[j][r];
       }
       descriptor.cell_type_vector[j] = CellType::Tetrahedron;
       descriptor.cell_number_vector[j] = __tetrahedra_number[j];
@@ -1277,134 +1385,248 @@ GmshReader::__proceedData()
     const size_t nb_hexahedra = __hexahedra.size();
     for (size_t j=0; j<nb_hexahedra; ++j) {
       const size_t jh = nb_tetrahedra+j;
-      descriptor.cell_by_node_vector[jh].resize(8);
+      descriptor.cell_to_node_vector[jh].resize(8);
       for (int r=0; r<8; ++r) {
-        descriptor.cell_by_node_vector[jh][r] = __hexahedra[j][r];
+        descriptor.cell_to_node_vector[jh][r] = __hexahedra[j][r];
       }
       descriptor.cell_type_vector[jh] = CellType::Hexahedron;
       descriptor.cell_number_vector[jh] = __hexahedra_number[j];
     }
 
-    __computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
+    std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
+    for (unsigned int r=0; r<__tetrahedra_ref.size(); ++r) {
+      const unsigned int elem_number = __tetrahedra_ref[r];
+      const unsigned int& ref = __tetrahedra_ref[r];
+      ref_cells_map[ref].push_back(elem_number);
+    }
 
-    descriptor.cell_owner_vector.resize(nb_cells);
-    std::fill(descriptor.cell_owner_vector.begin(),
-              descriptor.cell_owner_vector.end(),
-              parallel::rank());
+    for (unsigned int j=0; j<__hexahedra_ref.size(); ++j) {
+      const size_t elem_number = nb_tetrahedra+j;
+      const unsigned int& ref = __hexahedra_ref[j];
+      ref_cells_map[ref].push_back(elem_number);
+    }
 
-    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.node_owner_vector.resize(descriptor.node_number_vector.size());
-    std::fill(descriptor.node_owner_vector.begin(),
-              descriptor.node_owner_vector.end(),
-              parallel::rank());
+    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_physical_ref_map.at(ref_cell_list.first);
+      descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
+    }
 
-    const auto& node_number_vector = descriptor.node_number_vector;
+    this->__computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
 
-    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;
-          } ();
+    const auto& node_number_vector = descriptor.node_number_vector;
 
-    std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
-    for (unsigned int f=0; f<__triangles.size(); ++f) {
-      const unsigned int face_id
-          = [&]{
-              auto i = face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]},
-                                                node_number_vector));
-              if (i == face_to_id_map.end()) {
-                std::cerr << "face not found!\n";
-                std::terminate();
+    {
+      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 i->second;
-            }();
-
-      const unsigned int& ref = __triangles_ref[f];
-      ref_faces_map[ref].push_back(face_id);
-
-      if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
-        if (auto i_face = face_number_id_map.find(__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]);
+              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<__triangles.size(); ++f) {
+        const unsigned int face_id
+            = [&]{
+                auto i = face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]},
+                                                  node_number_vector));
+                if (i == face_to_id_map.end()) {
+                  std::cerr << "face not found!\n";
+                  std::terminate();
+                }
+                return i->second;
+              }();
+
+        const unsigned int& ref = __triangles_ref[f];
+        ref_faces_map[ref].push_back(face_id);
+
+        if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
+          if (auto i_face = face_number_id_map.find(__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] = __quadrangles_number[f];
+            face_number_id_map[descriptor.face_number_vector[face_id]] = 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]);
+      for (unsigned int f=0; f<__quadrangles.size(); ++f) {
+        const unsigned int face_id
+            = [&]{
+                auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
+                                                   __quadrangles[f][2], __quadrangles[f][3]}, node_number_vector));
+                if (i == face_to_id_map.end()) {
+                  std::cerr << "face not found!\n";
+                  std::terminate();
+                }
+                return i->second;
+              }();
+
+        const unsigned int& ref = __quadrangles_ref[f];
+        ref_faces_map[ref].push_back(face_id);
+
+        if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
+          if (auto i_face = face_number_id_map.find(__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] = __quadrangles_number[f];
+            face_number_id_map[descriptor.face_number_vector[face_id]] = 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] = __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_physical_ref_map.at(ref_face_list.first);
+        descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
       }
     }
-
-    for (unsigned int f=0; f<__quadrangles.size(); ++f) {
-      const unsigned int face_id
-          = [&]{
-              auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
-                                                 __quadrangles[f][2], __quadrangles[f][3]}, node_number_vector));
-              if (i == face_to_id_map.end()) {
-                std::cerr << "face not found!\n";
-                std::terminate();
+    this->__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 i->second;
-            }();
-
-      const unsigned int& ref = __quadrangles_ref[f];
-      ref_faces_map[ref].push_back(face_id);
-
-      if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
-        if (auto i_face = face_number_id_map.find(__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]);
+              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<__edges.size(); ++e) {
+        const unsigned int edge_id
+            = [&]{
+                auto i = edge_to_id_map.find(Edge({__edges[e][0], __edges[e][1]}, node_number_vector));
+                if (i == edge_to_id_map.end()) {
+                  std::cerr << "edge " << __edges[e][0] << " not found!\n";
+                  std::terminate();
+                }
+                return i->second;
+              }();
+        const unsigned int& ref = __edges_ref[e];
+        ref_edges_map[ref].push_back(edge_id);
+
+        if (descriptor.edge_number_vector[edge_id] != __edges_number[e]) {
+          if (auto i_edge = edge_number_id_map.find(__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] = __edges_number[e];
+            edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_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] = __quadrangles_number[f];
-          face_number_id_map[descriptor.face_number_vector[face_id]] = face_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_physical_ref_map.at(ref_edge_list.first);
+        descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list});
       }
     }
 
-    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];
+    std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+    for (unsigned int r=0; r<__points.size(); ++r) {
+      const unsigned int point_number = __points[r];
+      const unsigned int& ref = __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_physical_ref_map.at(ref_face_list.first);
-      descriptor.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
+      const PhysicalRefId& physical_ref_id = 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());
+
     std::shared_ptr p_connectivity = Connectivity3D::build(descriptor);
     Connectivity3D& connectivity = *p_connectivity;
 
@@ -1427,13 +1649,13 @@ GmshReader::__proceedData()
     descriptor.node_number_vector = __verticesNumbers;
     descriptor.cell_type_vector.resize(nb_cells);
     descriptor.cell_number_vector.resize(nb_cells);
-    descriptor.cell_by_node_vector.resize(nb_cells);
+    descriptor.cell_to_node_vector.resize(nb_cells);
 
     const size_t nb_triangles = __triangles.size();
     for (size_t j=0; j<nb_triangles; ++j) {
-      descriptor.cell_by_node_vector[j].resize(3);
+      descriptor.cell_to_node_vector[j].resize(3);
       for (int r=0; r<3; ++r) {
-        descriptor.cell_by_node_vector[j][r] = __triangles[j][r];
+        descriptor.cell_to_node_vector[j][r] = __triangles[j][r];
       }
       descriptor.cell_type_vector[j] = CellType::Triangle;
       descriptor.cell_number_vector[j] = __triangles_number[j];
@@ -1442,30 +1664,37 @@ GmshReader::__proceedData()
     const size_t nb_quadrangles = __quadrangles.size();
     for (size_t j=0; j<nb_quadrangles; ++j) {
       const size_t jq = j+nb_triangles;
-      descriptor.cell_by_node_vector[jq].resize(4);
+      descriptor.cell_to_node_vector[jq].resize(4);
       for (int r=0; r<4; ++r) {
-        descriptor.cell_by_node_vector[jq][r] = __quadrangles[j][r];
+        descriptor.cell_to_node_vector[jq][r] = __quadrangles[j][r];
       }
       descriptor.cell_type_vector[jq] = CellType::Quadrangle;
       descriptor.cell_number_vector[jq] = __quadrangles_number[j];
     }
 
-    this->__computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
+    std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
+    for (unsigned int r=0; r<__triangles_ref.size(); ++r) {
+      const unsigned int elem_number = __triangles_ref[r];
+      const unsigned int& ref = __triangles_ref[r];
+      ref_cells_map[ref].push_back(elem_number);
+    }
 
-    descriptor.cell_owner_vector.resize(nb_cells);
-    std::fill(descriptor.cell_owner_vector.begin(),
-              descriptor.cell_owner_vector.end(),
-              parallel::rank());
+    for (unsigned int j=0; j<__quadrangles_ref.size(); ++j) {
+      const size_t elem_number = nb_triangles+j;
+      const unsigned int& ref = __quadrangles_ref[j];
+      ref_cells_map[ref].push_back(elem_number);
+    }
 
-    descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
-    std::fill(descriptor.face_owner_vector.begin(),
-              descriptor.face_owner_vector.end(),
-              parallel::rank());
+    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_physical_ref_map.at(ref_cell_list.first);
+      descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
+    }
 
-    descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
-    std::fill(descriptor.node_owner_vector.begin(),
-              descriptor.node_owner_vector.end(),
-              parallel::rank());
+    this->__computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
 
     using Face = ConnectivityFace<2>;
     const auto& node_number_vector = descriptor.node_number_vector;
@@ -1529,12 +1758,9 @@ GmshReader::__proceedData()
         face_list[j]=ref_face_list.second[j];
       }
       const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
-      descriptor.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
+      descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
     }
 
-    std::shared_ptr p_connectivity = Connectivity2D::build(descriptor);
-    Connectivity2D& connectivity = *p_connectivity;
-
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
     for (unsigned int r=0; r<__points.size(); ++r) {
       const unsigned int point_number = __points[r];
@@ -1548,9 +1774,27 @@ GmshReader::__proceedData()
         point_list[j]=ref_point_list.second[j];
       }
       const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
-      connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
+      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.node_owner_vector.resize(descriptor.node_number_vector.size());
+    std::fill(descriptor.node_owner_vector.begin(),
+              descriptor.node_owner_vector.end(),
+              parallel::rank());
+
+    std::shared_ptr p_connectivity = Connectivity2D::build(descriptor);
+    Connectivity2D& connectivity = *p_connectivity;
+
     using MeshType = Mesh<Connectivity2D>;
     using Rd = TinyVector<2, double>;
 
@@ -1569,30 +1813,17 @@ GmshReader::__proceedData()
     descriptor.node_number_vector = __verticesNumbers;
     descriptor.cell_type_vector.resize(nb_cells);
     descriptor.cell_number_vector.resize(nb_cells);
-    descriptor.cell_by_node_vector.resize(nb_cells);
+    descriptor.cell_to_node_vector.resize(nb_cells);
 
     for (size_t j=0; j<nb_cells; ++j) {
-      descriptor.cell_by_node_vector[j].resize(2);
+      descriptor.cell_to_node_vector[j].resize(2);
       for (int r=0; r<2; ++r) {
-        descriptor.cell_by_node_vector[j][r] = __edges[j][r];
+        descriptor.cell_to_node_vector[j][r] = __edges[j][r];
       }
       descriptor.cell_type_vector[j] = CellType::Line;
       descriptor.cell_number_vector[j] =  __edges_number[j];
     }
 
-    descriptor.cell_owner_vector.resize(nb_cells);
-    std::fill(descriptor.cell_owner_vector.begin(),
-              descriptor.cell_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());
-
-
-    std::shared_ptr p_connectivity = Connectivity1D::build(descriptor);
-    Connectivity1D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
     for (unsigned int r=0; r<__points.size(); ++r) {
@@ -1607,9 +1838,22 @@ GmshReader::__proceedData()
         point_list[j]=ref_point_list.second[j];
       }
       const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
-      connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
+      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.node_owner_vector.resize(descriptor.node_number_vector.size());
+    std::fill(descriptor.node_owner_vector.begin(),
+              descriptor.node_owner_vector.end(),
+              parallel::rank());
+
+    std::shared_ptr p_connectivity = Connectivity1D::build(descriptor);
+    Connectivity1D& connectivity = *p_connectivity;
+
     using MeshType = Mesh<Connectivity1D>;
     using Rd = TinyVector<1, double>;
 
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 1ddb6c4db32dbb3deac2fe6f87364532f4da84aa..cd482e637c90a68c1bd208a51285db06d5f93741 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -262,6 +262,9 @@ private:
   template <size_t Dimension>
   void __computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
 
+  template <size_t Dimension>
+  void __computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
+
   template <int Dimension>
   void _dispatch();
 
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index d10aa5ac11e33c95aa54012828944c45a34024c1..bf62537d33b805aee07fddd9e886aa2ed48a922c 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -83,7 +83,7 @@ class ItemToItemMatrix
   PASTIS_INLINE
   const auto& operator[](const IndexType& source_id) const
   {
-    static_assert(std::is_same<IndexType, SourceItemId>(),
+    static_assert(std::is_same_v<IndexType, SourceItemId>,
                   "ItemToItemMatrix must be indexed using correct ItemId");
     using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
     return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
diff --git a/src/mesh/ItemType.hpp b/src/mesh/ItemType.hpp
index 1b5f3cb452e898dc957ff7f96fb07c4955312f16..e9de6153098342eeaa0306e3754fd8b03e45644c 100644
--- a/src/mesh/ItemType.hpp
+++ b/src/mesh/ItemType.hpp
@@ -118,4 +118,8 @@ struct ItemTypeId<3>
   }
 };
 
+template <ItemType item_type>
+PASTIS_INLINE
+constexpr bool is_false_item_type_v = false;
+
 #endif // ITEM_TYPE_HPP
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index ef10317521a3d26460dd0166de56679e84f11cd9..827daa32d0e912dd43c13f4113f5ceb2566f034d 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -95,13 +95,10 @@ class ItemValue
   }
 
   template <typename IndexType>
-  DataType& operator[](const IndexType& i) const noexcept(NO_ASSERT)
+  DataType& operator[](const IndexType&) const noexcept(NO_ASSERT)
   {
     static_assert(std::is_same_v<IndexType,ItemId>,
                   "ItemValue must be indexed by ItemId");
-    static_assert(not std::is_const_v<DataType>,
-                  "Cannot modify ItemValue of const");
-    return m_values[i];
   }
 
   PASTIS_INLINE
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index e111d14531c050961ee581841cc3ffedf936e966..b103e22946019e8ff2c50a056e8c86a75fbf4c05 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -11,7 +11,6 @@
 struct IMesh
 {
   virtual size_t dimension() const = 0;
-  // virtual CSRGraph cellToCellGraph() const = 0;
   ~IMesh() = default;
 };
 
@@ -66,7 +65,6 @@ public:
     return m_xr;
   }
 
-  [[deprecated("should rework this class: quite ugly")]]
   PASTIS_INLINE
   NodeValue<Rd> mutableXr() const
   {
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index dd8f33deac3c96a4fc97e7900cf4dbabd9e31ece..3d5cc9b753ae4835e78f70b729b831a1cffc06ba 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -178,13 +178,11 @@ class MeshData
               const FaceId& l = cell_faces[L];
               const auto& face_nodes = face_to_node_matrix[l];
 
-#warning should this lambda be replaced by a precomputed correspondance?
               auto local_node_number_in_cell
                   = [&](const NodeId& node_number) {
                       for (size_t i_node=0; i_node<cell_nodes.size(); ++i_node) {
                         if (node_number == cell_nodes[i_node]) {
                           return i_node;
-                          break;
                         }
                       }
                       return std::numeric_limits<size_t>::max();
@@ -227,31 +225,37 @@ class MeshData
   }
 
  public:
+  PASTIS_INLINE
   const MeshType& mesh() const
   {
     return m_mesh;
   }
 
+  PASTIS_INLINE
   const NodeValuePerCell<const Rd>& Cjr() const
   {
     return m_Cjr;
   }
 
+  PASTIS_INLINE
   const NodeValuePerCell<const double>& ljr() const
   {
     return m_ljr;
   }
 
+  PASTIS_INLINE
   const NodeValuePerCell<const Rd>& njr() const
   {
     return m_njr;
   }
 
+  PASTIS_INLINE
   const CellValue<const Rd>& xj() const
   {
     return m_xj;
   }
 
+  PASTIS_INLINE
   const CellValue<const double>& Vj() const
   {
     return m_Vj;
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 2f83c82dc3a1c5451d51bcd0b293e7d6090a7d4a..ff8ef4a920d0bc2a767b1a0d677c2157a1fe8a9a 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -7,8 +7,7 @@
 #include <Kokkos_Vector.hpp>
 #include <TinyVector.hpp>
 
-#include <RefNodeList.hpp>
-#include <RefFaceList.hpp>
+#include <RefItemList.hpp>
 
 #include <ConnectivityMatrix.hpp>
 #include <IConnectivity.hpp>
@@ -38,7 +37,7 @@ class MeshNodeBoundary
     const auto& face_to_cell_matrix
         = mesh.connectivity().faceToCellMatrix();
 
-    const Array<const FaceId>& face_list = ref_face_list.faceList();
+    const Array<const FaceId>& face_list = ref_face_list.list();
     parallel_for(face_list.size(), PASTIS_LAMBDA(const int& l){
         const auto& face_cells = face_to_cell_matrix[face_list[l]];
         if (face_cells.size()>1) {
@@ -74,7 +73,7 @@ class MeshNodeBoundary
 
   template <typename MeshType>
   MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
-      : m_node_list(ref_node_list.nodeList())
+      : m_node_list(ref_node_list.list())
   {
     static_assert(Dimension == MeshType::Dimension);
   }
@@ -288,7 +287,6 @@ _getNormal(const MeshType& mesh)
       zmax = x;
     }
   }
-#warning re work this part to avoir parallelism dependance
   Array<R3> xmin_array = parallel::allGather(xmin);
   Array<R3> xmax_array = parallel::allGather(xmax);
   Array<R3> ymin_array = parallel::allGather(ymin);
@@ -321,7 +319,6 @@ _getNormal(const MeshType& mesh)
     if (x[2] > zmax[2]) { zmax = x; }
   }
 
-
   const R3 u = xmax-xmin;
   const R3 v = ymax-ymin;
   const R3 w = zmax-zmin;
@@ -356,7 +353,6 @@ _getNormal(const MeshType& mesh)
 
   normal *= 1./sqrt(normal_l2);
 
-#warning Add flatness test
   // this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
 
   return normal;
diff --git a/src/mesh/RefFaceList.hpp b/src/mesh/RefFaceList.hpp
deleted file mode 100644
index cbc7d3b997983c0bdc03d9513684cde8bd669873..0000000000000000000000000000000000000000
--- a/src/mesh/RefFaceList.hpp
+++ /dev/null
@@ -1,40 +0,0 @@
-#ifndef REF_FACE_LIST_HPP
-#define REF_FACE_LIST_HPP
-
-#include <Array.hpp>
-#include <RefId.hpp>
-
-class RefFaceList
-{
- private:
-  RefId m_ref_id;
-  Array<const FaceId> m_face_id_list;
-
- public:
-  const RefId& refId() const
-  {
-    return m_ref_id;
-  }
-
-  const Array<const FaceId>& faceList() const
-  {
-    return m_face_id_list;
-  }
-
-  RefFaceList(const RefId& ref_id,
-              const Array<const FaceId>& face_id_list)
-      : m_ref_id(ref_id),
-        m_face_id_list(face_id_list)
-  {
-    ;
-  }
-
-  RefFaceList& operator=(const RefFaceList&) = default;
-  RefFaceList& operator=(RefFaceList&&) = default;
-
-  RefFaceList() = default;
-  RefFaceList(const RefFaceList&) = default;
-  ~RefFaceList() = default;
-};
-
-#endif // REF_FACE_LIST_HPP
diff --git a/src/mesh/RefItemList.hpp b/src/mesh/RefItemList.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dea69f7ec79e2ca139eeb9f01600e687f6bea3d4
--- /dev/null
+++ b/src/mesh/RefItemList.hpp
@@ -0,0 +1,50 @@
+#ifndef REF_ITEM_LIST_HPP
+#define REF_ITEM_LIST_HPP
+
+#include <Array.hpp>
+#include <RefId.hpp>
+#include <ItemId.hpp>
+
+template <ItemType item_type>
+class RefItemList
+{
+ public:
+  using ItemId = ItemIdT<item_type>;
+
+ private:
+  RefId m_ref_id;
+  Array<const ItemId> m_item_id_list;
+
+ public:
+  const RefId& refId() const
+  {
+    return m_ref_id;
+  }
+
+  const Array<const ItemId>& list() const
+  {
+    return m_item_id_list;
+  }
+
+  RefItemList(const RefId& ref_id,
+              const Array<const ItemId>& item_id_list)
+      : m_ref_id(ref_id),
+        m_item_id_list(item_id_list)
+  {
+    ;
+  }
+
+  RefItemList& operator=(const RefItemList&) = default;
+  RefItemList& operator=(RefItemList&&) = default;
+
+  RefItemList() = default;
+  RefItemList(const RefItemList&) = default;
+  ~RefItemList() = default;
+};
+
+using RefNodeList = RefItemList<ItemType::node>;
+using RefEdgeList = RefItemList<ItemType::edge>;
+using RefFaceList = RefItemList<ItemType::face>;
+using RefCellList = RefItemList<ItemType::cell>;
+
+#endif // REF_ITEM_LIST_HPP
diff --git a/src/mesh/RefNodeList.hpp b/src/mesh/RefNodeList.hpp
deleted file mode 100644
index 462d23bb3d35a6384fcd5176d0584419a586abc3..0000000000000000000000000000000000000000
--- a/src/mesh/RefNodeList.hpp
+++ /dev/null
@@ -1,40 +0,0 @@
-#ifndef REF_NODE_LIST_HPP
-#define REF_NODE_LIST_HPP
-
-#include <Array.hpp>
-#include <RefId.hpp>
-
-class RefNodeList
-{
- private:
-  RefId m_ref_id;
-  Array<const NodeId> m_node_id_list;
-
- public:
-  const RefId& refId() const
-  {
-    return m_ref_id;
-  }
-
-  const Array<const NodeId>& nodeList() const
-  {
-    return m_node_id_list;
-  }
-
-  RefNodeList(const RefId& ref_id,
-              const Array<const NodeId>& node_id_list)
-      : m_ref_id(ref_id),
-        m_node_id_list(node_id_list)
-  {
-    ;
-  }
-
-  RefNodeList& operator=(const RefNodeList&) = default;
-  RefNodeList& operator=(RefNodeList&&) = default;
-
-  RefNodeList() = default;
-  RefNodeList(const RefNodeList&) = default;
-  ~RefNodeList() = default;
-};
-
-#endif // REF_NODE_LIST_HPP
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 821acbd2b8df2a213b78af085d5115e26827a7f8..450554edfc2fc0d3729f9300bb0d5dc516f91f72 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -17,20 +17,10 @@
 
 #include <memory>
 
-// template <typename DataType,
-//           typename ItemOfItem,
-//           typename ConnectivityPtr = std::shared_ptr<const IConnectivity>,
-//           typename Allowed=void>
-// class SubItemValuePerItem;
-
 template <typename DataType,
           typename ItemOfItem,
           typename ConnectivityPtr = std::shared_ptr<const IConnectivity>>
-class SubItemValuePerItem// <DataType,
-                         //  ItemOfItem,
-                         //  ConnectivityPtr// ,
-                         //  // std::enable_if_t<sub_item_type != item_type>
-                         //  >
+class SubItemValuePerItem
 {
  public:
   static constexpr ItemType item_type{ItemOfItem::item_type};
diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
index 0380f97c021b637a623146cd1be7f18524913e2d..669f546c19faf2e0cf059acef0f90ee2083682a4 100644
--- a/src/scheme/BoundaryCondition.hpp
+++ b/src/scheme/BoundaryCondition.hpp
@@ -6,7 +6,7 @@
 
 #include <Array.hpp>
 
-#include <RefNodeList.hpp>
+#include <RefItemList.hpp>
 #include <MeshNodeBoundary.hpp>
 
 class BoundaryCondition
diff --git a/src/utils/Messenger.cpp b/src/utils/Messenger.cpp
index 2bf82efacf7fa363cdf4431b48d13c2906f4deba..6dc4056efd81957f3a5a89fefc328d3d86a1d118 100644
--- a/src/utils/Messenger.cpp
+++ b/src/utils/Messenger.cpp
@@ -26,7 +26,8 @@ void Messenger::destroy()
 }
 
 Messenger::
-Messenger(int& argc, char* argv[])
+Messenger([[maybe_unused]] int& argc,
+          [[maybe_unused]] char* argv[])
 {
 #ifdef PASTIS_HAS_MPI
   MPI_Init(&argc, &argv);
diff --git a/src/utils/Messenger.hpp b/src/utils/Messenger.hpp
index 1575c44410881fe76341c54c1346392321d14fd6..4a54994a273b55e74c253135951a59a95c186909 100644
--- a/src/utils/Messenger.hpp
+++ b/src/utils/Messenger.hpp
@@ -37,7 +37,7 @@ class Messenger
       } else {
         static_assert(std::is_arithmetic_v<DataType>,
                       "Unexpected arithmetic type! Should not occur!");
-        static_assert(not std::is_arithmetic_v<DataType>,
+        static_assert(is_false_v<DataType>,
                       "MPI_Datatype are only defined for arithmetic types!");
         return MPI_Datatype();
       }
@@ -138,7 +138,8 @@ class Messenger
   }
 
   template <typename DataType>
-  void _broadcast_value(DataType& data, const size_t& root_rank) const
+  void _broadcast_value([[maybe_unused]] DataType& data,
+                        [[maybe_unused]] const size_t& root_rank) const
   {
     static_assert(not std::is_const_v<DataType>);
     static_assert(std::is_arithmetic_v<DataType>);
@@ -152,7 +153,8 @@ class Messenger
   }
 
   template <typename ArrayType>
-  void _broadcast_array(ArrayType& array, const size_t& root_rank) const
+  void _broadcast_array([[maybe_unused]] ArrayType& array,
+                        [[maybe_unused]] const size_t& root_rank) const
   {
     using DataType = typename ArrayType::data_type;
     static_assert(not std::is_const_v<DataType>);
@@ -231,12 +233,14 @@ class Messenger
       }
     }
 
-    std::vector<MPI_Status> status_list(request_list.size());
-    if (MPI_SUCCESS != MPI_Waitall(request_list.size(), &(request_list[0]), &(status_list[0]))) {
-      // LCOV_EXCL_START
-      std::cerr << "Communication error!\n";
-      std::exit(1);
-      // LCOV_EXCL_STOP
+    if (request_list.size()>0) {
+      std::vector<MPI_Status> status_list(request_list.size());
+      if (MPI_SUCCESS != MPI_Waitall(request_list.size(), &(request_list[0]), &(status_list[0]))) {
+        // LCOV_EXCL_START
+        std::cerr << "Communication error!\n";
+        std::exit(1);
+        // LCOV_EXCL_STOP
+      }
     }
 
 #else // PASTIS_HAS_MPI
@@ -377,7 +381,7 @@ class Messenger
 
       _allGather(cast_value_array, cast_gather_array);
     } else {
-      static_assert(is_trivially_castable<DataType>, "unexpected type of data");
+      static_assert(is_false_v<DataType>, "unexpected type of data");
     }
     return gather_array;
   }
@@ -401,7 +405,7 @@ class Messenger
 
       _allGather(cast_array, cast_gather_array);
     } else {
-      static_assert(is_trivially_castable<DataType>, "unexpected type of data");
+      static_assert(is_false_v<DataType>, "unexpected type of data");
     }
     return gather_array;
   }
@@ -430,7 +434,7 @@ class Messenger
       auto recv_cast_array = cast_array_to<CastType>::from(recv_array);
       _allToAll(send_cast_array, recv_cast_array);
     } else {
-      static_assert(is_trivially_castable<DataType>, "unexpected type of data");
+      static_assert(is_false_v<DataType>, "unexpected type of data");
     }
     return recv_array;
   }
@@ -453,8 +457,7 @@ class Messenger
         _broadcast_array(cast_array, root_rank);
       }
     } else {
-      static_assert(is_trivially_castable<DataType>,
-                    "unexpected non trivial type of data");
+      static_assert(is_false_v<DataType>, "unexpected type of data");
     }
   }
 
@@ -483,8 +486,7 @@ class Messenger
       auto cast_array = cast_array_to<CastType>::from(array);
       _broadcast_array(cast_array, root_rank);
     } else{
-      static_assert(is_trivially_castable<DataType>,
-                    "unexpected non trivial type of data");
+      static_assert(is_false_v<DataType>, "unexpected type of data");
     }
   }
 
@@ -521,8 +523,7 @@ class Messenger
       using CastType = helper::split_cast_t<DataType>;
       _exchange_through_cast<SendDataType, CastType>(send_array_list, recv_array_list);
     } else {
-      static_assert(is_trivially_castable<RecvDataType>,
-                    "unexpected non trivial type of data");
+      static_assert(is_false_v<RecvDataType>, "unexpected type of data");
     }
   }
 
diff --git a/src/utils/Partitioner.cpp b/src/utils/Partitioner.cpp
index 77839a25c4d2fd91110bcb286dfb17c228024e9a..709038090a5454794bbff49ef9240f2ebd17a780 100644
--- a/src/utils/Partitioner.cpp
+++ b/src/utils/Partitioner.cpp
@@ -60,8 +60,6 @@ Array<int> Partitioner::partition(const CSRGraph& graph)
     part = Array<int>(local_number_of_nodes);
     std::vector<int> vtxdist{0,local_number_of_nodes};
 
-    static_assert(std::is_same<int, int>());
-
     const Array<int>& entries = graph.entries();
     const Array<int>& neighbors = graph.neighbors();
 
@@ -85,7 +83,7 @@ Array<int> Partitioner::partition(const CSRGraph& graph)
 
 #else // PASTIS_HAS_MPI
 
-Array<int> Partitioner::partition(const CSRGraph& graph)
+Array<int> Partitioner::partition(const CSRGraph&)
 {
   return Array<int>(0);
 }
diff --git a/src/utils/PastisTraits.hpp b/src/utils/PastisTraits.hpp
index 4f7877ac9716662363c7b0e1d8051fd78eaec2ac..871e30fd1c471ec67841b156ab89a624a41194b6 100644
--- a/src/utils/PastisTraits.hpp
+++ b/src/utils/PastisTraits.hpp
@@ -19,4 +19,7 @@ inline constexpr bool is_trivially_castable<TinyMatrix<N,T>> = is_trivially_cast
 template <size_t N, typename T>
 inline constexpr bool is_trivially_castable<const TinyMatrix<N,T>> = is_trivially_castable<T>;
 
+template <typename T>
+inline constexpr bool is_false_v = false;
+
 #endif // PASTIS_TRAITS_HPP
diff --git a/src/utils/SignalManager.cpp b/src/utils/SignalManager.cpp
index 319765461db6d8d6c356e21cd2942e8881976934..525db43c8d14a0191f8c9d18790def38051c9761 100644
--- a/src/utils/SignalManager.cpp
+++ b/src/utils/SignalManager.cpp
@@ -36,9 +36,7 @@ std::string SignalManager::signalName(const int& signal)
 void SignalManager::pauseForDebug(const int& signal)
 {
   if (std::string(PASTIS_BUILD_TYPE) != "Release") {
-#warning should try to detect if outputs to a terminal (buggy with mpi)
-    if (// (ConsoleManager::isTerminal(pout()) and (s_pause_on_error == "auto")) or
-        (s_pause_on_error == "yes")) {
+    if (s_pause_on_error == "yes") {
       std::cerr << "\n======================================\n"
 		<< rang::style::reset
 		<< rang::fg::reset