diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index cec70576f019b773fdbc11b3671f7b44885a2a8f..bc8e7093aabf110b3bcebab71fddd5f2585883e5 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -217,7 +217,8 @@ template<size_t Dimension>
 Connectivity<Dimension>::
 Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
              const std::vector<CellType>& cell_type_vector,
-             const std::vector<int>& cell_number_vector)
+             const std::vector<int>& cell_number_vector,
+             const std::vector<int>& node_number_vector)
 {
   Assert(cell_by_node_vector.size() == cell_type_vector.size());
   Assert(cell_number_vector.size() == cell_type_vector.size());
@@ -240,6 +241,12 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
     m_cell_number = cell_number;
   }
 
+  {
+    NodeValue<int> node_number(*this);
+    node_number = convert_to_array(node_number_vector);
+    m_node_number = node_number;
+  }
+
   {
     CellValue<int> cell_global_index(*this);
 #warning index must start accounting number of global indices of other procs
@@ -270,14 +277,17 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
 template Connectivity1D::
 Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
              const std::vector<CellType>& cell_type_vector,
-             const std::vector<int>& cell_number_vector);
+             const std::vector<int>& cell_number_vector,
+             const std::vector<int>& node_number_vector);
 
 template Connectivity2D::
 Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
              const std::vector<CellType>& cell_type_vector,
-             const std::vector<int>& cell_number_vector);
+             const std::vector<int>& cell_number_vector,
+             const std::vector<int>& node_number_vector);
 
 template Connectivity3D::
 Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
              const std::vector<CellType>& cell_type_vector,
-             const std::vector<int>& cell_number_vector);
+             const std::vector<int>& cell_number_vector,
+             const std::vector<int>& node_number_vector);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 8c1dfb6646877de0dbe49e6a03981bb42751192f..f78965cff84b0ba1de9c265ca60e430571036e7c 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -249,6 +249,8 @@ class Connectivity final
   CellValue<const int> m_cell_global_index;
   CellValue<const int> m_cell_number;
 
+  NodeValue<const int> m_node_number;
+
   FaceValuePerCell<const bool> m_cell_face_is_reversed;
 
   NodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
@@ -316,13 +318,19 @@ class Connectivity final
   const CellValue<const CellType>& cellType() const
   {
     return m_cell_type;
-  };
+  }
 
   PASTIS_INLINE
   const CellValue<const int>& cellNumber() const
   {
     return m_cell_number;
-  };
+  }
+
+  PASTIS_INLINE
+  const NodeValue<const int>& nodeNumber() const
+  {
+    return m_node_number;
+  }
 
   PASTIS_INLINE
   const bool& isConnectivityMatrixBuilt(const ItemType& item_type_0,
@@ -653,7 +661,8 @@ class Connectivity final
 
   Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
                const std::vector<CellType>& cell_type_vector,
-               const std::vector<int>& cell_number_vector);
+               const std::vector<int>& cell_number_vector,
+               const std::vector<int>& node_number_vector);
 
   ~Connectivity()
   {
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 4230ea5319f4ddc0208c52eacb67a9942af5adc6..51435e7294c7bb966ceb4b9dd6620bd27f2fab38 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -241,7 +241,7 @@ void GmshReader::_dispatch()
   using MeshType = Mesh<ConnectivityType>;
 
   if (not m_mesh) {
-    std::shared_ptr<ConnectivityType> connectivity(new ConnectivityType({},{},{}));
+    std::shared_ptr<ConnectivityType> connectivity(new ConnectivityType({},{},{},{}));
     NodeValue<Rd> xr;
     m_mesh = std::shared_ptr<IMesh>(new MeshType(connectivity, xr));
   }
@@ -249,11 +249,11 @@ void GmshReader::_dispatch()
 
   MeshDispatcher<Dimension> dispatcher(mesh);
 
-  MeshData data(mesh);
-  data.updateAllData();
+  // MeshData data(mesh);
+  // data.updateAllData();
 
-  std::vector<Array<Rd>> recv_cell_center_by_proc
-      = dispatcher.exchange(data.xj());
+  // std::vector<Array<Rd>> recv_cell_center_by_proc
+  //     = dispatcher.exchange(data.xj());
 
   std::vector<Array<int>> recv_cell_number_by_proc
       = dispatcher.exchange(mesh.connectivity().cellNumber());
@@ -261,17 +261,28 @@ void GmshReader::_dispatch()
   std::vector<Array<CellType>> recv_cell_type_by_proc
       = dispatcher.exchange(mesh.connectivity().cellType());
 
+  CellValue<int> number_of_node_per_cell(mesh.connectivity());
+
+  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+  parallel_for(mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+      number_of_node_per_cell[j] = cell_to_node_matrix[j].size();
+    });
+
+  std::vector<Array<int>> recv_number_of_node_per_cell_by_proc
+      = dispatcher.exchange(number_of_node_per_cell);
+
   for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
     if (parallel::rank() == i_rank) {
       std::cout << "----- rank=" << i_rank << " -----\n";
       for (size_t j_rank=0; j_rank < parallel::size(); ++j_rank) {
         std::cout << "recv from " << j_rank << ':';
         for (size_t i=0; i<recv_cell_number_by_proc[j_rank].size(); ++i) {
-          std::cout << ' ' << recv_cell_number_by_proc[j_rank][i] << '[' << name(recv_cell_type_by_proc[j_rank][i])<< "]:" << recv_cell_center_by_proc[j_rank][i];
+          std::cout << ' ' << recv_cell_number_by_proc[j_rank][i] << '[' << name(recv_cell_type_by_proc[j_rank][i])<< ':' << recv_number_of_node_per_cell_by_proc[j_rank][i] << "] ";
         }
         std::cout << '\n' << std::flush;
       }
     }
+    parallel::barrier();
   }
 
   parallel::Messenger::destroy();
@@ -965,6 +976,8 @@ GmshReader::__proceedData()
     }
   }
 
+  const std::vector<int>& node_number_vector = __verticesNumbers;
+
   TinyVector<15,int> dimension0_mask = zero;
   dimension0_mask[14]=1;
   TinyVector<15,int> dimension1_mask = zero;
@@ -1018,7 +1031,8 @@ GmshReader::__proceedData()
 
     std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector,
                                                                       cell_type_vector,
-                                                                      cell_number_vector));
+                                                                      cell_number_vector,
+                                                                      node_number_vector));
     Connectivity3D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -1085,7 +1099,8 @@ GmshReader::__proceedData()
 
     std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector,
                                                                       cell_type_vector,
-                                                                      cell_number_vector));
+                                                                      cell_number_vector,
+                                                                      node_number_vector));
     Connectivity2D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -1148,7 +1163,8 @@ GmshReader::__proceedData()
 
     std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector,
                                                                       cell_type_vector,
-                                                                      cell_number_vector));
+                                                                      cell_number_vector,
+                                                                      node_number_vector));
     Connectivity1D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;