diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 9e6483655340a0ab0e52dbc1fb46c2298ceddea5..ce4c90474faf032790c21539ad090ce29535002b 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -1,10 +1,210 @@
 #include <Connectivity.hpp>
 #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);
@@ -140,16 +340,6 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
     }
     face_to_node_matrix = face_to_node_vector;
   }
-
-  {
-    FaceId l=0;
-    for (const auto& face_cells_vector : face_cells_map) {
-      const Face& face = face_cells_vector.first;
-      m_face_number_map[face] = l;
-      ++l;
-    }
-  }
-#warning check that the number of cell per faces is <=2
 }
 
 template<>
@@ -157,6 +347,7 @@ 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>;
@@ -173,15 +364,6 @@ void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
     }
   }
 
-  {
-    FaceId l=0;
-    for (const auto& face_cells_vector : face_cells_map) {
-      const Face& face = face_cells_vector.first;
-      m_face_number_map[face] = l;
-      ++l;
-    }
-  }
-
   {
     std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
     int l=0;
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index d9789856f84df719d1144b03f33e7af37421d0b9..039b9300748f34da1ce16cd7828b4907c6bec8f2 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -57,203 +57,6 @@ class ConnectivityDescriptor
   ~ConnectivityDescriptor() = default;
 };
 
-template <size_t Dimension>
-class Connectivity;
-
-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 <size_t Dimension>
 class Connectivity final
     : public IConnectivity
@@ -302,10 +105,6 @@ class Connectivity final
 
   CellValue<const double> m_inv_cell_nb_nodes;
 
-  using Face = ConnectivityFace<Dimension>;
-
-  std::unordered_map<Face, FaceId, typename Face::Hash> m_face_number_map;
-
   void _computeCellFaceAndFaceNodeConnectivities();
 
   template <typename SubItemValuePerItemType>
@@ -352,12 +151,6 @@ class Connectivity final
     return m_cell_number;
   }
 
-  PASTIS_INLINE
-  const FaceValue<const int>& faceNumber() const
-  {
-    return m_face_number;
-  }
-
   PASTIS_INLINE
   const EdgeValue<const int>& edgeNumber() const
   {
@@ -695,19 +488,6 @@ class Connectivity final
     return m_inv_cell_nb_nodes;
   }
 
-#warning remove this member
-  unsigned int getFaceNumber(const std::vector<unsigned int>& face_nodes) const
-  {
-    const Face face(face_nodes);
-    auto i_face = m_face_number_map.find(face);
-    if (i_face == m_face_number_map.end()) {
-      perr() << "Face " << face << " not found!\n";
-      throw std::exception();
-      std::exit(0);
-    }
-    return i_face->second;
-  }
-
   Connectivity(const Connectivity&) = delete;
 
   Connectivity(const ConnectivityDescriptor& descriptor);
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index b218aa71a1dfe55e5151bf6c6a52884826c1c509..6412087cf5bc6b7d86a0ed6fd1750e07c5b9035f 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -18,6 +18,7 @@
 
 #include <ArrayUtils.hpp>
 
+#include <unordered_map>
 #include <map>
 #include <regex>
 #include <iomanip>
@@ -81,6 +82,7 @@ class ErrorHandler
     ;
   }
 };
+
 void ErrorHandler::writeErrorMessage() const
 {
   switch(__type) {
@@ -127,6 +129,199 @@ ErrorHandler(const std::string& filename,
   ;
 }
 
+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 GmshReader;
+  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 <int Dimension>
 class MeshDispatcher
 {
@@ -1293,17 +1488,49 @@ GmshReader::__proceedData()
     std::shared_ptr p_connectivity = std::make_shared<Connectivity3D>(descriptor);
     Connectivity3D& connectivity = *p_connectivity;
 
+    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;
+            const auto& face_node_list = connectivity.faceToNodeMatrix();
+            for (FaceId l=0; l<connectivity.numberOfFaces(); ++l) {
+              const auto& node_list = face_node_list[l];
+              std::vector<unsigned int> node_vector(node_list.size());
+              for (size_t r=0; r<node_list.size(); ++r) {
+                node_vector[r] = node_list[r];
+              }
+              face_to_id_map[Face(node_vector)] = l;
+            }
+            return face_to_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_number
-          = connectivity.getFaceNumber({__triangles[f][0], __triangles[f][1], __triangles[f][2]});
+          = [&]{
+              auto i = face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]}));
+              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_number);
     }
     for (unsigned int f=0; f<__quadrangles.size(); ++f) {
       const unsigned int face_number
-          = connectivity.getFaceNumber({__quadrangles[f][0], __quadrangles[f][1],
-                                        __quadrangles[f][2], __quadrangles[f][3]});
+          = [&]{
+              auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
+                                                 __quadrangles[f][2], __quadrangles[f][3]}));
+              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_number);
     }
@@ -1371,9 +1598,29 @@ GmshReader::__proceedData()
     std::shared_ptr p_connectivity = std::make_shared<Connectivity2D>(descriptor);
     Connectivity2D& connectivity = *p_connectivity;
 
+    using Face = ConnectivityFace<2>;
+    const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map
+        = [&]  {
+            std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
+            const auto& face_node_list = connectivity.faceToNodeMatrix();
+            for (FaceId l=0; l<connectivity.numberOfFaces(); ++l) {
+              const auto& node_list = face_node_list[l];
+              face_to_id_map[Face({node_list[0], node_list[1]})] = l;
+            }
+            return face_to_id_map;
+          } ();
+
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
     for (unsigned int e=0; e<__edges.size(); ++e) {
-      const unsigned int edge_number = connectivity.getFaceNumber({__edges[e][0], __edges[e][1]});
+      const unsigned int edge_number
+          = [&]{
+              auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}));
+              if (i == face_to_id_map.end()) {
+                std::cerr << "face " << __edges[e][0] << " not found!\n";
+                std::terminate();
+              }
+              return i->second;
+            }();
       const unsigned int& ref = __edges_ref[e];
       ref_faces_map[ref].push_back(edge_number);
     }