diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index b9a65729637128081a7206ad6fa5d98ad95e9978..0181e5b2cf87469c8481a8ec5b3b70414473a321 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -9,7 +9,7 @@
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/CartesianMeshBuilder.hpp>
 #include <mesh/Connectivity.hpp>
-#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
 #include <mesh/GmshReader.hpp>
 #include <mesh/Mesh.hpp>
 #include <utils/Exceptions.hpp>
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 46a16a505eb85683097e4f62de567ed4ef91633c..97b166608fae5c7bc454dfe063bf60c4e6f18ba2 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -3,7 +3,7 @@
 add_library(
   PugsMesh
   CartesianMeshBuilder.cpp
-  DiamondDualMeshBuilder.cpp
+  DiamondDualConnectivityBuilder.cpp
   Connectivity.cpp
   ConnectivityComputer.cpp
   ConnectivityDispatcher.cpp
diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
index a272dde4d93e7e164f129fcab543bf4ef03ed118..caff700fab42228f8aa2a1cd0b246e736881de0e 100644
--- a/src/mesh/CartesianMeshBuilder.cpp
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -7,16 +7,171 @@
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
 
+template <size_t Dimension>
+NodeValue<TinyVector<Dimension>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<Dimension>&,
+                                          const TinyVector<Dimension>&,
+                                          const TinyVector<Dimension, uint64_t>&,
+                                          const IConnectivity&) const
+{
+  static_assert(Dimension <= 3, "invalid dimension");
+}
+
+template <>
+NodeValue<TinyVector<1>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<1>& a,
+                                          const TinyVector<1>& b,
+                                          const TinyVector<1, uint64_t>& cell_size,
+                                          const IConnectivity& connectivity) const
+{
+  const double h = (b[0] - a[0]) / cell_size[0];
+  NodeValue<TinyVector<1>> xr(connectivity);
+  parallel_for(
+    connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId r) { xr[r] = a + r * h; });
+
+  return xr;
+}
+
+template <>
+NodeValue<TinyVector<2>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<2>& a,
+                                          const TinyVector<2>& b,
+                                          const TinyVector<2, uint64_t>& cell_size,
+                                          const IConnectivity& connectivity) const
+{
+  const TinyVector<2> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1]};
+
+  const TinyVector<2, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
+
+  const auto node_logic_id = [&](size_t r) {
+    const uint64_t r0 = r / node_size[1];
+    const uint64_t r1 = r % node_size[1];
+    return TinyVector<2, uint64_t>{r0, r1};
+  };
+
+  NodeValue<TinyVector<2>> xr(connectivity);
+
+  parallel_for(
+    connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+      const TinyVector<2, uint64_t> node_index = node_logic_id(r);
+      for (size_t i = 0; i < 2; ++i) {
+        xr[r][i] = a[i] + node_index[i] * h[i];
+      }
+    });
+
+  return xr;
+}
+
+template <>
+NodeValue<TinyVector<3>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<3>& a,
+                                          const TinyVector<3>& b,
+                                          const TinyVector<3, uint64_t>& cell_size,
+                                          const IConnectivity& connectivity) const
+{
+  const TinyVector<3, uint64_t> node_size = [&] {
+    TinyVector node_size{cell_size};
+    for (size_t i = 0; i < 3; ++i) {
+      node_size[i] += 1;
+    }
+    return node_size;
+  }();
+
+  const auto node_logic_id = [&](size_t r) {
+    const size_t slice1  = node_size[1] * node_size[2];
+    const size_t& slice2 = node_size[2];
+    const uint64_t r0    = r / slice1;
+    const uint64_t r1    = (r - r0 * slice1) / slice2;
+    const uint64_t r2    = r - (r0 * slice1 + r1 * slice2);
+    return TinyVector<3, uint64_t>{r0, r1, r2};
+  };
+
+  const TinyVector<3> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1], (b[2] - a[2]) / cell_size[2]};
+
+  NodeValue<TinyVector<3>> xr(connectivity);
+  parallel_for(
+    connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+      const TinyVector<3, uint64_t> node_index = node_logic_id(r);
+      for (size_t i = 0; i < 3; ++i) {
+        xr[r][i] = a[i] + node_index[i] * h[i];
+      }
+    });
+
+  return xr;
+}
+
+template <size_t Dimension>
+void
+CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<Dimension>& a,
+                                          const TinyVector<Dimension>& b,
+                                          const TinyVector<Dimension, uint64_t>& cell_size)
+{
+  static_assert(Dimension <= 3, "unexpected dimension");
+
+  LogicalConnectivityBuilder logical_connectivity_builder{cell_size};
+
+  using ConnectivityType = Connectivity<Dimension>;
+
+  std::shared_ptr<const ConnectivityType> p_connectivity =
+    std::dynamic_pointer_cast<const ConnectivityType>(logical_connectivity_builder.connectivity());
+  const ConnectivityType& connectivity = *p_connectivity;
+
+  NodeValue<TinyVector<Dimension>> xr = _getNodeCoordinates(a, b, cell_size, connectivity);
+
+  m_mesh = std::make_shared<Mesh<ConnectivityType>>(p_connectivity, xr);
+}
+
+template <size_t Dimension>
+CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<Dimension>& a,
+                                           const TinyVector<Dimension>& b,
+                                           const TinyVector<Dimension, uint64_t>& size)
+{
+  if (parallel::rank() == 0) {
+    TinyVector lenght = b - a;
+    for (size_t i = 0; i < Dimension; ++i) {
+      if (lenght[i] == 0) {
+        throw NormalError("invalid box definition corners share a component");
+      }
+    }
+
+    TinyVector<Dimension> corner0 = a;
+    TinyVector<Dimension> corner1 = b;
+
+    for (size_t i = 0; i < Dimension; ++i) {
+      if (corner0[i] > corner1[i]) {
+        std::swap(corner0[i], corner1[i]);
+      }
+    }
+
+    this->_buildCartesianMesh(corner0, corner1, size);
+  }
+  this->_dispatch<Dimension>();
+}
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<1>&,
+                                                    const TinyVector<1>&,
+                                                    const TinyVector<1, uint64_t>&);
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<2>&,
+                                                    const TinyVector<2>&,
+                                                    const TinyVector<2, uint64_t>&);
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<3>&,
+                                                    const TinyVector<3>&,
+                                                    const TinyVector<3, uint64_t>&);
+
+/// =====================
+
 template <size_t Dimension>
 inline void
-CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+LogicalConnectivityBuilder::_buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
 {
   static_assert(Dimension <= 3, "unexpected dimension");
 }
 
 template <>
 inline void
-CartesianMeshBuilder::_buildBoundaryNodeList(
+LogicalConnectivityBuilder::_buildBoundaryNodeList(
   const TinyVector<1, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
   ConnectivityDescriptor& descriptor)
 {
@@ -35,7 +190,7 @@ CartesianMeshBuilder::_buildBoundaryNodeList(
 
 template <>
 void
-CartesianMeshBuilder::_buildBoundaryNodeList(
+LogicalConnectivityBuilder::_buildBoundaryNodeList(
   const TinyVector<2, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
   ConnectivityDescriptor& descriptor)
 {
@@ -72,8 +227,8 @@ CartesianMeshBuilder::_buildBoundaryNodeList(
 
 template <>
 void
-CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<3, uint64_t>& cell_size,
-                                             ConnectivityDescriptor& descriptor)
+LogicalConnectivityBuilder::_buildBoundaryNodeList(const TinyVector<3, uint64_t>& cell_size,
+                                                   ConnectivityDescriptor& descriptor)
 {
   const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
 
@@ -132,15 +287,15 @@ CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<3, uint64_t>& cell
 
 template <size_t Dimension>
 void
-CartesianMeshBuilder::_buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
 {
   static_assert(Dimension == 3, "unexpected dimension");
 }
 
 template <>
 void
-CartesianMeshBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
-                                             ConnectivityDescriptor& descriptor)
+LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
+                                                   ConnectivityDescriptor& descriptor)
 {
   using Edge                                                                 = ConnectivityFace<2>;
   const auto& node_number_vector                                             = descriptor.node_number_vector;
@@ -242,14 +397,14 @@ CartesianMeshBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell
 
 template <size_t Dimension>
 void
-CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
 {
   static_assert(Dimension >= 2 and Dimension <= 3, "unexpected dimension");
 }
 
 template <>
 void
-CartesianMeshBuilder::_buildBoundaryFaceList(
+LogicalConnectivityBuilder::_buildBoundaryFaceList(
   const TinyVector<2, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
   ConnectivityDescriptor& descriptor)
 {
@@ -316,8 +471,8 @@ CartesianMeshBuilder::_buildBoundaryFaceList(
 
 template <>
 void
-CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
-                                             ConnectivityDescriptor& descriptor)
+LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
+                                                   ConnectivityDescriptor& descriptor)
 {
   using Face = ConnectivityFace<3>;
 
@@ -425,9 +580,7 @@ CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell
 
 template <>
 void
-CartesianMeshBuilder::_buildCartesianMesh(
-  const TinyVector<1>& a,                     // clazy:exclude=function-args-by-value
-  const TinyVector<1>& b,                     // clazy:exclude=function-args-by-value
+LogicalConnectivityBuilder::_buildConnectivity(
   const TinyVector<1, uint64_t>& cell_size)   // clazy:exclude=function-args-by-value
 {
   const size_t number_of_cells = cell_size[0];
@@ -465,24 +618,12 @@ CartesianMeshBuilder::_buildCartesianMesh(
   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;
-
-  const double h = (b[0] - a[0]) / cell_size[0];
-
-  NodeValue<TinyVector<1>> xr(connectivity);
-  for (NodeId r = 0; r < number_of_nodes; ++r) {
-    xr[r] = a + r * h;
-  }
-
-  m_mesh = std::make_shared<Mesh<Connectivity1D>>(p_connectivity, xr);
+  m_connectivity = Connectivity1D::build(descriptor);
 }
 
 template <>
 void
-CartesianMeshBuilder::_buildCartesianMesh(
-  const TinyVector<2>& a,                     // clazy:exclude=function-args-by-value
-  const TinyVector<2>& b,                     // clazy:exclude=function-args-by-value
+LogicalConnectivityBuilder::_buildConnectivity(
   const TinyVector<2, uint64_t>& cell_size)   // clazy:exclude=function-args-by-value
 {
   constexpr size_t Dimension = 2;
@@ -506,12 +647,6 @@ CartesianMeshBuilder::_buildCartesianMesh(
   descriptor.cell_type_vector.resize(number_of_cells);
   std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Quadrangle);
 
-  const auto node_logic_id = [&](size_t r) {
-    const uint64_t r0 = r / node_size[1];
-    const uint64_t r1 = r % node_size[1];
-    return TinyVector<Dimension, uint64_t>{r0, r1};
-  };
-
   const auto node_number = [&](const TinyVector<Dimension, uint64_t> node_logic_id) {
     return node_logic_id[0] * node_size[1] + node_logic_id[1];
   };
@@ -535,7 +670,7 @@ CartesianMeshBuilder::_buildCartesianMesh(
     }
   }
 
-  MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
 
   this->_buildBoundaryNodeList(cell_size, descriptor);
   this->_buildBoundaryFaceList(cell_size, descriptor);
@@ -549,27 +684,12 @@ CartesianMeshBuilder::_buildCartesianMesh(
   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        = Connectivity<Dimension>::build(descriptor);
-  Connectivity<Dimension>& connectivity = *p_connectivity;
-
-  const TinyVector<Dimension> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1]};
-
-  NodeValue<TinyVector<Dimension>> xr(connectivity);
-  for (NodeId r = 0; r < number_of_nodes; ++r) {
-    const TinyVector<Dimension, uint64_t> node_index = node_logic_id(r);
-    for (size_t i = 0; i < Dimension; ++i) {
-      xr[r][i] = a[i] + node_index[i] * h[i];
-    }
-  }
-
-  m_mesh = std::make_shared<Mesh<Connectivity<Dimension>>>(p_connectivity, xr);
+  m_connectivity = Connectivity<Dimension>::build(descriptor);
 }
 
 template <>
 void
-CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
-                                          const TinyVector<3>& b,
-                                          const TinyVector<3, uint64_t>& cell_size)
+LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& cell_size)
 {
   constexpr size_t Dimension = 3;
 
@@ -618,15 +738,6 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
     return TinyVector<Dimension, uint64_t>{j0, j1, j2};
   };
 
-  const auto node_logic_id = [&](size_t r) {
-    const size_t slice1  = node_size[1] * node_size[2];
-    const size_t& slice2 = node_size[2];
-    const uint64_t r0    = r / slice1;
-    const uint64_t r1    = (r - r0 * slice1) / slice2;
-    const uint64_t r2    = r - (r0 * slice1 + r1 * slice2);
-    return TinyVector<Dimension, uint64_t>{r0, r1, r2};
-  };
-
   const auto node_number = [&](const TinyVector<Dimension, uint64_t>& node_logic_id) {
     return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
   };
@@ -649,8 +760,8 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
     }
   }
 
-  MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
-  MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
+  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+  ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
 
   this->_buildBoundaryNodeList(cell_size, descriptor);
   this->_buildBoundaryEdgeList(cell_size, descriptor);
@@ -668,58 +779,19 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
   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        = Connectivity<Dimension>::build(descriptor);
-  Connectivity<Dimension>& connectivity = *p_connectivity;
-
-  const TinyVector<Dimension> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1],
-                                (b[2] - a[2]) / cell_size[2]};
-
-  NodeValue<TinyVector<Dimension>> xr(connectivity);
-  for (NodeId r = 0; r < number_of_nodes; ++r) {
-    const TinyVector<Dimension, uint64_t> node_index = node_logic_id(r);
-    for (size_t i = 0; i < Dimension; ++i) {
-      xr[r][i] = a[i] + node_index[i] * h[i];
-    }
-  }
-
-  m_mesh = std::make_shared<Mesh<Connectivity<Dimension>>>(p_connectivity, xr);
+  m_connectivity = Connectivity<Dimension>::build(descriptor);
 }
 
 template <size_t Dimension>
-CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<Dimension>& a,
-                                           const TinyVector<Dimension>& b,
-                                           const TinyVector<Dimension, uint64_t>& size)
+LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size)
 {
   if (parallel::rank() == 0) {
-    TinyVector lenght = b - a;
-    for (size_t i = 0; i < Dimension; ++i) {
-      if (lenght[i] == 0) {
-        throw NormalError("invalid box definition corners share a component");
-      }
-    }
-
-    TinyVector<Dimension> corner0 = a;
-    TinyVector<Dimension> corner1 = b;
-
-    for (size_t i = 0; i < Dimension; ++i) {
-      if (corner0[i] > corner1[i]) {
-        std::swap(corner0[i], corner1[i]);
-      }
-    }
-
-    this->_buildCartesianMesh(corner0, corner1, size);
+    this->_buildConnectivity(size);
   }
-  this->_dispatch<Dimension>();
 }
 
-template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<1>&,
-                                                    const TinyVector<1>&,
-                                                    const TinyVector<1, uint64_t>&);
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<1, uint64_t>&);
 
-template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<2>&,
-                                                    const TinyVector<2>&,
-                                                    const TinyVector<2, uint64_t>&);
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<2, uint64_t>&);
 
-template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<3>&,
-                                                    const TinyVector<3>&,
-                                                    const TinyVector<3, uint64_t>&);
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<3, uint64_t>&);
diff --git a/src/mesh/CartesianMeshBuilder.hpp b/src/mesh/CartesianMeshBuilder.hpp
index 712511583ab1e089c19ed29f6e6b50d2d3ab7be9..9ed1ea81cd4c3dbf0863fb50fe423cd49aee2d0b 100644
--- a/src/mesh/CartesianMeshBuilder.hpp
+++ b/src/mesh/CartesianMeshBuilder.hpp
@@ -8,6 +8,28 @@
 #include <memory>
 
 class CartesianMeshBuilder : public MeshBuilderBase
+{
+ private:
+  template <size_t Dimension>
+  void _buildCartesianMesh(const TinyVector<Dimension>& a,
+                           const TinyVector<Dimension>& b,
+                           const TinyVector<Dimension, uint64_t>& size);
+
+  template <size_t Dimension>
+  NodeValue<TinyVector<Dimension>> _getNodeCoordinates(const TinyVector<Dimension>& a,
+                                                       const TinyVector<Dimension>& b,
+                                                       const TinyVector<Dimension, uint64_t>& size,
+                                                       const IConnectivity& connectivity) const;
+
+ public:
+  template <size_t Dimension>
+  CartesianMeshBuilder(const TinyVector<Dimension>& a,
+                       const TinyVector<Dimension>& b,
+                       const TinyVector<Dimension, uint64_t>& size);
+  ~CartesianMeshBuilder() = default;
+};
+
+class LogicalConnectivityBuilder : public ConnectivityBuilderBase
 {
  private:
   template <size_t Dimension>
@@ -20,16 +42,12 @@ class CartesianMeshBuilder : public MeshBuilderBase
   void _buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
 
   template <size_t Dimension>
-  void _buildCartesianMesh(const TinyVector<Dimension>& a,
-                           const TinyVector<Dimension>& b,
-                           const TinyVector<Dimension, uint64_t>& size);
+  void _buildConnectivity(const TinyVector<Dimension, uint64_t>& size);
 
  public:
   template <size_t Dimension>
-  CartesianMeshBuilder(const TinyVector<Dimension>& a,
-                       const TinyVector<Dimension>& b,
-                       const TinyVector<Dimension, uint64_t>& size);
-  ~CartesianMeshBuilder() = default;
+  LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size);
+  ~LogicalConnectivityBuilder() = default;
 };
 
 #endif   // CARTESIAN_MESH_BUILDER_HPP
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
similarity index 83%
rename from src/mesh/DiamondDualMeshBuilder.cpp
rename to src/mesh/DiamondDualConnectivityBuilder.cpp
index 4c8c0658036bcd6f366897d641ba56bc784e046a..585982245f568a590f78217996c4fe64dc8eda5f 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -1,4 +1,4 @@
-#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
 
 #include <mesh/Connectivity.hpp>
 #include <mesh/ConnectivityDescriptor.hpp>
@@ -12,8 +12,8 @@
 
 template <size_t Dimension>
 void
-DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
-                                                            ConnectivityDescriptor& diamond_descriptor)
+DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
+                                                                    ConnectivityDescriptor& diamond_descriptor)
 {
   const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
   const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
@@ -154,8 +154,8 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
 
 template <>
 void
-DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivity<1>& primal_connectivity,
-                                                               ConnectivityDescriptor& diamond_descriptor)
+DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivity<1>& primal_connectivity,
+                                                                       ConnectivityDescriptor& diamond_descriptor)
 {
   const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
   const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
@@ -227,7 +227,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivit
 
 template <size_t Dimension>
 void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
+DiamondDualConnectivityBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
 {
   static_assert(Dimension <= 3, "invalid mesh dimension");
   using ConnectivityType = Connectivity<Dimension>;
@@ -242,9 +242,9 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
   this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor);
 
   if constexpr (Dimension > 1) {
-    MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
+    ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
     if constexpr (Dimension > 2) {
-      MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
+      ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
     }
   }
 
@@ -474,57 +474,128 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
     }
   }
 
-  std::shared_ptr p_diamond_connectivity = ConnectivityType::build(diamond_descriptor);
-  ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
+  m_connectivity = ConnectivityType::build(diamond_descriptor);
+}
+
+DiamondDualConnectivityBuilder::DiamondDualConnectivityBuilder(const std::shared_ptr<const IMesh>& p_mesh)
+{
+  switch (p_mesh->dimension()) {
+  case 1: {
+    this->_buildDiamondMeshFrom<1>(p_mesh);
+    break;
+  }
+  case 2: {
+    this->_buildDiamondMeshFrom<2>(p_mesh);
+    break;
+  }
+  case 3: {
+    this->_buildDiamondMeshFrom<3>(p_mesh);
+    break;
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension: " + std::to_string(p_mesh->dimension()));
+  }
+  }
+}
+
+template <size_t Dimension>
+void
+DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(
+  const Mesh<Connectivity<Dimension>>& primal_mesh,
+  const std::shared_ptr<const Connectivity<Dimension>>& p_diamond_connectivity)
+{
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<Connectivity<Dimension>>;
+
+  const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
 
   NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
 
-  const auto primal_xr = primal_mesh.xr();
+  const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
   MeshData<MeshType> primal_mesh_data{primal_mesh};
-  const auto primal_xj = primal_mesh_data.xj();
+  const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
 
-  {
-    if constexpr (Dimension == 1) {
-      const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
-      NodeId next_node_id             = 0;
-      for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
-        if (node_to_cell_matrix[primal_node_id].size() == 1) {
-          diamond_xr[next_node_id++] = primal_xr[primal_node_id];
-        }
-      }
+  NodeId i_node = 0;
+  for (; i_node < primal_mesh.numberOfNodes(); ++i_node) {
+    diamond_xr[i_node] = primal_xr[i_node];
+  }
 
-      for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) {
-        diamond_xr[next_node_id++] = primal_xj[i_cell];
-      }
+  for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
+    diamond_xr[i_node++] = primal_xj[i_cell];
+  }
 
-    } else {
-      NodeId i_node = 0;
-      for (; i_node < primal_number_of_nodes; ++i_node) {
-        diamond_xr[i_node] = primal_xr[i_node];
-      }
+  m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
+}
 
-      for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) {
-        diamond_xr[i_node++] = primal_xj[i_cell];
-      }
+template <>
+void
+DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const Mesh<Connectivity<1>>& primal_mesh,
+                                                  const std::shared_ptr<const Connectivity<1>>& p_diamond_connectivity)
+{
+  using ConnectivityType = Connectivity<1>;
+  using MeshType         = Mesh<Connectivity<1>>;
+
+  const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
+  const auto& node_to_cell_matrix             = primal_connectivity.nodeToCellMatrix();
+
+  const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
+
+  NodeValue<TinyVector<1>> diamond_xr{diamond_connectivity};
+
+  const NodeValue<const TinyVector<1>> primal_xr = primal_mesh.xr();
+  MeshData<MeshType> primal_mesh_data{primal_mesh};
+  const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
+
+  NodeId next_node_id = 0;
+  for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+    if (node_to_cell_matrix[primal_node_id].size() == 1) {
+      diamond_xr[next_node_id++] = primal_xr[primal_node_id];
     }
   }
 
+  for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
+    diamond_xr[next_node_id++] = primal_xj[i_cell];
+  }
+
   m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
 }
 
 DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
 {
+  DiamondDualConnectivityBuilder builder(p_mesh);
+
   switch (p_mesh->dimension()) {
   case 1: {
-    this->_buildDiamondMeshFrom<1>(p_mesh);
+    using ConnectivityType = Connectivity<1>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    std::shared_ptr p_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+    const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+    this->_buildDualDiamondMeshFrom(mesh, p_connectivity);
     break;
   }
   case 2: {
-    this->_buildDiamondMeshFrom<2>(p_mesh);
+    using ConnectivityType = Connectivity<2>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    std::shared_ptr p_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+    const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+    this->_buildDualDiamondMeshFrom(mesh, p_connectivity);
     break;
   }
   case 3: {
-    this->_buildDiamondMeshFrom<3>(p_mesh);
+    using ConnectivityType = Connectivity<3>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    std::shared_ptr p_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+    const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+    this->_buildDualDiamondMeshFrom(mesh, p_connectivity);
     break;
   }
   default: {
diff --git a/src/mesh/DiamondDualConnectivityBuilder.hpp b/src/mesh/DiamondDualConnectivityBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a2fe946b934d5a1e620023f14ffbbab643df3dd1
--- /dev/null
+++ b/src/mesh/DiamondDualConnectivityBuilder.hpp
@@ -0,0 +1,43 @@
+#ifndef DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
+#define DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
+
+#include <mesh/IMesh.hpp>
+#include <mesh/MeshBuilderBase.hpp>
+
+#include <memory>
+
+template <size_t>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+class ConnectivityDescriptor;
+
+class DiamondDualConnectivityBuilder : public ConnectivityBuilderBase
+{
+ private:
+  template <size_t Dimension>
+  void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
+
+  template <size_t Dimension>
+  void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
+
+ public:
+  DiamondDualConnectivityBuilder(const std::shared_ptr<const IMesh>&);
+  ~DiamondDualConnectivityBuilder() = default;
+};
+
+class DiamondDualMeshBuilder : public MeshBuilderBase
+{
+ private:
+  template <size_t Dimension>
+  void _buildDualDiamondMeshFrom(const Mesh<Connectivity<Dimension>>&,
+                                 const std::shared_ptr<const Connectivity<Dimension>>&);
+
+ public:
+  DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
+  ~DiamondDualMeshBuilder() = default;
+};
+
+#endif   // DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
deleted file mode 100644
index a6af1b830c3d86738eb18e3d4c4b73a76a287247..0000000000000000000000000000000000000000
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ /dev/null
@@ -1,28 +0,0 @@
-#ifndef DIAMOND_DUAL_MESH_BUILDER_HPP
-#define DIAMOND_DUAL_MESH_BUILDER_HPP
-
-#include <mesh/IMesh.hpp>
-#include <mesh/MeshBuilderBase.hpp>
-
-#include <memory>
-
-template <size_t>
-class Connectivity;
-
-class ConnectivityDescriptor;
-
-class DiamondDualMeshBuilder : public MeshBuilderBase
-{
- private:
-  template <size_t Dimension>
-  void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
-
-  template <size_t Dimension>
-  void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
-
- public:
-  DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
-  ~DiamondDualMeshBuilder() = default;
-};
-
-#endif   // DIAMOND_DUAL_MESH_BUILDER_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index e99c04722bff61b143d7fbf7c6b09ee881136320..6ca41953f80673e9784ef6f8e238ddaa1d8e101f 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -23,6 +23,450 @@
 #include <sstream>
 #include <unordered_map>
 
+template <size_t Dimension>
+class GmshConnectivityBuilder : public ConnectivityBuilderBase
+{
+ public:
+  GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells);
+};
+
+template <>
+GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells)
+{
+  ConnectivityDescriptor descriptor;
+
+  descriptor.node_number_vector = gmsh_data.__verticesNumbers;
+  descriptor.cell_type_vector.resize(nb_cells);
+  descriptor.cell_number_vector.resize(nb_cells);
+  descriptor.cell_to_node_vector.resize(nb_cells);
+
+  for (size_t j = 0; j < nb_cells; ++j) {
+    descriptor.cell_to_node_vector[j].resize(2);
+    for (int r = 0; r < 2; ++r) {
+      descriptor.cell_to_node_vector[j][r] = gmsh_data.__edges[j][r];
+    }
+    descriptor.cell_type_vector[j]   = CellType::Line;
+    descriptor.cell_number_vector[j] = gmsh_data.__edges_number[j];
+  }
+
+  std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+  for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
+    const unsigned int point_number = gmsh_data.__points[r];
+    const unsigned int& ref         = gmsh_data.__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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.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.node_owner_vector.resize(descriptor.node_number_vector.size());
+  std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+
+  m_connectivity = Connectivity1D::build(descriptor);
+}
+
+template <>
+GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells)
+{
+  ConnectivityDescriptor descriptor;
+
+  descriptor.node_number_vector = gmsh_data.__verticesNumbers;
+  descriptor.cell_type_vector.resize(nb_cells);
+  descriptor.cell_number_vector.resize(nb_cells);
+  descriptor.cell_to_node_vector.resize(nb_cells);
+
+  const size_t nb_triangles = gmsh_data.__triangles.size();
+  for (size_t j = 0; j < nb_triangles; ++j) {
+    descriptor.cell_to_node_vector[j].resize(3);
+    for (int r = 0; r < 3; ++r) {
+      descriptor.cell_to_node_vector[j][r] = gmsh_data.__triangles[j][r];
+    }
+    descriptor.cell_type_vector[j]   = CellType::Triangle;
+    descriptor.cell_number_vector[j] = gmsh_data.__triangles_number[j];
+  }
+
+  const size_t nb_quadrangles = gmsh_data.__quadrangles.size();
+  for (size_t j = 0; j < nb_quadrangles; ++j) {
+    const size_t jq = j + nb_triangles;
+    descriptor.cell_to_node_vector[jq].resize(4);
+    for (int r = 0; r < 4; ++r) {
+      descriptor.cell_to_node_vector[jq][r] = gmsh_data.__quadrangles[j][r];
+    }
+    descriptor.cell_type_vector[jq]   = CellType::Quadrangle;
+    descriptor.cell_number_vector[jq] = gmsh_data.__quadrangles_number[j];
+  }
+
+  std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
+  for (unsigned int r = 0; r < gmsh_data.__triangles_ref.size(); ++r) {
+    const unsigned int elem_number = gmsh_data.__triangles_ref[r];
+    const unsigned int& ref        = gmsh_data.__triangles_ref[r];
+    ref_cells_map[ref].push_back(elem_number);
+  }
+
+  for (unsigned int j = 0; j < gmsh_data.__quadrangles_ref.size(); ++j) {
+    const size_t elem_number = nb_triangles + j;
+    const unsigned int& ref  = gmsh_data.__quadrangles_ref[j];
+    ref_cells_map[ref].push_back(elem_number);
+  }
+
+  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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
+    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
+  }
+
+  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
+
+  using Face                                                                 = ConnectivityFace<2>;
+  const auto& node_number_vector                                             = descriptor.node_number_vector;
+  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;
+  }();
+
+  std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
+  for (unsigned int e = 0; e < gmsh_data.__edges.size(); ++e) {
+    const unsigned int edge_id = [&] {
+      auto i = face_to_id_map.find(Face({gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]}, node_number_vector));
+      if (i == face_to_id_map.end()) {
+        std::stringstream error_msg;
+        error_msg << "face " << gmsh_data.__edges[e][0] << " not found";
+        throw NormalError(error_msg.str());
+      }
+      return i->second;
+    }();
+    const unsigned int& ref = gmsh_data.__edges_ref[e];
+    ref_faces_map[ref].push_back(edge_id);
+
+    if (descriptor.face_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
+      if (auto i_face = face_number_id_map.find(gmsh_data.__edges_number[e]); i_face != face_number_id_map.end()) {
+        const int other_edge_id = i_face->second;
+        std::swap(descriptor.face_number_vector[edge_id], descriptor.face_number_vector[other_edge_id]);
+
+        face_number_id_map.erase(descriptor.face_number_vector[edge_id]);
+        face_number_id_map.erase(descriptor.face_number_vector[other_edge_id]);
+
+        face_number_id_map[descriptor.face_number_vector[edge_id]]       = edge_id;
+        face_number_id_map[descriptor.face_number_vector[other_edge_id]] = other_edge_id;
+      } else {
+        face_number_id_map.erase(descriptor.face_number_vector[edge_id]);
+        descriptor.face_number_vector[edge_id]                     = gmsh_data.__edges_number[e];
+        face_number_id_map[descriptor.face_number_vector[edge_id]] = edge_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_face_list.first);
+    descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
+  }
+
+  std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+  for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
+    const unsigned int point_number = gmsh_data.__points[r];
+    const unsigned int& ref         = gmsh_data.__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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.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.node_owner_vector.resize(descriptor.node_number_vector.size());
+  std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+
+  m_connectivity = Connectivity2D::build(descriptor);
+}
+
+template <>
+GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells)
+{
+  ConnectivityDescriptor descriptor;
+
+  descriptor.node_number_vector = gmsh_data.__verticesNumbers;
+  descriptor.cell_type_vector.resize(nb_cells);
+  descriptor.cell_number_vector.resize(nb_cells);
+  descriptor.cell_to_node_vector.resize(nb_cells);
+
+  const size_t nb_tetrahedra = gmsh_data.__tetrahedra.size();
+  for (size_t j = 0; j < nb_tetrahedra; ++j) {
+    descriptor.cell_to_node_vector[j].resize(4);
+    for (int r = 0; r < 4; ++r) {
+      descriptor.cell_to_node_vector[j][r] = gmsh_data.__tetrahedra[j][r];
+    }
+    descriptor.cell_type_vector[j]   = CellType::Tetrahedron;
+    descriptor.cell_number_vector[j] = gmsh_data.__tetrahedra_number[j];
+  }
+  const size_t nb_hexahedra = gmsh_data.__hexahedra.size();
+  for (size_t j = 0; j < nb_hexahedra; ++j) {
+    const size_t jh = nb_tetrahedra + j;
+    descriptor.cell_to_node_vector[jh].resize(8);
+    for (int r = 0; r < 8; ++r) {
+      descriptor.cell_to_node_vector[jh][r] = gmsh_data.__hexahedra[j][r];
+    }
+    descriptor.cell_type_vector[jh]   = CellType::Hexahedron;
+    descriptor.cell_number_vector[jh] = gmsh_data.__hexahedra_number[j];
+  }
+
+  std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
+  for (unsigned int r = 0; r < gmsh_data.__tetrahedra_ref.size(); ++r) {
+    const unsigned int elem_number = gmsh_data.__tetrahedra_ref[r];
+    const unsigned int& ref        = gmsh_data.__tetrahedra_ref[r];
+    ref_cells_map[ref].push_back(elem_number);
+  }
+
+  for (unsigned int j = 0; j < gmsh_data.__hexahedra_ref.size(); ++j) {
+    const size_t elem_number = nb_tetrahedra + j;
+    const unsigned int& ref  = gmsh_data.__hexahedra_ref[j];
+    ref_cells_map[ref].push_back(elem_number);
+  }
+
+  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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
+    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
+  }
+
+  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
+
+  const auto& node_number_vector = descriptor.node_number_vector;
+
+  {
+    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;
+    }();
+
+    std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
+    for (unsigned int f = 0; f < gmsh_data.__triangles.size(); ++f) {
+      const unsigned int face_id = [&] {
+        auto i = face_to_id_map.find(
+          Face({gmsh_data.__triangles[f][0], gmsh_data.__triangles[f][1], gmsh_data.__triangles[f][2]},
+               node_number_vector));
+        if (i == face_to_id_map.end()) {
+          throw NormalError("face not found");
+        }
+        return i->second;
+      }();
+
+      const unsigned int& ref = gmsh_data.__triangles_ref[f];
+      ref_faces_map[ref].push_back(face_id);
+
+      if (descriptor.face_number_vector[face_id] != gmsh_data.__quadrangles_number[f]) {
+        if (auto i_face = face_number_id_map.find(gmsh_data.__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]                     = gmsh_data.__quadrangles_number[f];
+          face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+        }
+      }
+    }
+
+    for (unsigned int f = 0; f < gmsh_data.__quadrangles.size(); ++f) {
+      const unsigned int face_id = [&] {
+        auto i = face_to_id_map.find(Face({gmsh_data.__quadrangles[f][0], gmsh_data.__quadrangles[f][1],
+                                           gmsh_data.__quadrangles[f][2], gmsh_data.__quadrangles[f][3]},
+                                          node_number_vector));
+        if (i == face_to_id_map.end()) {
+          throw NormalError("face not found");
+        }
+        return i->second;
+      }();
+
+      const unsigned int& ref = gmsh_data.__quadrangles_ref[f];
+      ref_faces_map[ref].push_back(face_id);
+
+      if (descriptor.face_number_vector[face_id] != gmsh_data.__quadrangles_number[f]) {
+        if (auto i_face = face_number_id_map.find(gmsh_data.__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]                     = gmsh_data.__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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_face_list.first);
+      descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
+    }
+  }
+
+  ConnectivityBuilderBase::_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 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 < gmsh_data.__edges.size(); ++e) {
+      const unsigned int edge_id = [&] {
+        auto i = edge_to_id_map.find(Edge({gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]}, node_number_vector));
+        if (i == edge_to_id_map.end()) {
+          std::stringstream error_msg;
+          error_msg << "edge " << gmsh_data.__edges[e][0] << " not found";
+          throw NormalError(error_msg.str());
+        }
+        return i->second;
+      }();
+      const unsigned int& ref = gmsh_data.__edges_ref[e];
+      ref_edges_map[ref].push_back(edge_id);
+
+      if (descriptor.edge_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
+        if (auto i_edge = edge_number_id_map.find(gmsh_data.__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]                     = gmsh_data.__edges_number[e];
+          edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_edge_list.first);
+      descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list});
+    }
+  }
+
+  std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+  for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
+    const unsigned int point_number = gmsh_data.__points[r];
+    const unsigned int& ref         = gmsh_data.__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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.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());
+
+  m_connectivity = Connectivity3D::build(descriptor);
+}
+
 GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
 {
   if (parallel::rank() == 0) {
@@ -106,7 +550,6 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
         throw NormalError("Cannot read Gmsh format '" + std::to_string(fileVersion) + "'");
       }
       int fileType = this->_getInteger();
-      __binary     = (fileType == 1);
       if ((fileType < 0) or (fileType > 1)) {
         throw NormalError("Cannot read Gmsh file type '" + std::to_string(fileType) + "'");
       }
@@ -116,10 +559,6 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
         throw NormalError("Data size not supported '" + std::to_string(dataSize) + "'");
       }
 
-      if (__binary) {
-        throw NotImplementedError("cannot read binary files");
-      }
-
       kw = this->__nextKeyword();
       if (kw.second != ENDMESHFORMAT) {
         throw NormalError("reading file '" + m_filename + "': expecting $EndMeshFormat, '" + kw.first + "' was found");
@@ -189,71 +628,18 @@ GmshReader::__readVertices()
     throw NormalError("reading file '" + this->m_filename + "': number of vertices is negative");
   }
 
-  __verticesNumbers.resize(numberOfVerices);
-  __vertices = Array<R3>(numberOfVerices);
+  m_mesh_data.__verticesNumbers.resize(numberOfVerices);
+  m_mesh_data.__vertices = Array<R3>(numberOfVerices);
 
-  if (not __binary) {
-    for (int i = 0; i < numberOfVerices; ++i) {
-      __verticesNumbers[i] = this->_getInteger();
-      const double x       = this->_getReal();
-      const double y       = this->_getReal();
-      const double z       = this->_getReal();
-      __vertices[i]        = TinyVector<3, double>(x, y, z);
-    }
-  } else {
-    throw NotImplementedError("cannot read binary format");
+  for (int i = 0; i < numberOfVerices; ++i) {
+    m_mesh_data.__verticesNumbers[i] = this->_getInteger();
+    const double x                   = this->_getReal();
+    const double y                   = this->_getReal();
+    const double z                   = this->_getReal();
+    m_mesh_data.__vertices[i]        = TinyVector<3, double>(x, y, z);
   }
 }
 
-// void
-// GmshReader::__readElements1()
-// {
-//   const int numberOfElements = this->_getInteger();
-//   std::cout << "- Number Of Elements: " << numberOfElements << '\n';
-//   if (numberOfElements<0) {
-//     throw ErrorHandler(__FILE__,__LINE__,
-// 		       "reading file '"+m_filename
-// 		       +"': number of elements is negative",
-// 		       ErrorHandler::normal);
-//   }
-
-//   __elementType.resize(numberOfElements);
-//   __references.resize(numberOfElements);
-//   __elementVertices.resize(numberOfElements);
-
-//   for (int i=0; i<numberOfElements; ++i) {
-//     this->_getInteger(); // drops element number
-//     const int elementType = this->_getInteger()-1;
-
-//     if ((elementType < 0) or (elementType > 14)) {
-//       throw ErrorHandler(__FILE__,__LINE__,
-// 			 "reading file '"+m_filename
-// 			 +"': unknown element type
-// '"+std::to_string(elementType)+"'", 			 ErrorHandler::normal);
-//     }
-//     __elementType[i] = elementType;
-//     __references[i] = this->_getInteger(); // physical reference
-//     this->_getInteger(); // drops entity reference
-
-//     const int numberOfVertices = this->_getInteger();
-//     if (numberOfVertices != __numberOfPrimitiveNodes[elementType]) {
-//       std::stringstream errorMsg;
-//       errorMsg << "reading file '" <<m_filename
-// 	       << "':number of vertices (" << numberOfVertices
-// 	       << ") does not match expected ("
-// 	       << __numberOfPrimitiveNodes[elementType] << ")" << std::ends;
-
-//       throw ErrorHandler(__FILE__,__LINE__,
-// 			 errorMsg.str(),
-// 			 ErrorHandler::normal);
-//     }
-//     __elementVertices[i].resize(numberOfVertices);
-//     for (int j=0; j<numberOfVertices; ++j) {
-//       __elementVertices[i][j] = this->_getInteger();
-//     }
-//   }
-// }
-
 void
 GmshReader::__readElements2_2()
 {
@@ -263,96 +649,31 @@ GmshReader::__readElements2_2()
     throw NormalError("reading file '" + m_filename + "': number of elements is negative");
   }
 
-  __elementNumber.resize(numberOfElements);
-  __elementType.resize(numberOfElements);
-  __references.resize(numberOfElements);
-  __elementVertices.resize(numberOfElements);
-
-  if (not __binary) {
-    for (int i = 0; i < numberOfElements; ++i) {
-      __elementNumber[i]    = this->_getInteger();
-      const int elementType = this->_getInteger() - 1;
+  m_mesh_data.__elementNumber.resize(numberOfElements);
+  m_mesh_data.__elementType.resize(numberOfElements);
+  m_mesh_data.__references.resize(numberOfElements);
+  m_mesh_data.__elementVertices.resize(numberOfElements);
 
-      if ((elementType < 0) or (elementType > 14)) {
-        throw NormalError("reading file '" + m_filename + "': unknown element type '" + std::to_string(elementType) +
-                          "'");
-      }
-      __elementType[i]       = elementType;
-      const int numberOfTags = this->_getInteger();
-      __references[i]        = this->_getInteger();   // physical reference
-      for (int tag = 1; tag < numberOfTags; ++tag) {
-        this->_getInteger();   // drops remaining tags
-      }
+  for (int i = 0; i < numberOfElements; ++i) {
+    m_mesh_data.__elementNumber[i] = this->_getInteger();
+    const int elementType          = this->_getInteger() - 1;
 
-      const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
-      __elementVertices[i].resize(numberOfVertices);
-      for (int j = 0; j < numberOfVertices; ++j) {
-        __elementVertices[i][j] = this->_getInteger();
-      }
+    if ((elementType < 0) or (elementType > 14)) {
+      throw NormalError("reading file '" + m_filename + "': unknown element type '" + std::to_string(elementType) +
+                        "'");
+    }
+    m_mesh_data.__elementType[i] = elementType;
+    const int numberOfTags       = this->_getInteger();
+    m_mesh_data.__references[i]  = this->_getInteger();   // physical reference
+    for (int tag = 1; tag < numberOfTags; ++tag) {
+      this->_getInteger();   // drops remaining tags
     }
-  } else {
-    // fseek(__ifh,1L,SEEK_CUR);
-    // int i=0;
-    // for (;i<numberOfElements;) {
-    //   int elementType = 0;
-    //   fread(reinterpret_cast<char*>(&elementType),sizeof(int),1,__ifh);
-    //   if (__convertEndian) {
-    // 	invertEndianess(elementType);
-    //   }
-    //   --elementType;
-
-    //   if ((elementType < 0) or (elementType > 14)) {
-    // 	 throw ErrorHandler(__FILE__,__LINE__,
-    // 	 		   "reading file '"+m_filename
-    // 	 		   +"': unknown element type
-    // '"+std::to_string(elementType)+"'",
-    // ErrorHandler::normal);
-    //   }
-
-    //   int elementTypeNumber = 0;
-    //   fread(reinterpret_cast<char*>(&elementTypeNumber),sizeof(int),1,__ifh);
-    //   if (__convertEndian) {
-    // 	invertEndianess(elementTypeNumber);
-    //   }
-
-    //   int numberOfTags = 0;
-    //   fread(reinterpret_cast<char*>(&numberOfTags),sizeof(int),1,__ifh);
-    //   if (__convertEndian) {
-    // 	invertEndianess(numberOfTags);
-    //   }
-
-    //   for (int k=0; k<elementTypeNumber; ++k) {
-    // 	__elementType[i] = elementType;
-
-    // 	int elementNumber = 0;
-    // 	fread(reinterpret_cast<char*>(&elementNumber),sizeof(int),1,__ifh);
-
-    // 	int reference = 0;
-    // 	fread(reinterpret_cast<char*>(&reference),sizeof(int),1,__ifh);
-
-    // 	__references[i] = reference; // physical reference
-    // 	for (int tag=1; tag<numberOfTags; ++tag) {
-    // 	  int t;
-    // 	  fread(reinterpret_cast<char*>(&t),sizeof(int),1,__ifh);
-    // 	}
-
-    // 	const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
-    // 	__elementVertices[i].resize(numberOfVertices);
-    // 	fread(reinterpret_cast<char*>(&(__elementVertices[i][0])),sizeof(int),numberOfVertices,__ifh);
-    // 	++i;
-    //   }
-    // }
 
-    // if (__convertEndian) {
-    //   for (size_t i=0; i<__references.size(); ++i) {
-    // 	invertEndianess(__references[i]);
-    //   }
-    //   for (size_t i=0; i<__elementVertices.size(); ++i) {
-    // 	for (size_t j=0; j<__elementVertices[i].size(); ++i) {
-    // 	  invertEndianess(__elementVertices[i][j]);
-    // 	}
-    //   }
-    // }
+    const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
+    m_mesh_data.__elementVertices[i].resize(numberOfVertices);
+    for (int j = 0; j < numberOfVertices; ++j) {
+      m_mesh_data.__elementVertices[i][j] = this->_getInteger();
+    }
   }
 }
 
@@ -369,27 +690,277 @@ GmshReader::__readPhysicalNames2_2()
 
     PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
 
-    if (auto i_searched_physical_ref_id = m_physical_ref_map.find(physical_number);
-        i_searched_physical_ref_id != m_physical_ref_map.end()) {
+    if (auto i_searched_physical_ref_id = m_mesh_data.m_physical_ref_map.find(physical_number);
+        i_searched_physical_ref_id != m_mesh_data.m_physical_ref_map.end()) {
       std::stringstream os;
       os << "Physical reference id '" << physical_ref_id << "' already defined as '"
          << i_searched_physical_ref_id->second << "'!";
       throw NormalError(os.str());
     }
-    m_physical_ref_map[physical_number] = physical_ref_id;
+    m_mesh_data.m_physical_ref_map[physical_number] = physical_ref_id;
   }
 }
 
+// std::shared_ptr<IConnectivity>
+// GmshReader::_buildConnectivity3D(const size_t nb_cells)
+// {
+//   ConnectivityDescriptor descriptor;
+
+//   descriptor.node_number_vector = m_mesh_data.__verticesNumbers;
+//   descriptor.cell_type_vector.resize(nb_cells);
+//   descriptor.cell_number_vector.resize(nb_cells);
+//   descriptor.cell_to_node_vector.resize(nb_cells);
+
+//   const size_t nb_tetrahedra = m_mesh_data.__tetrahedra.size();
+//   for (size_t j = 0; j < nb_tetrahedra; ++j) {
+//     descriptor.cell_to_node_vector[j].resize(4);
+//     for (int r = 0; r < 4; ++r) {
+//       descriptor.cell_to_node_vector[j][r] = m_mesh_data.__tetrahedra[j][r];
+//     }
+//     descriptor.cell_type_vector[j]   = CellType::Tetrahedron;
+//     descriptor.cell_number_vector[j] = m_mesh_data.__tetrahedra_number[j];
+//   }
+//   const size_t nb_hexahedra = m_mesh_data.__hexahedra.size();
+//   for (size_t j = 0; j < nb_hexahedra; ++j) {
+//     const size_t jh = nb_tetrahedra + j;
+//     descriptor.cell_to_node_vector[jh].resize(8);
+//     for (int r = 0; r < 8; ++r) {
+//       descriptor.cell_to_node_vector[jh][r] = m_mesh_data.__hexahedra[j][r];
+//     }
+//     descriptor.cell_type_vector[jh]   = CellType::Hexahedron;
+//     descriptor.cell_number_vector[jh] = m_mesh_data.__hexahedra_number[j];
+//   }
+
+//   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
+//   for (unsigned int r = 0; r < m_mesh_data.__tetrahedra_ref.size(); ++r) {
+//     const unsigned int elem_number = m_mesh_data.__tetrahedra_ref[r];
+//     const unsigned int& ref        = m_mesh_data.__tetrahedra_ref[r];
+//     ref_cells_map[ref].push_back(elem_number);
+//   }
+
+//   for (unsigned int j = 0; j < m_mesh_data.__hexahedra_ref.size(); ++j) {
+//     const size_t elem_number = nb_tetrahedra + j;
+//     const unsigned int& ref  = m_mesh_data.__hexahedra_ref[j];
+//     ref_cells_map[ref].push_back(elem_number);
+//   }
+
+//   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_mesh_data.m_physical_ref_map.at(ref_cell_list.first);
+//     descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
+//   }
+
+//   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
+
+//   const auto& node_number_vector = descriptor.node_number_vector;
+
+//   {
+//     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;
+//     }();
+
+//     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
+//     for (unsigned int f = 0; f < m_mesh_data.__triangles.size(); ++f) {
+//       const unsigned int face_id = [&] {
+//         auto i = face_to_id_map.find(
+//           Face({m_mesh_data.__triangles[f][0], m_mesh_data.__triangles[f][1], m_mesh_data.__triangles[f][2]},
+//                node_number_vector));
+//         if (i == face_to_id_map.end()) {
+//           throw NormalError("face not found");
+//         }
+//         return i->second;
+//       }();
+
+//       const unsigned int& ref = m_mesh_data.__triangles_ref[f];
+//       ref_faces_map[ref].push_back(face_id);
+
+//       if (descriptor.face_number_vector[face_id] != m_mesh_data.__quadrangles_number[f]) {
+//         if (auto i_face = face_number_id_map.find(m_mesh_data.__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]                     = m_mesh_data.__quadrangles_number[f];
+//           face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+//         }
+//       }
+//     }
+
+//     for (unsigned int f = 0; f < m_mesh_data.__quadrangles.size(); ++f) {
+//       const unsigned int face_id = [&] {
+//         auto i = face_to_id_map.find(Face({m_mesh_data.__quadrangles[f][0], m_mesh_data.__quadrangles[f][1],
+//                                            m_mesh_data.__quadrangles[f][2], m_mesh_data.__quadrangles[f][3]},
+//                                           node_number_vector));
+//         if (i == face_to_id_map.end()) {
+//           throw NormalError("face not found");
+//         }
+//         return i->second;
+//       }();
+
+//       const unsigned int& ref = m_mesh_data.__quadrangles_ref[f];
+//       ref_faces_map[ref].push_back(face_id);
+
+//       if (descriptor.face_number_vector[face_id] != m_mesh_data.__quadrangles_number[f]) {
+//         if (auto i_face = face_number_id_map.find(m_mesh_data.__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]                     = m_mesh_data.__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_mesh_data.m_physical_ref_map.at(ref_face_list.first);
+//       descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
+//     }
+//   }
+
+//   ConnectivityBuilderBase::_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 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 < m_mesh_data.__edges.size(); ++e) {
+//       const unsigned int edge_id = [&] {
+//         auto i = edge_to_id_map.find(Edge({m_mesh_data.__edges[e][0], m_mesh_data.__edges[e][1]},
+//         node_number_vector)); if (i == edge_to_id_map.end()) {
+//           std::stringstream error_msg;
+//           error_msg << "edge " << m_mesh_data.__edges[e][0] << " not found";
+//           throw NormalError(error_msg.str());
+//         }
+//         return i->second;
+//       }();
+//       const unsigned int& ref = m_mesh_data.__edges_ref[e];
+//       ref_edges_map[ref].push_back(edge_id);
+
+//       if (descriptor.edge_number_vector[edge_id] != m_mesh_data.__edges_number[e]) {
+//         if (auto i_edge = edge_number_id_map.find(m_mesh_data.__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]                     = m_mesh_data.__edges_number[e];
+//           edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_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_mesh_data.m_physical_ref_map.at(ref_edge_list.first);
+//       descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list});
+//     }
+//   }
+
+//   std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+//   for (unsigned int r = 0; r < m_mesh_data.__points.size(); ++r) {
+//     const unsigned int point_number = m_mesh_data.__points[r];
+//     const unsigned int& ref         = m_mesh_data.__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_mesh_data.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());
+
+//   return Connectivity3D::build(descriptor);
+// }
+
 void
 GmshReader::__proceedData()
 {
   TinyVector<15, int> elementNumber = zero;
   std::vector<std::set<size_t>> elementReferences(15);
 
-  for (size_t i = 0; i < __elementType.size(); ++i) {
-    const int elementType = __elementType[i];
+  for (size_t i = 0; i < m_mesh_data.__elementType.size(); ++i) {
+    const int elementType = m_mesh_data.__elementType[i];
     elementNumber[elementType]++;
-    elementReferences[elementType].insert(__references[i]);
+    elementReferences[elementType].insert(m_mesh_data.__references[i]);
   }
 
   for (size_t i = 0; i < elementNumber.dimension(); ++i) {
@@ -409,39 +980,39 @@ GmshReader::__proceedData()
         switch (i) {
         // Supported entities
         case 0: {   // edges
-          __edges = Array<Edge>(elementNumber[i]);
-          __edges_ref.resize(elementNumber[i]);
-          __edges_number.resize(elementNumber[i]);
+          m_mesh_data.__edges = Array<GmshData::Edge>(elementNumber[i]);
+          m_mesh_data.__edges_ref.resize(elementNumber[i]);
+          m_mesh_data.__edges_number.resize(elementNumber[i]);
           break;
         }
         case 1: {   // triangles
-          __triangles = Array<Triangle>(elementNumber[i]);
-          __triangles_ref.resize(elementNumber[i]);
-          __triangles_number.resize(elementNumber[i]);
+          m_mesh_data.__triangles = Array<GmshData::Triangle>(elementNumber[i]);
+          m_mesh_data.__triangles_ref.resize(elementNumber[i]);
+          m_mesh_data.__triangles_number.resize(elementNumber[i]);
           break;
         }
         case 2: {   // quadrangles
-          __quadrangles = Array<Quadrangle>(elementNumber[i]);
-          __quadrangles_ref.resize(elementNumber[i]);
-          __quadrangles_number.resize(elementNumber[i]);
+          m_mesh_data.__quadrangles = Array<GmshData::Quadrangle>(elementNumber[i]);
+          m_mesh_data.__quadrangles_ref.resize(elementNumber[i]);
+          m_mesh_data.__quadrangles_number.resize(elementNumber[i]);
           break;
         }
         case 3: {   // tetrahedra
-          __tetrahedra = Array<Tetrahedron>(elementNumber[i]);
-          __tetrahedra_ref.resize(elementNumber[i]);
-          __tetrahedra_number.resize(elementNumber[i]);
+          m_mesh_data.__tetrahedra = Array<GmshData::Tetrahedron>(elementNumber[i]);
+          m_mesh_data.__tetrahedra_ref.resize(elementNumber[i]);
+          m_mesh_data.__tetrahedra_number.resize(elementNumber[i]);
           break;
         }
         case 4: {   // hexahedra
-          __hexahedra = Array<Hexahedron>(elementNumber[i]);
-          __hexahedra_ref.resize(elementNumber[i]);
-          __hexahedra_number.resize(elementNumber[i]);
+          m_mesh_data.__hexahedra = Array<GmshData::Hexahedron>(elementNumber[i]);
+          m_mesh_data.__hexahedra_ref.resize(elementNumber[i]);
+          m_mesh_data.__hexahedra_number.resize(elementNumber[i]);
           break;
         }
         case 14: {   // point
-          __points = Array<Point>(elementNumber[i]);
-          __points_ref.resize(elementNumber[i]);
-          __points_number.resize(elementNumber[i]);
+          m_mesh_data.__points = Array<GmshData::Point>(elementNumber[i]);
+          m_mesh_data.__points_ref.resize(elementNumber[i]);
+          m_mesh_data.__points_number.resize(elementNumber[i]);
           break;
         }
           // Unsupported entities
@@ -466,8 +1037,8 @@ GmshReader::__proceedData()
   std::cout << "- Building correspondance table\n";
   int minNumber = std::numeric_limits<int>::max();
   int maxNumber = std::numeric_limits<int>::min();
-  for (size_t i = 0; i < __verticesNumbers.size(); ++i) {
-    const int& vertexNumber = __verticesNumbers[i];
+  for (size_t i = 0; i < m_mesh_data.__verticesNumbers.size(); ++i) {
+    const int& vertexNumber = m_mesh_data.__verticesNumbers[i];
     minNumber               = std::min(minNumber, vertexNumber);
     maxNumber               = std::max(maxNumber, vertexNumber);
   }
@@ -477,112 +1048,112 @@ GmshReader::__proceedData()
   }
 
   // A value of -1 means that the vertex is unknown
-  __verticesCorrepondance.resize(maxNumber + 1, -1);
+  m_mesh_data.__verticesCorrepondance.resize(maxNumber + 1, -1);
 
-  for (size_t i = 0; i < __verticesNumbers.size(); ++i) {
-    __verticesCorrepondance[__verticesNumbers[i]] = i;
+  for (size_t i = 0; i < m_mesh_data.__verticesNumbers.size(); ++i) {
+    m_mesh_data.__verticesCorrepondance[m_mesh_data.__verticesNumbers[i]] = i;
   }
 
   // reset element number to count them while filling structures
   elementNumber = zero;
 
   // Element structures filling
-  for (size_t i = 0; i < __elementType.size(); ++i) {
-    std::vector<int>& elementVertices = __elementVertices[i];
-    switch (__elementType[i]) {
+  for (size_t i = 0; i < m_mesh_data.__elementType.size(); ++i) {
+    std::vector<int>& elementVertices = m_mesh_data.__elementVertices[i];
+    switch (m_mesh_data.__elementType[i]) {
     // Supported entities
     case 0: {   // edge
       int& edgeNumber = elementNumber[0];
-      const int a     = __verticesCorrepondance[elementVertices[0]];
-      const int b     = __verticesCorrepondance[elementVertices[1]];
+      const int a     = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+      const int b     = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
       if ((a < 0) or (b < 0)) {
         throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
                           " [bad vertices definition]");
       }
-      __edges[edgeNumber]        = Edge(a, b);
-      __edges_ref[edgeNumber]    = __references[i];
-      __edges_number[edgeNumber] = __elementNumber[i];
+      m_mesh_data.__edges[edgeNumber]        = GmshData::Edge(a, b);
+      m_mesh_data.__edges_ref[edgeNumber]    = m_mesh_data.__references[i];
+      m_mesh_data.__edges_number[edgeNumber] = m_mesh_data.__elementNumber[i];
       edgeNumber++;   // one more edge
       break;
     }
     case 1: {   // triangles
       int& triangleNumber = elementNumber[1];
 
-      const int a = __verticesCorrepondance[elementVertices[0]];
-      const int b = __verticesCorrepondance[elementVertices[1]];
-      const int c = __verticesCorrepondance[elementVertices[2]];
+      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
       if ((a < 0) or (b < 0) or (c < 0)) {
         throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
                           " [bad vertices definition]");
       }
-      __triangles[triangleNumber]        = Triangle(a, b, c);
-      __triangles_ref[triangleNumber]    = __references[i];
-      __triangles_number[triangleNumber] = __elementNumber[i];
+      m_mesh_data.__triangles[triangleNumber]        = GmshData::Triangle(a, b, c);
+      m_mesh_data.__triangles_ref[triangleNumber]    = m_mesh_data.__references[i];
+      m_mesh_data.__triangles_number[triangleNumber] = m_mesh_data.__elementNumber[i];
       triangleNumber++;   // one more triangle
       break;
     }
     case 2: {   // quadrangle
       int& quadrilateralNumber = elementNumber[2];
 
-      const int a = __verticesCorrepondance[elementVertices[0]];
-      const int b = __verticesCorrepondance[elementVertices[1]];
-      const int c = __verticesCorrepondance[elementVertices[2]];
-      const int d = __verticesCorrepondance[elementVertices[3]];
+      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
+      const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
         throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
                           " [bad vertices definition]");
       }
-      __quadrangles[quadrilateralNumber]        = Quadrangle(a, b, c, d);
-      __quadrangles_ref[quadrilateralNumber]    = __references[i];
-      __quadrangles_number[quadrilateralNumber] = __elementNumber[i];
+      m_mesh_data.__quadrangles[quadrilateralNumber]        = GmshData::Quadrangle(a, b, c, d);
+      m_mesh_data.__quadrangles_ref[quadrilateralNumber]    = m_mesh_data.__references[i];
+      m_mesh_data.__quadrangles_number[quadrilateralNumber] = m_mesh_data.__elementNumber[i];
       quadrilateralNumber++;   // one more quadrangle
       break;
     }
     case 3: {   // tetrahedra
       int& tetrahedronNumber = elementNumber[3];
 
-      const int a = __verticesCorrepondance[elementVertices[0]];
-      const int b = __verticesCorrepondance[elementVertices[1]];
-      const int c = __verticesCorrepondance[elementVertices[2]];
-      const int d = __verticesCorrepondance[elementVertices[3]];
+      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
+      const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
         throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
                           " [bad vertices definition]");
       }
-      __tetrahedra[tetrahedronNumber]        = Tetrahedron(a, b, c, d);
-      __tetrahedra_ref[tetrahedronNumber]    = __references[i];
-      __tetrahedra_number[tetrahedronNumber] = __elementNumber[i];
+      m_mesh_data.__tetrahedra[tetrahedronNumber]        = GmshData::Tetrahedron(a, b, c, d);
+      m_mesh_data.__tetrahedra_ref[tetrahedronNumber]    = m_mesh_data.__references[i];
+      m_mesh_data.__tetrahedra_number[tetrahedronNumber] = m_mesh_data.__elementNumber[i];
       tetrahedronNumber++;   // one more tetrahedron
       break;
     }
     case 4: {   // hexaredron
       int& hexahedronNumber = elementNumber[4];
 
-      const int a = __verticesCorrepondance[elementVertices[0]];
-      const int b = __verticesCorrepondance[elementVertices[1]];
-      const int c = __verticesCorrepondance[elementVertices[2]];
-      const int d = __verticesCorrepondance[elementVertices[3]];
-      const int e = __verticesCorrepondance[elementVertices[4]];
-      const int f = __verticesCorrepondance[elementVertices[5]];
-      const int g = __verticesCorrepondance[elementVertices[6]];
-      const int h = __verticesCorrepondance[elementVertices[7]];
+      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
+      const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
+      const int e = m_mesh_data.__verticesCorrepondance[elementVertices[4]];
+      const int f = m_mesh_data.__verticesCorrepondance[elementVertices[5]];
+      const int g = m_mesh_data.__verticesCorrepondance[elementVertices[6]];
+      const int h = m_mesh_data.__verticesCorrepondance[elementVertices[7]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0) or (f < 0) or (g < 0) or (h < 0)) {
         throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
                           " [bad vertices definition]");
       }
-      __hexahedra[hexahedronNumber]        = Hexahedron(a, b, c, d, e, f, g, h);
-      __hexahedra_ref[hexahedronNumber]    = __references[i];
-      __hexahedra_number[hexahedronNumber] = __elementNumber[i];
+      m_mesh_data.__hexahedra[hexahedronNumber]        = GmshData::Hexahedron(a, b, c, d, e, f, g, h);
+      m_mesh_data.__hexahedra_ref[hexahedronNumber]    = m_mesh_data.__references[i];
+      m_mesh_data.__hexahedra_number[hexahedronNumber] = m_mesh_data.__elementNumber[i];
       hexahedronNumber++;   // one more hexahedron
       break;
     }
       // Unsupported entities
     case 14: {   // point
-      int& point_number             = elementNumber[14];
-      const int a                   = __verticesCorrepondance[elementVertices[0]];
-      __points[point_number]        = a;
-      __points_ref[point_number]    = __references[i];
-      __points_number[point_number] = __elementNumber[i];
+      int& point_number                         = elementNumber[14];
+      const int a                               = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+      m_mesh_data.__points[point_number]        = a;
+      m_mesh_data.__points_ref[point_number]    = m_mesh_data.__references[i];
+      m_mesh_data.__points_number[point_number] = m_mesh_data.__elementNumber[i];
       point_number++;
       break;
     }
@@ -629,463 +1200,57 @@ GmshReader::__proceedData()
   if ((dimension3_mask, elementNumber) > 0) {
     const size_t nb_cells = (dimension3_mask, elementNumber);
 
-    ConnectivityDescriptor descriptor;
-
-    descriptor.node_number_vector = __verticesNumbers;
-    descriptor.cell_type_vector.resize(nb_cells);
-    descriptor.cell_number_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_to_node_vector[j].resize(4);
-      for (int r = 0; r < 4; ++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];
-    }
-    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_to_node_vector[jh].resize(8);
-      for (int r = 0; r < 8; ++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];
-    }
-
-    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);
-    }
-
-    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);
-    }
-
-    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));
-    }
-
-    MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
-
-    const auto& node_number_vector = descriptor.node_number_vector;
-
-    {
-      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;
-      }();
-
-      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()) {
-            throw NormalError("face not found");
-          }
-          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;
-          }
-        }
-      }
-
-      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()) {
-            throw NormalError("face not found");
-          }
-          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;
-          }
-        }
-      }
-
-      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});
-      }
-    }
-    MeshBuilderBase::_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 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::stringstream error_msg;
-            error_msg << "edge " << __edges[e][0] << " not found";
-            throw NormalError(error_msg.str());
-          }
-          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;
-          }
-        }
-      }
-
-      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});
-      }
-    }
-
-    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_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());
+    GmshConnectivityBuilder<3> connectivity_builder(m_mesh_data, nb_cells);
 
-    std::shared_ptr p_connectivity = Connectivity3D::build(descriptor);
-    Connectivity3D& connectivity   = *p_connectivity;
+    std::shared_ptr p_connectivity =
+      std::dynamic_pointer_cast<const Connectivity3D>(connectivity_builder.connectivity());
+    const Connectivity3D& connectivity = *p_connectivity;
 
     using MeshType = Mesh<Connectivity3D>;
     using Rd       = TinyVector<3, double>;
 
     NodeValue<Rd> xr(connectivity);
-    for (NodeId i = 0; i < __vertices.size(); ++i) {
-      xr[i][0] = __vertices[i][0];
-      xr[i][1] = __vertices[i][1];
-      xr[i][2] = __vertices[i][2];
+    for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
+      xr[i][0] = m_mesh_data.__vertices[i][0];
+      xr[i][1] = m_mesh_data.__vertices[i][1];
+      xr[i][2] = m_mesh_data.__vertices[i][2];
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
 
   } else if ((dimension2_mask, elementNumber) > 0) {
     const size_t nb_cells = (dimension2_mask, elementNumber);
 
-    ConnectivityDescriptor descriptor;
-
-    descriptor.node_number_vector = __verticesNumbers;
-    descriptor.cell_type_vector.resize(nb_cells);
-    descriptor.cell_number_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_to_node_vector[j].resize(3);
-      for (int r = 0; r < 3; ++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];
-    }
-
-    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_to_node_vector[jq].resize(4);
-      for (int r = 0; r < 4; ++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];
-    }
-
-    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);
-    }
-
-    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);
-    }
-
-    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));
-    }
-
-    MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
-
-    using Face                                                                 = ConnectivityFace<2>;
-    const auto& node_number_vector                                             = descriptor.node_number_vector;
-    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;
-    }();
-
-    std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
-    for (unsigned int e = 0; e < __edges.size(); ++e) {
-      const unsigned int edge_id = [&] {
-        auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}, node_number_vector));
-        if (i == face_to_id_map.end()) {
-          std::stringstream error_msg;
-          error_msg << "face " << __edges[e][0] << " not found";
-          throw NormalError(error_msg.str());
-        }
-        return i->second;
-      }();
-      const unsigned int& ref = __edges_ref[e];
-      ref_faces_map[ref].push_back(edge_id);
-
-      if (descriptor.face_number_vector[edge_id] != __edges_number[e]) {
-        if (auto i_face = face_number_id_map.find(__edges_number[e]); i_face != face_number_id_map.end()) {
-          const int other_edge_id = i_face->second;
-          std::swap(descriptor.face_number_vector[edge_id], descriptor.face_number_vector[other_edge_id]);
-
-          face_number_id_map.erase(descriptor.face_number_vector[edge_id]);
-          face_number_id_map.erase(descriptor.face_number_vector[other_edge_id]);
-
-          face_number_id_map[descriptor.face_number_vector[edge_id]]       = edge_id;
-          face_number_id_map[descriptor.face_number_vector[other_edge_id]] = other_edge_id;
-        } else {
-          face_number_id_map.erase(descriptor.face_number_vector[edge_id]);
-          descriptor.face_number_vector[edge_id]                     = __edges_number[e];
-          face_number_id_map[descriptor.face_number_vector[edge_id]] = edge_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});
-    }
-
-    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_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.node_owner_vector.resize(descriptor.node_number_vector.size());
-    std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+    GmshConnectivityBuilder<2> connectivity_builder(m_mesh_data, nb_cells);
 
-    std::shared_ptr p_connectivity = Connectivity2D::build(descriptor);
-    Connectivity2D& connectivity   = *p_connectivity;
+    std::shared_ptr p_connectivity =
+      std::dynamic_pointer_cast<const Connectivity2D>(connectivity_builder.connectivity());
+    const Connectivity2D& connectivity = *p_connectivity;
 
     using MeshType = Mesh<Connectivity2D>;
     using Rd       = TinyVector<2, double>;
 
     NodeValue<Rd> xr(connectivity);
-    for (NodeId i = 0; i < __vertices.size(); ++i) {
-      xr[i][0] = __vertices[i][0];
-      xr[i][1] = __vertices[i][1];
+    for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
+      xr[i][0] = m_mesh_data.__vertices[i][0];
+      xr[i][1] = m_mesh_data.__vertices[i][1];
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
 
   } else if ((dimension1_mask, elementNumber) > 0) {
     const size_t nb_cells = (dimension1_mask, elementNumber);
 
-    ConnectivityDescriptor descriptor;
-
-    descriptor.node_number_vector = __verticesNumbers;
-    descriptor.cell_type_vector.resize(nb_cells);
-    descriptor.cell_number_vector.resize(nb_cells);
-    descriptor.cell_to_node_vector.resize(nb_cells);
-
-    for (size_t j = 0; j < nb_cells; ++j) {
-      descriptor.cell_to_node_vector[j].resize(2);
-      for (int r = 0; r < 2; ++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];
-    }
-
-    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_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.node_owner_vector.resize(descriptor.node_number_vector.size());
-    std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+    GmshConnectivityBuilder<1> connectivity_builder(m_mesh_data, nb_cells);
 
-    std::shared_ptr p_connectivity = Connectivity1D::build(descriptor);
-    Connectivity1D& connectivity   = *p_connectivity;
+    std::shared_ptr p_connectivity =
+      std::dynamic_pointer_cast<const Connectivity1D>(connectivity_builder.connectivity());
+    const Connectivity1D& connectivity = *p_connectivity;
 
     using MeshType = Mesh<Connectivity1D>;
     using Rd       = TinyVector<1, double>;
 
     NodeValue<Rd> xr(connectivity);
-    for (NodeId i = 0; i < __vertices.size(); ++i) {
-      xr[i][0] = __vertices[i][0];
+    for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
+      xr[i][0] = m_mesh_data.__vertices[i][0];
     }
 
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 4edd1fc7b081d67bf74fa132ed4407f43896346d..ae18d30ca3acd989f7a48d610afb458de8c4342d 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -69,73 +69,80 @@ class GmshReader : public MeshBuilderBase
     ~PhysicalRefId()                    = default;
   };
 
+  struct GmshData
+  {
+    /**
+     * Gmsh format provides a numbered, none ordrered and none dense
+     * vertices list, this stores the number of read vertices.
+     */
+    std::vector<int> __verticesNumbers;
+
+    Array<R3> __vertices;
+
+    using Point = unsigned int;
+    Array<Point> __points;
+    std::vector<int> __points_ref;
+    std::vector<int> __points_number;
+
+    using Edge = TinyVector<2, unsigned int>;
+    Array<Edge> __edges;
+    std::vector<int> __edges_ref;
+    std::vector<int> __edges_number;
+
+    using Triangle = TinyVector<3, unsigned int>;
+    Array<Triangle> __triangles;
+    std::vector<int> __triangles_ref;
+    std::vector<int> __triangles_number;
+
+    using Quadrangle = TinyVector<4, unsigned int>;
+    Array<Quadrangle> __quadrangles;
+    std::vector<int> __quadrangles_ref;
+    std::vector<int> __quadrangles_number;
+
+    using Tetrahedron = TinyVector<4, unsigned int>;
+    Array<Tetrahedron> __tetrahedra;
+    std::vector<int> __tetrahedra_ref;
+    std::vector<int> __tetrahedra_number;
+
+    using Hexahedron = TinyVector<8, unsigned int>;
+    Array<Hexahedron> __hexahedra;
+    std::vector<int> __hexahedra_ref;
+    std::vector<int> __hexahedra_number;
+
+    /**
+     * Gmsh format provides a numbered, none ordrered and none dense
+     * vertices list, this provides vertices renumbering correspondence
+     */
+    std::vector<int> __verticesCorrepondance;
+
+    /**
+     * elements types
+     */
+    std::vector<int> __elementNumber;
+
+    /**
+     * elements types
+     */
+    std::vector<short> __elementType;
+
+    /**
+     * References
+     */
+    std::vector<int> __references;
+
+    /**
+     * References
+     */
+    std::vector<std::vector<int> > __elementVertices;
+
+    std::map<unsigned int, PhysicalRefId> m_physical_ref_map;
+  };
+
  private:
   std::ifstream m_fin;
   const std::string m_filename;
 
-  /**
-   * Gmsh format provides a numbered, none ordrered and none dense
-   * vertices list, this stores the number of read vertices.
-   */
-  std::vector<int> __verticesNumbers;
-
-  Array<R3> __vertices;
-
-  using Point = unsigned int;
-  Array<Point> __points;
-  std::vector<int> __points_ref;
-  std::vector<int> __points_number;
-
-  using Edge = TinyVector<2, unsigned int>;
-  Array<Edge> __edges;
-  std::vector<int> __edges_ref;
-  std::vector<int> __edges_number;
-
-  using Triangle = TinyVector<3, unsigned int>;
-  Array<Triangle> __triangles;
-  std::vector<int> __triangles_ref;
-  std::vector<int> __triangles_number;
-
-  using Quadrangle = TinyVector<4, unsigned int>;
-  Array<Quadrangle> __quadrangles;
-  std::vector<int> __quadrangles_ref;
-  std::vector<int> __quadrangles_number;
-
-  using Tetrahedron = TinyVector<4, unsigned int>;
-  Array<Tetrahedron> __tetrahedra;
-  std::vector<int> __tetrahedra_ref;
-  std::vector<int> __tetrahedra_number;
-
-  using Hexahedron = TinyVector<8, unsigned int>;
-  Array<Hexahedron> __hexahedra;
-  std::vector<int> __hexahedra_ref;
-  std::vector<int> __hexahedra_number;
-
-  /**
-   * Gmsh format provides a numbered, none ordrered and none dense
-   * vertices list, this provides vertices renumbering correspondance
-   */
-  std::vector<int> __verticesCorrepondance;
-
-  /**
-   * elements types
-   */
-  std::vector<int> __elementNumber;
-
-  /**
-   * elements types
-   */
-  std::vector<short> __elementType;
-
-  /**
-   * References
-   */
-  std::vector<int> __references;
-
-  /**
-   * References
-   */
-  std::vector<std::vector<int> > __elementVertices;
+  GmshData m_mesh_data;
 
   /**
    * Stores the number of nodes associated to each primitive
@@ -155,11 +162,6 @@ class GmshReader : public MeshBuilderBase
    */
   std::map<int, std::string> __primitivesNames;
 
-  bool __binary;        /**< true if binary format */
-  bool __convertEndian; /**< true if needs to adapt endianess */
-
-  std::map<unsigned int, PhysicalRefId> m_physical_ref_map;
-
   int
   _getInteger()
   {
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index 7f7221b049872fac2d78b936093be706f03ea080..fc9020325fe64cb2d27dd3feab97ffcf778cdb3c 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -38,9 +38,15 @@ MeshBuilderBase::_dispatch()
   m_mesh = std::make_shared<MeshType>(dispatched_connectivity, dispatched_xr);
 }
 
+template void MeshBuilderBase::_dispatch<1>();
+template void MeshBuilderBase::_dispatch<2>();
+template void MeshBuilderBase::_dispatch<3>();
+
+/// ============
+
 template <size_t Dimension>
 void
-MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
+ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
 {
   static_assert((Dimension == 2) or (Dimension == 3), "Invalid dimension to compute cell-face connectivities");
   using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
@@ -247,7 +253,7 @@ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescripto
 
 template <size_t Dimension>
 void
-MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
+ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
 {
   static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
   using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
@@ -437,12 +443,8 @@ MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Connectivi
   }
 }
 
-template void MeshBuilderBase::_dispatch<1>();
-template void MeshBuilderBase::_dispatch<2>();
-template void MeshBuilderBase::_dispatch<3>();
-
-template void MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
-template void MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
+template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
+template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
 
-template void MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
+template void ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
   ConnectivityDescriptor& descriptor);
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
index 03eae11d52dda3e593e6e8dba7e589e3c5f40ab1..24bd055824223a604d1a913d1e85af0a757b1965 100644
--- a/src/mesh/MeshBuilderBase.hpp
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -2,6 +2,9 @@
 #define MESH_BUILDER_BASE_HPP
 
 #include <mesh/IMesh.hpp>
+
+#include <mesh/IConnectivity.hpp>
+
 #include <mesh/ItemId.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
@@ -11,13 +14,13 @@
 
 class ConnectivityDescriptor;
 
-class MeshBuilderBase
+class ConnectivityBuilderBase
 {
  protected:
   template <size_t Dimension>
   class ConnectivityFace;
 
-  std::shared_ptr<const IMesh> m_mesh;
+  std::shared_ptr<const IConnectivity> m_connectivity;
 
   template <size_t Dimension>
   static void _computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
@@ -25,22 +28,19 @@ class MeshBuilderBase
   template <size_t Dimension>
   static void _computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
 
-  template <int Dimension>
-  void _dispatch();
-
  public:
-  std::shared_ptr<const IMesh>
-  mesh() const
+  std::shared_ptr<const IConnectivity>
+  connectivity() const
   {
-    return m_mesh;
+    return m_connectivity;
   }
 
-  MeshBuilderBase()  = default;
-  ~MeshBuilderBase() = default;
+  ConnectivityBuilderBase()  = default;
+  ~ConnectivityBuilderBase() = default;
 };
 
 template <>
-class MeshBuilderBase::ConnectivityFace<2>
+class ConnectivityBuilderBase::ConnectivityFace<2>
 {
  public:
   friend struct Hash;
@@ -119,7 +119,7 @@ class MeshBuilderBase::ConnectivityFace<2>
 };
 
 template <>
-class MeshBuilderBase::ConnectivityFace<3>
+class ConnectivityBuilderBase::ConnectivityFace<3>
 {
  public:
   friend struct Hash;
@@ -225,4 +225,23 @@ class MeshBuilderBase::ConnectivityFace<3>
   ~ConnectivityFace() = default;
 };
 
+class MeshBuilderBase
+{
+ protected:
+  std::shared_ptr<const IMesh> m_mesh;
+
+  template <int Dimension>
+  void _dispatch();
+
+ public:
+  std::shared_ptr<const IMesh>
+  mesh() const
+  {
+    return m_mesh;
+  }
+
+  MeshBuilderBase()  = default;
+  ~MeshBuilderBase() = default;
+};
+
 #endif   // MESH_BUILDER_BASE_HPP