diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index ce4c90474faf032790c21539ad090ce29535002b..7c7f7c72669a95bd096c697c8aa1929f134dbcda 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -2,402 +2,11 @@
 #include <map>
 
 
-namespace temporary_to_remove
-{
-
-template <size_t Dimension>
-class ConnectivityFace;
-
-template<>
-class ConnectivityFace<1>
-{
- public:
-  friend struct Hash;
-  struct Hash
-  {
-    size_t operator()(const ConnectivityFace& f) const;
-  };
-};
-
-template<>
-class ConnectivityFace<2>
-{
- public:
-  friend struct Hash;
-  struct Hash
-  {
-    size_t operator()(const ConnectivityFace& f) const {
-      size_t hash = 0;
-      hash ^= std::hash<unsigned int>()(f.m_node0_id);
-      hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
-      return hash;
-    }
-  };
-
-  unsigned int m_node0_id;
-  unsigned int m_node1_id;
-
-  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
-  {
-    os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
-    return os;
-  }
-
-  PASTIS_INLINE
-  bool operator==(const ConnectivityFace& f) const
-  {
-    return ((m_node0_id == f.m_node0_id) and
-            (m_node1_id == f.m_node1_id));
-  }
-
-  PASTIS_INLINE
-  bool operator<(const ConnectivityFace& f) const
-  {
-    return ((m_node0_id<f.m_node0_id) or
-            ((m_node0_id == f.m_node0_id) and
-             (m_node1_id<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>& given_node_id_list)
-  {
-    Assert(given_node_id_list.size()==2);
-#warning rework this dirty constructor
-    const auto& [min, max] = std::minmax(given_node_id_list[0], given_node_id_list[1]);
-    m_node0_id = min;
-    m_node1_id = max;
-  }
-
-  PASTIS_INLINE
-  ConnectivityFace(const ConnectivityFace&) = default;
-
-  PASTIS_INLINE
-  ConnectivityFace(ConnectivityFace&&) = default;
-
-  PASTIS_INLINE
-  ~ConnectivityFace() = default;
-};
-
-template <>
-class ConnectivityFace<3>
-{
- private:
-  friend class Connectivity<3>;
-  friend struct Hash;
-  struct Hash
-  {
-    size_t operator()(const ConnectivityFace& f) const {
-      size_t hash = 0;
-      for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
-        hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
-      }
-      return hash;
-    }
-  };
-
-  bool m_reversed;
-  std::vector<NodeId::base_type> m_node_id_list;
-
-  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
-  {
-    for (auto id : f.m_node_id_list) {
-      os << id << ' ';
-    }
-    return os;
-  }
-
-  PASTIS_INLINE
-  const bool& reversed() const
-  {
-    return m_reversed;
-  }
-
-  PASTIS_INLINE
-  const std::vector<unsigned int>& nodeIdList() const
-  {
-    return m_node_id_list;
-  }
-
-  PASTIS_INLINE
-  std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
-  {
-    const auto min_id = std::min_element(node_list.begin(), node_list.end());
-    const int shift = std::distance(node_list.begin(), min_id);
-
-    std::vector<unsigned int> rotated_node_list(node_list.size());
-    if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
-      for (size_t i=0; i<node_list.size(); ++i) {
-        rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
-        m_reversed = true;
-      }
-    } else {
-      for (size_t i=0; i<node_list.size(); ++i) {
-        rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
-      }
-    }
-
-    return rotated_node_list;
-  }
-
-  PASTIS_INLINE
-  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
-      : m_reversed(false),
-        m_node_id_list(_sort(given_node_id_list))
-  {
-    ;
-  }
-
- public:
-  bool operator==(const ConnectivityFace& f) const
-  {
-    if (m_node_id_list.size() == f.nodeIdList().size()) {
-      for (size_t j=0; j<m_node_id_list.size(); ++j) {
-        if (m_node_id_list[j] != f.nodeIdList()[j]) {
-          return false;
-        }
-      }
-      return true;
-    }
-    return false;
-  }
-
-  PASTIS_INLINE
-  bool operator<(const ConnectivityFace& f) const
-  {
-    const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
-    for (size_t i=0; i<min_nb_nodes; ++i) {
-      if (m_node_id_list[i] <  f.m_node_id_list[i]) return true;
-      if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
-    }
-    return m_node_id_list.size() < f.m_node_id_list.size();
-  }
-
-  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;
-
-  PASTIS_INLINE
-  ~ConnectivityFace() = default;
-};
-}
-
-
-template<>
-void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
-{
-  using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
-  using Face = temporary_to_remove::ConnectivityFace<3>;
-
-  const auto& cell_to_node_matrix
-      = this->_getMatrix(ItemType::cell, ItemType::node);
-
-  CellValue<unsigned short> cell_nb_faces(*this);
-  std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
-  for (CellId j=0; j<this->numberOfCells(); ++j) {
-    const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
-
-    switch (m_cell_type[j]) {
-      case CellType::Tetrahedron: {
-        cell_nb_faces[j] = 4;
-        // face 0
-        Face f0({cell_nodes(1),
-                 cell_nodes(2),
-                 cell_nodes(3)});
-        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
-        // face 1
-        Face f1({cell_nodes(0),
-                 cell_nodes(3),
-                 cell_nodes(2)});
-        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
-        // face 2
-        Face f2({cell_nodes(0),
-                 cell_nodes(1),
-                 cell_nodes(3)});
-        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
-        // face 3
-        Face f3({cell_nodes(0),
-                 cell_nodes(2),
-                 cell_nodes(1)});
-        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
-        break;
-      }
-      case CellType::Hexahedron: {
-        // face 0
-        Face f0({cell_nodes(3),
-                 cell_nodes(2),
-                 cell_nodes(1),
-                 cell_nodes(0)});
-        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
-        // face 1
-        Face f1({cell_nodes(4),
-                 cell_nodes(5),
-                 cell_nodes(6),
-                 cell_nodes(7)});
-        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
-        // face 2
-        Face f2({cell_nodes(0),
-                 cell_nodes(4),
-                 cell_nodes(7),
-                 cell_nodes(3)});
-        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
-        // face 3
-        Face f3({cell_nodes(1),
-                 cell_nodes(2),
-                 cell_nodes(6),
-                 cell_nodes(5)});
-        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
-
-        // face 4
-        Face f4({cell_nodes(0),
-                 cell_nodes(1),
-                 cell_nodes(5),
-                 cell_nodes(4)});
-        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
-
-        // face 5
-        Face f5({cell_nodes(3),
-                 cell_nodes(7),
-                 cell_nodes(6),
-                 cell_nodes(2)});
-        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
-
-        cell_nb_faces[j] = 6;
-        break;
-      }
-      default: {
-        perr() << "unexpected cell type!\n";
-        std::exit(0);
-      }
-    }
-  }
-
-  {
-    auto& cell_to_face_matrix
-        = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
-    std::vector<std::vector<unsigned int>> cell_to_face_vector(this->numberOfCells());
-    for (CellId j=0; j<cell_to_face_vector.size(); ++j) {
-      cell_to_face_vector[j].resize(cell_nb_faces[j]);
-    }
-    FaceId l=0;
-    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];
-        cell_to_face_vector[cell_number][cell_local_face] = l;
-      }
-      ++l;
-    }
-    cell_to_face_matrix = cell_to_face_vector;
-  }
-
-  FaceValuePerCell<bool> cell_face_is_reversed(*this);
-  {
-    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];
-        cell_face_is_reversed(cell_number, cell_local_face) = reversed;
-      }
-    }
-
-    m_cell_face_is_reversed = cell_face_is_reversed;
-  }
-
-  {
-    auto& face_to_node_matrix
-        = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
-
-    std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
-    int l=0;
-    for (const auto& face_info : face_cells_map) {
-      const Face& face = face_info.first;
-      face_to_node_vector[l] = face.nodeIdList();
-      ++l;
-    }
-    face_to_node_matrix = face_to_node_vector;
-  }
-}
-
-template<>
-void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
-{
-  const auto& cell_to_node_matrix
-      = this->_getMatrix(ItemType::cell, ItemType::node);
-  using Face = temporary_to_remove::ConnectivityFace<2>;
-
-  // In 2D faces are simply define
-  using CellFaceId = std::pair<CellId, unsigned short>;
-  std::map<Face, std::vector<CellFaceId>> face_cells_map;
-  for (CellId j=0; j<this->numberOfCells(); ++j) {
-    const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
-    for (unsigned short r=0; r<cell_nodes.length; ++r) {
-      NodeId node0_id = cell_nodes(r);
-      NodeId node1_id = cell_nodes((r+1)%cell_nodes.length);
-      if (node1_id<node0_id) {
-        std::swap(node0_id, node1_id);
-      }
-      face_cells_map[Face({node0_id, node1_id})].push_back(std::make_pair(j, r));
-    }
-  }
-
-  {
-    std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
-    int l=0;
-    for (const auto& face_info : face_cells_map) {
-      const Face& face = face_info.first;
-      face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
-      ++l;
-    }
-    auto& face_to_node_matrix
-        = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
-    face_to_node_matrix = face_to_node_vector;
-  }
-
-  {
-    std::vector<std::vector<unsigned int>> face_to_cell_vector(face_cells_map.size());
-    int l=0;
-    for (const auto& face_cells_vector : face_cells_map) {
-      const auto& [face, cell_info_vector] = face_cells_vector;
-      for (const auto& cell_info : cell_info_vector) {
-        face_to_cell_vector[l].push_back(cell_info.first);
-      }
-      ++l;
-    }
-    auto& face_to_cell_matrix
-        = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::cell)];
-    face_to_cell_matrix = face_to_cell_vector;
-  }
-}
-
 template<size_t Dimension>
 Connectivity<Dimension>::
 Connectivity(const ConnectivityDescriptor& descriptor)
 {
-#warning should be checked by
+#warning should be checked by ConnectivityDescriptor
   Assert(descriptor.cell_by_node_vector.size() == descriptor.cell_type_vector.size());
   Assert(descriptor.cell_number_vector.size() == descriptor.cell_type_vector.size());
 
@@ -459,17 +68,43 @@ Connectivity(const ConnectivityDescriptor& descriptor)
   }
 
   if constexpr (Dimension>1) {
-    this->_computeCellFaceAndFaceNodeConnectivities();
+    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;
+
+    if constexpr (Dimension==2) {
+      auto& face_to_cell_matrix
+          = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::cell)];
+      face_to_cell_matrix = descriptor.face_to_cell_vector;
+
+    } else if constexpr (Dimension==3) {
+      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;
+
+      {
+        FaceValuePerCell<bool> cell_face_is_reversed(*this);
+        for (CellId j=0; j<descriptor.cell_face_is_reversed_vector.size(); ++j) {
+          const auto& face_cells_vector = descriptor.cell_face_is_reversed_vector[j];
+          for (unsigned short lj=0; lj<face_cells_vector.size(); ++lj) {
+            cell_face_is_reversed(j,lj) = face_cells_vector[lj];
+          }
+        }
+
+        m_cell_face_is_reversed = cell_face_is_reversed;
+      }
+    }
 
     {
-#warning incorrect in parallel
       FaceValue<int> face_number(*this);
-
-      parallel_for (face_number.size(), PASTIS_LAMBDA(const FaceId& l) {
-          face_number[l] = l;
-        });
+      face_number = convert_to_array(descriptor.face_number_vector);
       m_face_number = face_number;
     }
+    {
+      FaceValue<int> face_owner(*this);
+      face_owner = convert_to_array(descriptor.face_owner_vector);
+      m_face_owner = face_owner;
+    }
   }
 }
 
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 039b9300748f34da1ce16cd7828b4907c6bec8f2..4531a89f2caa20ecdf44f4faff47700ac02242c1 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -42,10 +42,22 @@ class ConnectivityDescriptor
 
  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;
 
   ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
@@ -70,14 +82,19 @@ class Connectivity final
  private:
   ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
   CellValue<const CellType> m_cell_type;
+
+#warning is m_cell_global_index really necessary? should it be computed on demand instead?
   CellValue<const int> m_cell_global_index;
 
   CellValue<const int> m_cell_number;
   FaceValue<const int> m_face_number;
+#warning check that m_edge_number is filled correctly
   EdgeValue<const int> m_edge_number;
   NodeValue<const int> m_node_number;
 
   CellValue<const int> m_cell_owner;
+  FaceValue<const int> m_face_owner;
+#warning Missing EdgeValue<const int> m_edge_owner;
   NodeValue<const int> m_node_owner;
 
   FaceValuePerCell<const bool> m_cell_face_is_reversed;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 6412087cf5bc6b7d86a0ed6fd1750e07c5b9035f..fda25db04737cab3e10d6047870e6ebc377d11c2 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -1168,6 +1168,204 @@ __readPhysicalNames2_2()
   }
 }
 
+void
+GmshReader::__compute2DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
+{
+  // const auto& cell_to_node_matrix
+  //     = this->_getMatrix(ItemType::cell, ItemType::node);
+  using Face = ConnectivityFace<2>;
+
+  // In 2D faces are simply define
+  using CellFaceId = std::pair<CellId, unsigned short>;
+  std::map<Face, std::vector<CellFaceId>> 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 (unsigned short r=0; r<cell_nodes.size(); ++r) {
+      NodeId node0_id = cell_nodes[r];
+      NodeId node1_id = cell_nodes[(r+1)%cell_nodes.size()];
+      if (node1_id<node0_id) {
+        std::swap(node0_id, node1_id);
+      }
+      face_cells_map[Face({node0_id, node1_id})].push_back(std::make_pair(j, r));
+    }
+  }
+
+  {
+    descriptor.face_to_node_vector.resize(face_cells_map.size());
+    int l=0;
+    for (const auto& face_info : face_cells_map) {
+      const Face& face = face_info.first;
+      descriptor.face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
+      ++l;
+    }
+  }
+
+  {
+    descriptor.face_to_cell_vector.resize(face_cells_map.size());
+    int l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& [face, cell_info_vector] = face_cells_vector;
+      for (const auto& cell_info : cell_info_vector) {
+        descriptor.face_to_cell_vector[l].push_back(cell_info.first);
+      }
+      ++l;
+    }
+  }
+
+  {
+    descriptor.face_number_vector.resize(face_cells_map.size());
+    for (size_t l=0; l<face_cells_map.size(); ++l) {
+      descriptor.face_number_vector[l] = l;
+    }
+    perr() << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red << "must use provided faces number" << rang::fg::reset << '\n';
+  }
+}
+
+void
+GmshReader::__compute3DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
+{
+  using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
+  using Face = ConnectivityFace<3>;
+
+  // const auto& cell_to_node_matrix
+  //     = this->_getMatrix(ItemType::cell, ItemType::node);
+
+  Array<unsigned short> cell_nb_faces(descriptor.cell_by_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];
+
+    switch (descriptor.cell_type_vector[j]) {
+      case CellType::Tetrahedron: {
+        cell_nb_faces[j] = 4;
+        // face 0
+        Face f0({cell_nodes[1],
+                 cell_nodes[2],
+                 cell_nodes[3]});
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes[0],
+                 cell_nodes[3],
+                 cell_nodes[2]});
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes[0],
+                 cell_nodes[1],
+                 cell_nodes[3]});
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+        // face 3
+        Face f3({cell_nodes[0],
+                 cell_nodes[2],
+                 cell_nodes[1]});
+        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+        break;
+      }
+      case CellType::Hexahedron: {
+        // face 0
+        Face f0({cell_nodes[3],
+                 cell_nodes[2],
+                 cell_nodes[1],
+                 cell_nodes[0]});
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes[4],
+                 cell_nodes[5],
+                 cell_nodes[6],
+                 cell_nodes[7]});
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes[0],
+                 cell_nodes[4],
+                 cell_nodes[7],
+                 cell_nodes[3]});
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+        // face 3
+        Face f3({cell_nodes[1],
+                 cell_nodes[2],
+                 cell_nodes[6],
+                 cell_nodes[5]});
+        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+
+        // face 4
+        Face f4({cell_nodes[0],
+                 cell_nodes[1],
+                 cell_nodes[5],
+                 cell_nodes[4]});
+        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+
+        // face 5
+        Face f5({cell_nodes[3],
+                 cell_nodes[7],
+                 cell_nodes[6],
+                 cell_nodes[2]});
+        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+
+        cell_nb_faces[j] = 6;
+        break;
+      }
+      default: {
+        perr() << "unexpected cell type!\n";
+        std::exit(0);
+      }
+    }
+  }
+
+  {
+    descriptor.cell_to_face_vector.resize(descriptor.cell_by_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]);
+    }
+    FaceId l=0;
+    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];
+        descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
+      }
+      ++l;
+    }
+  }
+
+  {
+    descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_by_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]);
+    }
+    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];
+        descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
+      }
+    }
+  }
+
+  {
+    descriptor.face_to_node_vector.resize(face_cells_map.size());
+    int l=0;
+    for (const auto& face_info : face_cells_map) {
+      const Face& face = face_info.first;
+      descriptor.face_to_node_vector[l] = face.nodeIdList();
+      ++l;
+    }
+  }
+
+  {
+    descriptor.face_number_vector.resize(face_cells_map.size());
+    for (size_t l=0; l<face_cells_map.size(); ++l) {
+      descriptor.face_number_vector[l] = l;
+    }
+    perr() << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red << "must use provided faces number" << rang::fg::reset << '\n';
+  }
+}
+
+
 void
 GmshReader::__proceedData()
 {
@@ -1277,7 +1475,6 @@ GmshReader::__proceedData()
   __verticesCorrepondance.resize(maxNumber+1,-1);
 
   for (size_t i=0; i<__verticesNumbers.size(); ++i) {
-
     __verticesCorrepondance[__verticesNumbers[i]] = i;
   }
 
@@ -1475,11 +1672,18 @@ GmshReader::__proceedData()
       descriptor.cell_number_vector[jh] = __hexahedra_number[j];
     }
 
+    __compute3DCellFaceAndFaceNodeConnectivities(descriptor);
+
     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(),
@@ -1585,11 +1789,18 @@ GmshReader::__proceedData()
       descriptor.cell_number_vector[jq] = __quadrangles_number[j];
     }
 
+    this->__compute2DCellFaceAndFaceNodeConnectivities(descriptor);
+
     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(),
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index dc1aa9914c114574de4f435b572384c12e360c20..cd7d89f236272cdf6c584a270587677c1339f379 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -14,6 +14,8 @@
 #include <map>
 #include <fstream>
 
+class ConnectivityDescriptor;
+
 class GmshReader
 {
 public:
@@ -256,6 +258,8 @@ private:
    *
    */
   void __readPhysicalNames2_2();
+  void __compute2DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
+  void __compute3DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
 
   template <int Dimension>
   void _dispatch();