diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 0181e5b2cf87469c8481a8ec5b3b70414473a321..b9a65729637128081a7206ad6fa5d98ad95e9978 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/DiamondDualConnectivityBuilder.hpp>
+#include <mesh/DiamondDualMeshBuilder.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 97b166608fae5c7bc454dfe063bf60c4e6f18ba2..657b3b37500304446627038514ede886ad205da4 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -3,11 +3,14 @@
 add_library(
   PugsMesh
   CartesianMeshBuilder.cpp
-  DiamondDualConnectivityBuilder.cpp
   Connectivity.cpp
+  ConnectivityBuilderBase.cpp
   ConnectivityComputer.cpp
   ConnectivityDispatcher.cpp
+  DiamondDualConnectivityBuilder.cpp
+  DiamondDualMeshBuilder.cpp
   GmshReader.cpp
+  LogicalConnectivityBuilder.cpp
   MeshBuilderBase.cpp
   SynchronizerManager.cpp)
 
diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
index caff700fab42228f8aa2a1cd0b246e736881de0e..d6a4589cced3958ca7401a5a80d881d23fe87749 100644
--- a/src/mesh/CartesianMeshBuilder.cpp
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -1,8 +1,7 @@
 #include <mesh/CartesianMeshBuilder.hpp>
 
 #include <mesh/Connectivity.hpp>
-#include <mesh/ConnectivityDescriptor.hpp>
-#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/LogicalConnectivityBuilder.hpp>
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
@@ -159,639 +158,3 @@ template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<2>&,
 template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<3>&,
                                                     const TinyVector<3>&,
                                                     const TinyVector<3, uint64_t>&);
-
-/// =====================
-
-template <size_t Dimension>
-inline void
-LogicalConnectivityBuilder::_buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
-{
-  static_assert(Dimension <= 3, "unexpected dimension");
-}
-
-template <>
-inline void
-LogicalConnectivityBuilder::_buildBoundaryNodeList(
-  const TinyVector<1, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
-  ConnectivityDescriptor& descriptor)
-{
-  {   // xmin
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = 0;
-    descriptor.addRefItemList(RefNodeList{RefId{0, "XMIN"}, boundary_nodes});
-  }
-
-  {   // xmax
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = cell_size[0];
-    descriptor.addRefItemList(RefNodeList{RefId{1, "XMAX"}, boundary_nodes});
-  }
-}
-
-template <>
-void
-LogicalConnectivityBuilder::_buildBoundaryNodeList(
-  const TinyVector<2, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
-  ConnectivityDescriptor& descriptor)
-{
-  const TinyVector<2, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
-
-  const auto node_number = [&](const TinyVector<2, uint64_t> node_logic_id) {
-    return node_logic_id[0] * node_size[1] + node_logic_id[1];
-  };
-
-  {   // xminymin
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, 0});
-    descriptor.addRefItemList(RefNodeList{RefId{10, "XMINYMIN"}, boundary_nodes});
-  }
-
-  {   // xmaxymin
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], 0});
-    descriptor.addRefItemList(RefNodeList{RefId{11, "XMAXYMIN"}, boundary_nodes});
-  }
-
-  {   // xmaxymax
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], cell_size[1]});
-    descriptor.addRefItemList(RefNodeList{RefId{12, "XMAXYMAX"}, boundary_nodes});
-  }
-
-  {   // xminymax
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, cell_size[1]});
-    descriptor.addRefItemList(RefNodeList{RefId{13, "XMINYMAX"}, boundary_nodes});
-  }
-}
-
-template <>
-void
-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};
-
-  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
-    return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
-  };
-
-  {   // xminyminzmin
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, 0, 0});
-    descriptor.addRefItemList(RefNodeList{RefId{10, "XMINYMINZMIN"}, boundary_nodes});
-  }
-
-  {   // xmaxyminzmin
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], 0, 0});
-    descriptor.addRefItemList(RefNodeList{RefId{11, "XMAXYMINZMIN"}, boundary_nodes});
-  }
-
-  {   // xmaxymaxzmin
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], cell_size[1], 0});
-    descriptor.addRefItemList(RefNodeList{RefId{12, "XMAXYMAXZMIN"}, boundary_nodes});
-  }
-
-  {   // xminymaxzmin
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, cell_size[1], 0});
-    descriptor.addRefItemList(RefNodeList{RefId{13, "XMINYMAXZMIN"}, boundary_nodes});
-  }
-
-  {   // xminyminzmax
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, 0, cell_size[2]});
-    descriptor.addRefItemList(RefNodeList{RefId{14, "XMINYMINZMAX"}, boundary_nodes});
-  }
-
-  {   // xmaxyminzmax
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], 0, cell_size[2]});
-    descriptor.addRefItemList(RefNodeList{RefId{15, "XMAXYMINZMAX"}, boundary_nodes});
-  }
-
-  {   // xmaxymaxzmax
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], cell_size[1], cell_size[2]});
-    descriptor.addRefItemList(RefNodeList{RefId{16, "XMAXYMAXZMAX"}, boundary_nodes});
-  }
-
-  {   // xminymaxzmax
-    Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, cell_size[1], cell_size[2]});
-    descriptor.addRefItemList(RefNodeList{RefId{17, "XMINYMAXZMAX"}, boundary_nodes});
-  }
-}
-
-template <size_t Dimension>
-void
-LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
-{
-  static_assert(Dimension == 3, "unexpected dimension");
-}
-
-template <>
-void
-LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
-                                                   ConnectivityDescriptor& 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;
-  }();
-
-  const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
-  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
-    return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
-  };
-
-  {   // Ridge edges in direction of Z axis
-    const auto add_ref_item_list_along_z = [&](const size_t i, const size_t j, const unsigned ref_id,
-                                               const std::string& ref_name) {
-      Array<EdgeId> boundary_edges(cell_size[2]);
-      size_t l = 0;
-      for (size_t k = 0; k < cell_size[2]; ++k) {
-        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
-
-        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
-        Assert(i_edge != edge_to_id_map.end());
-
-        boundary_edges[l++] = i_edge->second;
-      }
-      Assert(l == cell_size[2]);
-      descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges});
-    };
-
-    add_ref_item_list_along_z(0, 0, 20, "XMINYMIN");
-    add_ref_item_list_along_z(0, cell_size[1], 21, "XMINYMAX");
-    add_ref_item_list_along_z(cell_size[0], cell_size[1], 22, "XMAXYMAX");
-    add_ref_item_list_along_z(cell_size[0], 0, 23, "XMAXYMIN");
-  }
-
-  {   // Ridge edges in direction of Y axis
-    const auto add_ref_item_list_along_y = [&](const size_t i, const size_t k, const unsigned ref_id,
-                                               const std::string& ref_name) {
-      Array<EdgeId> boundary_edges(cell_size[1]);
-      size_t l = 0;
-      for (size_t j = 0; j < cell_size[1]; ++j) {
-        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
-
-        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
-        Assert(i_edge != edge_to_id_map.end());
-
-        boundary_edges[l++] = i_edge->second;
-      }
-      Assert(l == cell_size[1]);
-      descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges});
-    };
-
-    add_ref_item_list_along_y(0, 0, 24, "XMINZMIN");
-    add_ref_item_list_along_y(0, cell_size[2], 25, "XMINZMAX");
-    add_ref_item_list_along_y(cell_size[0], cell_size[2], 26, "XMAXZMAX");
-    add_ref_item_list_along_y(cell_size[0], 0, 27, "XMAXZMIN");
-  }
-
-  {   // Ridge edges in direction of X axis
-    const auto add_ref_item_list_along_x = [&](const size_t j, const size_t k, const unsigned ref_id,
-                                               const std::string& ref_name) {
-      Array<EdgeId> boundary_edges(cell_size[0]);
-      size_t l = 0;
-      for (size_t i = 0; i < cell_size[0]; ++i) {
-        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
-
-        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
-        Assert(i_edge != edge_to_id_map.end());
-
-        boundary_edges[l++] = i_edge->second;
-      }
-      Assert(l == cell_size[0]);
-      descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges});
-    };
-
-    add_ref_item_list_along_x(0, 0, 28, "YMINZMIN");
-    add_ref_item_list_along_x(0, cell_size[2], 29, "YMINZMAX");
-    add_ref_item_list_along_x(cell_size[1], cell_size[2], 30, "YMAXZMAX");
-    add_ref_item_list_along_x(cell_size[1], 0, 31, "YMAXZMIN");
-  }
-}
-
-template <size_t Dimension>
-void
-LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
-{
-  static_assert(Dimension >= 2 and Dimension <= 3, "unexpected dimension");
-}
-
-template <>
-void
-LogicalConnectivityBuilder::_buildBoundaryFaceList(
-  const TinyVector<2, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
-  ConnectivityDescriptor& descriptor)
-{
-  const auto cell_number = [&](const TinyVector<2, uint64_t> cell_logic_id) {
-    return cell_logic_id[0] * cell_size[1] + cell_logic_id[1];
-  };
-
-  {   // xmin
-    const size_t i = 0;
-    Array<FaceId> boundary_faces(cell_size[1]);
-    for (size_t j = 0; j < cell_size[1]; ++j) {
-      constexpr size_t left_face = 3;
-
-      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][left_face];
-
-      boundary_faces[j] = face_id;
-    }
-    descriptor.addRefItemList(RefFaceList{RefId{0, "XMIN"}, boundary_faces});
-  }
-
-  {   // xmax
-    const size_t i = cell_size[0] - 1;
-    Array<FaceId> boundary_faces(cell_size[1]);
-    for (size_t j = 0; j < cell_size[1]; ++j) {
-      constexpr size_t right_face = 1;
-
-      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][right_face];
-
-      boundary_faces[j] = face_id;
-    }
-    descriptor.addRefItemList(RefFaceList{RefId{1, "XMAX"}, boundary_faces});
-  }
-
-  {   // ymin
-    const size_t j = 0;
-    Array<FaceId> boundary_faces(cell_size[0]);
-    for (size_t i = 0; i < cell_size[0]; ++i) {
-      constexpr size_t bottom_face = 0;
-
-      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][bottom_face];
-
-      boundary_faces[i] = face_id;
-    }
-    descriptor.addRefItemList(RefFaceList{RefId{2, "YMIN"}, boundary_faces});
-  }
-
-  {   // ymax
-    const size_t j = cell_size[1] - 1;
-    Array<FaceId> boundary_faces(cell_size[0]);
-    for (size_t i = 0; i < cell_size[0]; ++i) {
-      constexpr size_t top_face = 2;
-
-      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
-      const size_t face_id = descriptor.cell_to_face_vector[cell_id][top_face];
-
-      boundary_faces[i] = face_id;
-    }
-    descriptor.addRefItemList(RefFaceList{RefId{3, "YMAX"}, boundary_faces});
-  }
-}
-
-template <>
-void
-LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
-                                                   ConnectivityDescriptor& descriptor)
-{
-  using Face = ConnectivityFace<3>;
-
-  const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
-    std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
-    for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-      const auto& node_vector                                          = descriptor.face_to_node_vector[l];
-      face_to_id_map[Face(node_vector, descriptor.node_number_vector)] = l;
-    }
-    return face_to_id_map;
-  }();
-
-  const std::unordered_map<int, FaceId> face_number_id_map = [&] {
-    std::unordered_map<int, FaceId> face_number_id_map;
-    for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
-      face_number_id_map[descriptor.face_number_vector[l]] = l;
-    }
-    Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
-    return face_number_id_map;
-  }();
-
-  const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
-  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
-    return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
-  };
-
-  {   // faces orthogonal to X
-    const auto add_ref_item_list_for_x = [&](const size_t i, const unsigned ref_id, const std::string& ref_name) {
-      Array<FaceId> boundary_faces(cell_size[1] * cell_size[2]);
-      size_t l = 0;
-      for (size_t j = 0; j < cell_size[1]; ++j) {
-        for (size_t k = 0; k < cell_size[2]; ++k) {
-          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
-          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k + 1});
-          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
-
-          auto i_face =
-            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
-          Assert(i_face != face_to_id_map.end());
-
-          boundary_faces[l++] = i_face->second;
-        }
-      }
-      Assert(l == cell_size[1] * cell_size[2]);
-      descriptor.addRefItemList(RefFaceList{RefId{ref_id, ref_name}, boundary_faces});
-    };
-
-    add_ref_item_list_for_x(0, 0, "XMIN");
-    add_ref_item_list_for_x(cell_size[0], 1, "XMAX");
-  }
-
-  {   // faces orthogonal to Y
-    const auto add_ref_item_list_for_y = [&](const size_t j, const unsigned ref_id, const std::string& ref_name) {
-      Array<FaceId> boundary_faces(cell_size[0] * cell_size[2]);
-      size_t l = 0;
-      for (size_t i = 0; i < cell_size[0]; ++i) {
-        for (size_t k = 0; k < cell_size[2]; ++k) {
-          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
-          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k + 1});
-          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
-
-          auto i_face =
-            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
-          Assert(i_face != face_to_id_map.end());
-
-          boundary_faces[l++] = i_face->second;
-        }
-      }
-      Assert(l == cell_size[0] * cell_size[2]);
-      descriptor.addRefItemList(RefFaceList{RefId{ref_id, ref_name}, boundary_faces});
-    };
-
-    add_ref_item_list_for_y(0, 2, "YMIN");
-    add_ref_item_list_for_y(cell_size[1], 3, "YMAX");
-  }
-
-  {   // faces orthogonal to Z
-    const auto add_ref_item_list_for_z = [&](const size_t k, const unsigned ref_id, const std::string& ref_name) {
-      Array<FaceId> boundary_faces(cell_size[0] * cell_size[1]);
-      size_t l = 0;
-      for (size_t i = 0; i < cell_size[0]; ++i) {
-        for (size_t j = 0; j < cell_size[1]; ++j) {
-          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
-          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
-          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k});
-          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
-
-          auto i_face =
-            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
-          Assert(i_face != face_to_id_map.end());
-
-          boundary_faces[l++] = i_face->second;
-        }
-      }
-      Assert(l == cell_size[0] * cell_size[1]);
-      descriptor.addRefItemList(RefFaceList{RefId{ref_id, ref_name}, boundary_faces});
-    };
-
-    add_ref_item_list_for_z(0, 4, "ZMIN");
-    add_ref_item_list_for_z(cell_size[2], 5, "ZMAX");
-  }
-}
-
-template <>
-void
-LogicalConnectivityBuilder::_buildConnectivity(
-  const TinyVector<1, uint64_t>& cell_size)   // clazy:exclude=function-args-by-value
-{
-  const size_t number_of_cells = cell_size[0];
-  const size_t number_of_nodes = cell_size[0] + 1;
-
-  ConnectivityDescriptor descriptor;
-  descriptor.node_number_vector.resize(number_of_nodes);
-  for (size_t i = 0; i < number_of_nodes; ++i) {
-    descriptor.node_number_vector[i] = i;
-  }
-
-  descriptor.cell_number_vector.resize(number_of_cells);
-  for (size_t i = 0; i < number_of_cells; ++i) {
-    descriptor.cell_number_vector[i] = i;
-  }
-
-  descriptor.cell_type_vector.resize(number_of_cells);
-  std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Line);
-
-  descriptor.cell_to_node_vector.resize(number_of_cells);
-  constexpr size_t nb_node_per_cell = 2;
-  for (size_t j = 0; j < number_of_cells; ++j) {
-    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
-    for (size_t r = 0; r < nb_node_per_cell; ++r) {
-      descriptor.cell_to_node_vector[j][0] = j;
-      descriptor.cell_to_node_vector[j][1] = j + 1;
-    }
-  }
-
-  this->_buildBoundaryNodeList(cell_size, descriptor);
-
-  descriptor.cell_owner_vector.resize(number_of_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 <>
-void
-LogicalConnectivityBuilder::_buildConnectivity(
-  const TinyVector<2, uint64_t>& cell_size)   // clazy:exclude=function-args-by-value
-{
-  constexpr size_t Dimension = 2;
-
-  const TinyVector<Dimension, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
-
-  const size_t number_of_cells = cell_size[0] * cell_size[1];
-  const size_t number_of_nodes = node_size[0] * node_size[1];
-
-  ConnectivityDescriptor descriptor;
-  descriptor.node_number_vector.resize(number_of_nodes);
-  for (size_t i = 0; i < number_of_nodes; ++i) {
-    descriptor.node_number_vector[i] = i;
-  }
-
-  descriptor.cell_number_vector.resize(number_of_cells);
-  for (size_t i = 0; i < number_of_cells; ++i) {
-    descriptor.cell_number_vector[i] = i;
-  }
-
-  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_number = [&](const TinyVector<Dimension, uint64_t> node_logic_id) {
-    return node_logic_id[0] * node_size[1] + node_logic_id[1];
-  };
-
-  const auto cell_logic_id = [&](size_t j) {
-    const uint64_t j0 = j / cell_size[1];
-    const uint64_t j1 = j % cell_size[1];
-    return TinyVector<Dimension, uint64_t>{j0, j1};
-  };
-
-  descriptor.cell_to_node_vector.resize(number_of_cells);
-  constexpr size_t nb_node_per_cell = 1 << Dimension;
-  for (size_t j = 0; j < number_of_cells; ++j) {
-    TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
-    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
-    for (size_t r = 0; r < nb_node_per_cell; ++r) {
-      descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0});
-      descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0});
-      descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1});
-      descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1});
-    }
-  }
-
-  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
-
-  this->_buildBoundaryNodeList(cell_size, descriptor);
-  this->_buildBoundaryFaceList(cell_size, descriptor);
-
-  descriptor.cell_owner_vector.resize(number_of_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 = Connectivity<Dimension>::build(descriptor);
-}
-
-template <>
-void
-LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& cell_size)
-{
-  constexpr size_t Dimension = 3;
-
-  const TinyVector<Dimension, uint64_t> node_size = [&] {
-    TinyVector node_size{cell_size};
-    for (size_t i = 0; i < Dimension; ++i) {
-      node_size[i] += 1;
-    }
-    return node_size;
-  }();
-
-  auto count_items = [](auto&& v) {
-    using V_Type = std::decay_t<decltype(v)>;
-    static_assert(is_tiny_vector_v<V_Type>);
-
-    size_t sum = v[0];
-    for (size_t i = 1; i < Dimension; ++i) {
-      sum *= v[i];
-    }
-    return sum;
-  };
-
-  const size_t number_of_cells = count_items(cell_size);
-  const size_t number_of_nodes = count_items(node_size);
-
-  ConnectivityDescriptor descriptor;
-  descriptor.node_number_vector.resize(number_of_nodes);
-  for (size_t i = 0; i < number_of_nodes; ++i) {
-    descriptor.node_number_vector[i] = i;
-  }
-
-  descriptor.cell_number_vector.resize(number_of_cells);
-  for (size_t i = 0; i < number_of_cells; ++i) {
-    descriptor.cell_number_vector[i] = i;
-  }
-
-  descriptor.cell_type_vector.resize(number_of_cells);
-  std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Hexahedron);
-
-  const auto cell_logic_id = [&](size_t j) {
-    const size_t slice1  = cell_size[1] * cell_size[2];
-    const size_t& slice2 = cell_size[2];
-    const uint64_t j0    = j / slice1;
-    const uint64_t j1    = (j - j0 * slice1) / slice2;
-    const uint64_t j2    = j - (j0 * slice1 + j1 * slice2);
-    return TinyVector<Dimension, uint64_t>{j0, j1, j2};
-  };
-
-  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];
-  };
-
-  descriptor.cell_to_node_vector.resize(number_of_cells);
-  constexpr size_t nb_node_per_cell = 1 << Dimension;
-  for (size_t j = 0; j < number_of_cells; ++j) {
-    TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
-    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
-    for (size_t r = 0; r < nb_node_per_cell; ++r) {
-      static_assert(Dimension == 3, "unexpected dimension");
-      descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 0});
-      descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 0});
-      descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 0});
-      descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 0});
-      descriptor.cell_to_node_vector[j][4] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 1});
-      descriptor.cell_to_node_vector[j][5] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 1});
-      descriptor.cell_to_node_vector[j][6] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 1});
-      descriptor.cell_to_node_vector[j][7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
-    }
-  }
-
-  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
-  ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
-
-  this->_buildBoundaryNodeList(cell_size, descriptor);
-  this->_buildBoundaryEdgeList(cell_size, descriptor);
-  this->_buildBoundaryFaceList(cell_size, descriptor);
-
-  descriptor.cell_owner_vector.resize(number_of_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 = Connectivity<Dimension>::build(descriptor);
-}
-
-template <size_t Dimension>
-LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size)
-{
-  if (parallel::rank() == 0) {
-    this->_buildConnectivity(size);
-  }
-}
-
-template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<1, uint64_t>&);
-
-template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<2, uint64_t>&);
-
-template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<3, uint64_t>&);
diff --git a/src/mesh/CartesianMeshBuilder.hpp b/src/mesh/CartesianMeshBuilder.hpp
index 9ed1ea81cd4c3dbf0863fb50fe423cd49aee2d0b..7b32b55b64d021c655ff9c35245ece4b31b2ed76 100644
--- a/src/mesh/CartesianMeshBuilder.hpp
+++ b/src/mesh/CartesianMeshBuilder.hpp
@@ -29,25 +29,4 @@ class CartesianMeshBuilder : public MeshBuilderBase
   ~CartesianMeshBuilder() = default;
 };
 
-class LogicalConnectivityBuilder : public ConnectivityBuilderBase
-{
- private:
-  template <size_t Dimension>
-  void _buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
-
-  template <size_t Dimension>
-  void _buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
-
-  template <size_t Dimension>
-  void _buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
-
-  template <size_t Dimension>
-  void _buildConnectivity(const TinyVector<Dimension, uint64_t>& size);
-
- public:
-  template <size_t Dimension>
-  LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size);
-  ~LogicalConnectivityBuilder() = default;
-};
-
 #endif   // CARTESIAN_MESH_BUILDER_HPP
diff --git a/src/mesh/ConnectivityBuilderBase.cpp b/src/mesh/ConnectivityBuilderBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f5d47e25003a710e1e2c18ba6b772ffd04bb532f
--- /dev/null
+++ b/src/mesh/ConnectivityBuilderBase.cpp
@@ -0,0 +1,416 @@
+#include <mesh/ConnectivityBuilderBase.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityDescriptor.hpp>
+#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ItemId.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <vector>
+
+template <size_t Dimension>
+void
+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>;
+  using Face         = ConnectivityFace<Dimension>;
+
+  const auto& node_number_vector = descriptor.node_number_vector;
+  Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
+  std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
+  for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
+    const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+
+    if constexpr (Dimension == 2) {
+      switch (descriptor.cell_type_vector[j]) {
+      case CellType::Triangle: {
+        cell_nb_faces[j] = 3;
+        // face 0
+        Face f0({cell_nodes[1], cell_nodes[2]}, node_number_vector);
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes[2], cell_nodes[0]}, node_number_vector);
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes[0], cell_nodes[1]}, node_number_vector);
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+        break;
+      }
+      case CellType::Quadrangle: {
+        cell_nb_faces[j] = 4;
+        // face 0
+        Face f0({cell_nodes[0], cell_nodes[1]}, node_number_vector);
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes[1], cell_nodes[2]}, node_number_vector);
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes[2], cell_nodes[3]}, node_number_vector);
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+        // face 3
+        Face f3({cell_nodes[3], cell_nodes[0]}, node_number_vector);
+        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+        break;
+      }
+      default: {
+        std::ostringstream error_msg;
+        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
+        throw UnexpectedError(error_msg.str());
+      }
+      }
+    } else if constexpr (Dimension == 3) {
+      switch (descriptor.cell_type_vector[j]) {
+      case CellType::Hexahedron: {
+        // face 0
+        Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes[4], cell_nodes[5], cell_nodes[6], cell_nodes[7]}, node_number_vector);
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+        // face 3
+        Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
+        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+
+        // face 4
+        Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
+        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+
+        // face 5
+        Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
+        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+
+        cell_nb_faces[j] = 6;
+        break;
+      }
+      case CellType::Tetrahedron: {
+        cell_nb_faces[j] = 4;
+        // face 0
+        Face f0({cell_nodes[1], cell_nodes[2], cell_nodes[3]}, node_number_vector);
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes[0], cell_nodes[3], cell_nodes[2]}, node_number_vector);
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes[0], cell_nodes[1], cell_nodes[3]}, node_number_vector);
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+        // face 3
+        Face f3({cell_nodes[0], cell_nodes[2], cell_nodes[1]}, node_number_vector);
+        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+        break;
+      }
+      case CellType::Pyramid: {
+        cell_nb_faces[j] = cell_nodes.size();
+        std::vector<unsigned int> base_nodes;
+        std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
+
+        // base face
+        {
+          Face base_face(base_nodes, node_number_vector);
+          face_cells_map[base_face].emplace_back(std::make_tuple(j, 0, base_face.reversed()));
+        }
+        // side faces
+        const auto pyramid_vertex = cell_nodes[cell_nodes.size() - 1];
+        for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
+          Face side_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], pyramid_vertex},
+                         node_number_vector);
+          face_cells_map[side_face].emplace_back(std::make_tuple(j, i_node + 1, side_face.reversed()));
+        }
+        break;
+      }
+      case CellType::Diamond: {
+        cell_nb_faces[j] = 2 * (cell_nodes.size() - 2);
+        std::vector<unsigned int> base_nodes;
+        std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
+
+        {   // top faces
+          const auto top_vertex = cell_nodes[cell_nodes.size() - 1];
+          for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
+            Face top_face({base_nodes[i_node], base_nodes[(i_node + 1) % base_nodes.size()], top_vertex},
+                          node_number_vector);
+            face_cells_map[top_face].emplace_back(std::make_tuple(j, i_node, top_face.reversed()));
+          }
+        }
+
+        {   // bottom faces
+          const auto bottom_vertex = cell_nodes[0];
+          for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
+            Face bottom_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], bottom_vertex},
+                             node_number_vector);
+            face_cells_map[bottom_face].emplace_back(
+              std::make_tuple(j, i_node + base_nodes.size(), bottom_face.reversed()));
+          }
+        }
+        break;
+      }
+      default: {
+        std::ostringstream error_msg;
+        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+        throw UnexpectedError(error_msg.str());
+      }
+      }
+    }
+  }
+
+  {
+    descriptor.cell_to_face_vector.resize(descriptor.cell_to_node_vector.size());
+    for (CellId j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
+      descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
+    }
+    FaceId l = 0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& cells_vector = face_cells_vector.second;
+      for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
+        const auto& [cell_number, cell_local_face, reversed]         = cells_vector[lj];
+        descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
+      }
+      ++l;
+    }
+  }
+
+  {
+    descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
+    for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
+      descriptor.cell_face_is_reversed_vector[j] = Array<bool>(cell_nb_faces[j]);
+    }
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& cells_vector = face_cells_vector.second;
+      for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
+        const auto& [cell_number, cell_local_face, reversed]                  = cells_vector[lj];
+        descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
+      }
+    }
+  }
+
+  {
+    descriptor.face_to_node_vector.resize(face_cells_map.size());
+    int l = 0;
+    for (const auto& face_info : face_cells_map) {
+      const Face& face                  = face_info.first;
+      descriptor.face_to_node_vector[l] = face.nodeIdList();
+      ++l;
+    }
+  }
+
+  {
+    // Face numbers may change if numbers are provided in the file
+    descriptor.face_number_vector.resize(face_cells_map.size());
+    for (size_t l = 0; l < face_cells_map.size(); ++l) {
+      descriptor.face_number_vector[l] = l;
+    }
+  }
+}
+
+template <size_t Dimension>
+void
+ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
+{
+  static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
+  using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
+  using Edge         = ConnectivityFace<2>;
+
+  const auto& node_number_vector = descriptor.node_number_vector;
+  Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
+  std::map<Edge, std::vector<FaceEdgeInfo>> edge_faces_map;
+  for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
+    const auto& face_nodes = descriptor.face_to_node_vector[l];
+
+    face_nb_edges[l] = face_nodes.size();
+    for (size_t r = 0; r < face_nodes.size() - 1; ++r) {
+      Edge e({face_nodes[r], face_nodes[r + 1]}, node_number_vector);
+      edge_faces_map[e].emplace_back(std::make_tuple(l, r, e.reversed()));
+    }
+    {
+      Edge e({face_nodes[face_nodes.size() - 1], face_nodes[0]}, node_number_vector);
+      edge_faces_map[e].emplace_back(std::make_tuple(l, face_nodes.size() - 1, e.reversed()));
+    }
+  }
+
+  std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_id_map;
+  {
+    descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
+    for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
+      descriptor.face_to_edge_vector[l].resize(face_nb_edges[l]);
+    }
+    EdgeId e = 0;
+    for (const auto& edge_faces_vector : edge_faces_map) {
+      const auto& faces_vector = edge_faces_vector.second;
+      for (unsigned short l = 0; l < faces_vector.size(); ++l) {
+        const auto& [face_number, face_local_edge, reversed]         = faces_vector[l];
+        descriptor.face_to_edge_vector[face_number][face_local_edge] = e;
+      }
+      edge_id_map[edge_faces_vector.first] = e;
+      ++e;
+    }
+  }
+
+  {
+    descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
+    for (FaceId j = 0; j < descriptor.face_edge_is_reversed_vector.size(); ++j) {
+      descriptor.face_edge_is_reversed_vector[j] = Array<bool>(face_nb_edges[j]);
+    }
+    for (const auto& edge_faces_vector : edge_faces_map) {
+      const auto& faces_vector = edge_faces_vector.second;
+      for (unsigned short lj = 0; lj < faces_vector.size(); ++lj) {
+        const auto& [face_number, face_local_edge, reversed]                  = faces_vector[lj];
+        descriptor.face_edge_is_reversed_vector[face_number][face_local_edge] = reversed;
+      }
+    }
+  }
+
+  {
+    descriptor.edge_to_node_vector.resize(edge_faces_map.size());
+    int e = 0;
+    for (const auto& edge_info : edge_faces_map) {
+      const Edge& edge                  = edge_info.first;
+      descriptor.edge_to_node_vector[e] = edge.nodeIdList();
+      ++e;
+    }
+  }
+
+  {
+    // Edge numbers may change if numbers are provided in the file
+    descriptor.edge_number_vector.resize(edge_faces_map.size());
+    for (size_t e = 0; e < edge_faces_map.size(); ++e) {
+      descriptor.edge_number_vector[e] = e;
+    }
+  }
+
+  {
+    descriptor.cell_to_edge_vector.reserve(descriptor.cell_to_node_vector.size());
+    for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
+      const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+
+      switch (descriptor.cell_type_vector[j]) {
+      case CellType::Tetrahedron: {
+        constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
+        std::vector<unsigned int> cell_edge_vector;
+        cell_edge_vector.reserve(6);
+        for (int i_edge = 0; i_edge < 6; ++i_edge) {
+          const auto e = local_edge[i_edge];
+          Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            throw NormalError("could not find this edge");
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+        break;
+      }
+      case CellType::Hexahedron: {
+        constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
+                                           {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
+        std::vector<unsigned int> cell_edge_vector;
+        cell_edge_vector.reserve(12);
+        for (int i_edge = 0; i_edge < 12; ++i_edge) {
+          const auto e = local_edge[i_edge];
+          Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            throw NormalError("could not find this edge");
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+        break;
+      }
+      case CellType::Pyramid: {
+        const size_t number_of_edges = 2 * cell_nodes.size();
+        std::vector<unsigned int> base_nodes;
+        std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
+
+        std::vector<unsigned int> cell_edge_vector;
+        cell_edge_vector.reserve(number_of_edges);
+        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+          Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            throw NormalError("could not find this edge");
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+
+        const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+          Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            throw NormalError("could not find this edge");
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+        break;
+      }
+      case CellType::Diamond: {
+        const size_t number_of_edges = 3 * cell_nodes.size();
+        std::vector<unsigned int> base_nodes;
+        std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
+
+        std::vector<unsigned int> cell_edge_vector;
+        cell_edge_vector.reserve(number_of_edges);
+        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+          Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            throw NormalError("could not find this edge");
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+
+        const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+          Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            throw NormalError("could not find this edge");
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+
+        const unsigned int bottom_vertex = cell_nodes[0];
+        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+          Edge edge{{base_nodes[i_edge], bottom_vertex}, node_number_vector};
+          auto i = edge_id_map.find(edge);
+          if (i == edge_id_map.end()) {
+            throw NormalError("could not find this edge");
+          }
+          cell_edge_vector.push_back(i->second);
+        }
+
+        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+
+        break;
+      }
+      default: {
+        std::stringstream error_msg;
+        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+        throw NotImplementedError(error_msg.str());
+      }
+      }
+    }
+  }
+}
+
+template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
+template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
+
+template void ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
+  ConnectivityDescriptor& descriptor);
diff --git a/src/mesh/ConnectivityBuilderBase.hpp b/src/mesh/ConnectivityBuilderBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9b0d28a7a9389093348648573c85875f1e4a7de
--- /dev/null
+++ b/src/mesh/ConnectivityBuilderBase.hpp
@@ -0,0 +1,226 @@
+#ifndef CONNECTIVITY_BUILDER_BASE_HPP
+#define CONNECTIVITY_BUILDER_BASE_HPP
+
+#include <mesh/IConnectivity.hpp>
+
+#include <mesh/ItemId.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <memory>
+#include <vector>
+
+class ConnectivityDescriptor;
+
+class ConnectivityBuilderBase
+{
+ protected:
+  template <size_t Dimension>
+  class ConnectivityFace;
+
+  std::shared_ptr<const IConnectivity> m_connectivity;
+
+  template <size_t Dimension>
+  static void _computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
+
+  template <size_t Dimension>
+  static void _computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
+
+ public:
+  std::shared_ptr<const IConnectivity>
+  connectivity() const
+  {
+    return m_connectivity;
+  }
+
+  ConnectivityBuilderBase()  = default;
+  ~ConnectivityBuilderBase() = default;
+};
+
+template <>
+class ConnectivityBuilderBase::ConnectivityFace<2>
+{
+ public:
+  friend struct Hash;
+
+  struct Hash
+  {
+    size_t
+    operator()(const ConnectivityFace& f) const
+    {
+      size_t hash = 0;
+      hash ^= std::hash<unsigned int>()(f.m_node0_id);
+      hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
+      return hash;
+    }
+  };
+
+ private:
+  const std::vector<int>& m_node_number_vector;
+
+  unsigned int m_node0_id;
+  unsigned int m_node1_id;
+
+  bool m_reversed;
+
+ public:
+  std::vector<unsigned int>
+  nodeIdList() const
+  {
+    return {m_node0_id, m_node1_id};
+  }
+
+  bool
+  reversed() const
+  {
+    return m_reversed;
+  }
+
+  PUGS_INLINE
+  bool
+  operator==(const ConnectivityFace& f) const
+  {
+    return ((m_node0_id == f.m_node0_id) and (m_node1_id == f.m_node1_id));
+  }
+
+  PUGS_INLINE
+  bool
+  operator<(const ConnectivityFace& f) const
+  {
+    return ((m_node_number_vector[m_node0_id] < m_node_number_vector[f.m_node0_id]) or
+            ((m_node_number_vector[m_node0_id] == m_node_number_vector[f.m_node0_id]) and
+             (m_node_number_vector[m_node1_id] < m_node_number_vector[f.m_node1_id])));
+  }
+
+  PUGS_INLINE
+  ConnectivityFace(const std::vector<unsigned int>& node_id_list, const std::vector<int>& node_number_vector)
+    : m_node_number_vector(node_number_vector)
+  {
+    Assert(node_id_list.size() == 2);
+
+    if (m_node_number_vector[node_id_list[0]] < m_node_number_vector[node_id_list[1]]) {
+      m_node0_id = node_id_list[0];
+      m_node1_id = node_id_list[1];
+      m_reversed = false;
+    } else {
+      m_node0_id = node_id_list[1];
+      m_node1_id = node_id_list[0];
+      m_reversed = true;
+    }
+  }
+
+  PUGS_INLINE
+  ConnectivityFace(const ConnectivityFace&) = default;
+
+  PUGS_INLINE
+  ~ConnectivityFace() = default;
+};
+
+template <>
+class ConnectivityBuilderBase::ConnectivityFace<3>
+{
+ public:
+  friend struct Hash;
+
+  struct Hash
+  {
+    size_t
+    operator()(const ConnectivityFace& f) const
+    {
+      size_t hash = 0;
+      for (size_t i = 0; i < f.m_node_id_list.size(); ++i) {
+        hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
+      }
+      return hash;
+    }
+  };
+
+ private:
+  bool m_reversed;
+  std::vector<NodeId::base_type> m_node_id_list;
+  const std::vector<int>& m_node_number_vector;
+
+  PUGS_INLINE
+  std::vector<unsigned int>
+  _sort(const std::vector<unsigned int>& node_list)
+  {
+    const auto min_id = std::min_element(node_list.begin(), node_list.end());
+    const int shift   = std::distance(node_list.begin(), min_id);
+
+    std::vector<unsigned int> rotated_node_list(node_list.size());
+    if (node_list[(shift + 1) % node_list.size()] > node_list[(shift + node_list.size() - 1) % node_list.size()]) {
+      for (size_t i = 0; i < node_list.size(); ++i) {
+        rotated_node_list[i] = node_list[(shift + node_list.size() - i) % node_list.size()];
+        m_reversed           = true;
+      }
+    } else {
+      for (size_t i = 0; i < node_list.size(); ++i) {
+        rotated_node_list[i] = node_list[(shift + i) % node_list.size()];
+      }
+    }
+
+    return rotated_node_list;
+  }
+
+ public:
+  PUGS_INLINE
+  const bool&
+  reversed() const
+  {
+    return m_reversed;
+  }
+
+  PUGS_INLINE
+  const std::vector<unsigned int>&
+  nodeIdList() const
+  {
+    return m_node_id_list;
+  }
+
+  PUGS_INLINE
+  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list, const std::vector<int>& node_number_vector)
+    : m_reversed(false), m_node_id_list(_sort(given_node_id_list)), m_node_number_vector(node_number_vector)
+  {
+    ;
+  }
+
+ public:
+  bool
+  operator==(const ConnectivityFace& f) const
+  {
+    if (m_node_id_list.size() == f.nodeIdList().size()) {
+      for (size_t j = 0; j < m_node_id_list.size(); ++j) {
+        if (m_node_id_list[j] != f.nodeIdList()[j]) {
+          return false;
+        }
+      }
+      return true;
+    }
+    return false;
+  }
+
+  PUGS_INLINE
+  bool
+  operator<(const ConnectivityFace& f) const
+  {
+    const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
+    for (size_t i = 0; i < min_nb_nodes; ++i) {
+      if (m_node_id_list[i] < f.m_node_id_list[i])
+        return true;
+      if (m_node_id_list[i] != f.m_node_id_list[i])
+        return false;
+    }
+    return m_node_id_list.size() < f.m_node_id_list.size();
+  }
+
+  PUGS_INLINE
+  ConnectivityFace(const ConnectivityFace&) = default;
+
+  PUGS_INLINE
+  ConnectivityFace() = delete;
+
+  PUGS_INLINE
+  ~ConnectivityFace() = default;
+};
+
+#endif   // CONNECTIVITY_BUILDER_BASE_HPP
diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 585982245f568a590f78217996c4fe64dc8eda5f..8bcbe9bad8e15cbef1dfa449c0d6defc115f126f 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -5,7 +5,6 @@
 #include <mesh/ConnectivityDispatcher.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
@@ -227,15 +226,12 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Con
 
 template <size_t Dimension>
 void
-DiamondDualConnectivityBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
+DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivity& i_primal_connectivity)
 {
-  static_assert(Dimension <= 3, "invalid mesh dimension");
+  static_assert(Dimension <= 3, "invalid connectivity dimension");
   using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
 
-  const MeshType& primal_mesh = dynamic_cast<const MeshType&>(*p_mesh);
-
-  const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
+  const ConnectivityType& primal_connectivity = dynamic_cast<const ConnectivityType&>(i_primal_connectivity);
 
   ConnectivityDescriptor diamond_descriptor;
 
@@ -477,129 +473,23 @@ DiamondDualConnectivityBuilder::_buildDiamondMeshFrom(const std::shared_ptr<cons
   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 NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
-  MeshData<MeshType> primal_mesh_data{primal_mesh};
-  const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
-
-  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_mesh.numberOfCells(); ++i_cell) {
-    diamond_xr[i_node++] = primal_xj[i_cell];
-  }
-
-  m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
-}
-
-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::DiamondDualConnectivityBuilder(const IConnectivity& connectivity)
 {
-  DiamondDualConnectivityBuilder builder(p_mesh);
-
-  switch (p_mesh->dimension()) {
+  switch (connectivity.dimension()) {
   case 1: {
-    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);
+    this->_buildDiamondConnectivityFrom<1>(connectivity);
     break;
   }
   case 2: {
-    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);
+    this->_buildDiamondConnectivityFrom<2>(connectivity);
     break;
   }
   case 3: {
-    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);
+    this->_buildDiamondConnectivityFrom<3>(connectivity);
     break;
   }
   default: {
-    throw UnexpectedError("invalid mesh dimension: " + std::to_string(p_mesh->dimension()));
+    throw UnexpectedError("invalid connectivity dimension: " + std::to_string(connectivity.dimension()));
   }
   }
 }
diff --git a/src/mesh/DiamondDualConnectivityBuilder.hpp b/src/mesh/DiamondDualConnectivityBuilder.hpp
index a2fe946b934d5a1e620023f14ffbbab643df3dd1..670372df9e53d3bfa8dd3f1b03545d9264b41dc4 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.hpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.hpp
@@ -1,17 +1,12 @@
 #ifndef DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
 #define DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
 
-#include <mesh/IMesh.hpp>
-#include <mesh/MeshBuilderBase.hpp>
+#include <mesh/ConnectivityBuilderBase.hpp>
 
 #include <memory>
 
 template <size_t>
 class Connectivity;
-
-template <typename ConnectivityType>
-class Mesh;
-
 class ConnectivityDescriptor;
 
 class DiamondDualConnectivityBuilder : public ConnectivityBuilderBase
@@ -21,23 +16,11 @@ class DiamondDualConnectivityBuilder : public ConnectivityBuilderBase
   void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
 
   template <size_t Dimension>
-  void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
+  void _buildDiamondConnectivityFrom(const IConnectivity&);
 
  public:
-  DiamondDualConnectivityBuilder(const std::shared_ptr<const IMesh>&);
+  DiamondDualConnectivityBuilder(const IConnectivity&);
   ~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.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c5416939c4c81dd67114b30a720ee8a2a258a0c8
--- /dev/null
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -0,0 +1,120 @@
+#include <mesh/DiamondDualMeshBuilder.hpp>
+
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityDescriptor.hpp>
+#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/RefId.hpp>
+#include <utils/Array.hpp>
+#include <utils/Messenger.hpp>
+
+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 NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
+  MeshData<MeshType> primal_mesh_data{primal_mesh};
+  const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
+
+  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_mesh.numberOfCells(); ++i_cell) {
+    diamond_xr[i_node++] = primal_xj[i_cell];
+  }
+
+  m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
+}
+
+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)
+{
+  switch (p_mesh->dimension()) {
+  case 1: {
+    using ConnectivityType = Connectivity<1>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+    DiamondDualConnectivityBuilder builder(mesh->connectivity());
+
+    std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+    this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+    break;
+  }
+  case 2: {
+    using ConnectivityType = Connectivity<2>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+    DiamondDualConnectivityBuilder builder(mesh->connectivity());
+
+    std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+    this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+    break;
+  }
+  case 3: {
+    using ConnectivityType = Connectivity<3>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+    DiamondDualConnectivityBuilder builder(mesh->connectivity());
+
+    std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+    this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+    break;
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension: " + std::to_string(p_mesh->dimension()));
+  }
+  }
+}
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ded3cfb0f40ca4d1b1bd9dc06f7981391f9401c6
--- /dev/null
+++ b/src/mesh/DiamondDualMeshBuilder.hpp
@@ -0,0 +1,26 @@
+#ifndef DIAMOND_DUAL_MESH_BUILDER_HPP
+#define DIAMOND_DUAL_MESH_BUILDER_HPP
+
+#include <mesh/MeshBuilderBase.hpp>
+
+#include <memory>
+
+template <size_t>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+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_MESH_BUILDER_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 6ca41953f80673e9784ef6f8e238ddaa1d8e101f..5e65fcbb74b655acae2cbd428df35f35fc675973 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -4,6 +4,7 @@
 
 #include <mesh/CellType.hpp>
 #include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityBuilderBase.hpp>
 #include <mesh/ConnectivityDispatcher.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
diff --git a/src/mesh/LogicalConnectivityBuilder.cpp b/src/mesh/LogicalConnectivityBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c44f374fd5b2a64ed52ae05e3ede2aa012532319
--- /dev/null
+++ b/src/mesh/LogicalConnectivityBuilder.cpp
@@ -0,0 +1,641 @@
+#include <mesh/LogicalConnectivityBuilder.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityDescriptor.hpp>
+#include <mesh/RefId.hpp>
+#include <utils/Array.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension>
+inline void
+LogicalConnectivityBuilder::_buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+{
+  static_assert(Dimension <= 3, "unexpected dimension");
+}
+
+template <>
+inline void
+LogicalConnectivityBuilder::_buildBoundaryNodeList(
+  const TinyVector<1, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
+  ConnectivityDescriptor& descriptor)
+{
+  {   // xmin
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = 0;
+    descriptor.addRefItemList(RefNodeList{RefId{0, "XMIN"}, boundary_nodes});
+  }
+
+  {   // xmax
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = cell_size[0];
+    descriptor.addRefItemList(RefNodeList{RefId{1, "XMAX"}, boundary_nodes});
+  }
+}
+
+template <>
+void
+LogicalConnectivityBuilder::_buildBoundaryNodeList(
+  const TinyVector<2, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
+  ConnectivityDescriptor& descriptor)
+{
+  const TinyVector<2, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
+
+  const auto node_number = [&](const TinyVector<2, uint64_t> node_logic_id) {
+    return node_logic_id[0] * node_size[1] + node_logic_id[1];
+  };
+
+  {   // xminymin
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({0, 0});
+    descriptor.addRefItemList(RefNodeList{RefId{10, "XMINYMIN"}, boundary_nodes});
+  }
+
+  {   // xmaxymin
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({cell_size[0], 0});
+    descriptor.addRefItemList(RefNodeList{RefId{11, "XMAXYMIN"}, boundary_nodes});
+  }
+
+  {   // xmaxymax
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({cell_size[0], cell_size[1]});
+    descriptor.addRefItemList(RefNodeList{RefId{12, "XMAXYMAX"}, boundary_nodes});
+  }
+
+  {   // xminymax
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({0, cell_size[1]});
+    descriptor.addRefItemList(RefNodeList{RefId{13, "XMINYMAX"}, boundary_nodes});
+  }
+}
+
+template <>
+void
+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};
+
+  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+    return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
+  };
+
+  {   // xminyminzmin
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({0, 0, 0});
+    descriptor.addRefItemList(RefNodeList{RefId{10, "XMINYMINZMIN"}, boundary_nodes});
+  }
+
+  {   // xmaxyminzmin
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({cell_size[0], 0, 0});
+    descriptor.addRefItemList(RefNodeList{RefId{11, "XMAXYMINZMIN"}, boundary_nodes});
+  }
+
+  {   // xmaxymaxzmin
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({cell_size[0], cell_size[1], 0});
+    descriptor.addRefItemList(RefNodeList{RefId{12, "XMAXYMAXZMIN"}, boundary_nodes});
+  }
+
+  {   // xminymaxzmin
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({0, cell_size[1], 0});
+    descriptor.addRefItemList(RefNodeList{RefId{13, "XMINYMAXZMIN"}, boundary_nodes});
+  }
+
+  {   // xminyminzmax
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({0, 0, cell_size[2]});
+    descriptor.addRefItemList(RefNodeList{RefId{14, "XMINYMINZMAX"}, boundary_nodes});
+  }
+
+  {   // xmaxyminzmax
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({cell_size[0], 0, cell_size[2]});
+    descriptor.addRefItemList(RefNodeList{RefId{15, "XMAXYMINZMAX"}, boundary_nodes});
+  }
+
+  {   // xmaxymaxzmax
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({cell_size[0], cell_size[1], cell_size[2]});
+    descriptor.addRefItemList(RefNodeList{RefId{16, "XMAXYMAXZMAX"}, boundary_nodes});
+  }
+
+  {   // xminymaxzmax
+    Array<NodeId> boundary_nodes(1);
+    boundary_nodes[0] = node_number({0, cell_size[1], cell_size[2]});
+    descriptor.addRefItemList(RefNodeList{RefId{17, "XMINYMAXZMAX"}, boundary_nodes});
+  }
+}
+
+template <size_t Dimension>
+void
+LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+{
+  static_assert(Dimension == 3, "unexpected dimension");
+}
+
+template <>
+void
+LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
+                                                   ConnectivityDescriptor& 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;
+  }();
+
+  const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
+  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+    return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
+  };
+
+  {   // Ridge edges in direction of Z axis
+    const auto add_ref_item_list_along_z = [&](const size_t i, const size_t j, const unsigned ref_id,
+                                               const std::string& ref_name) {
+      Array<EdgeId> boundary_edges(cell_size[2]);
+      size_t l = 0;
+      for (size_t k = 0; k < cell_size[2]; ++k) {
+        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+
+        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
+        Assert(i_edge != edge_to_id_map.end());
+
+        boundary_edges[l++] = i_edge->second;
+      }
+      Assert(l == cell_size[2]);
+      descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges});
+    };
+
+    add_ref_item_list_along_z(0, 0, 20, "XMINYMIN");
+    add_ref_item_list_along_z(0, cell_size[1], 21, "XMINYMAX");
+    add_ref_item_list_along_z(cell_size[0], cell_size[1], 22, "XMAXYMAX");
+    add_ref_item_list_along_z(cell_size[0], 0, 23, "XMAXYMIN");
+  }
+
+  {   // Ridge edges in direction of Y axis
+    const auto add_ref_item_list_along_y = [&](const size_t i, const size_t k, const unsigned ref_id,
+                                               const std::string& ref_name) {
+      Array<EdgeId> boundary_edges(cell_size[1]);
+      size_t l = 0;
+      for (size_t j = 0; j < cell_size[1]; ++j) {
+        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+
+        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
+        Assert(i_edge != edge_to_id_map.end());
+
+        boundary_edges[l++] = i_edge->second;
+      }
+      Assert(l == cell_size[1]);
+      descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges});
+    };
+
+    add_ref_item_list_along_y(0, 0, 24, "XMINZMIN");
+    add_ref_item_list_along_y(0, cell_size[2], 25, "XMINZMAX");
+    add_ref_item_list_along_y(cell_size[0], cell_size[2], 26, "XMAXZMAX");
+    add_ref_item_list_along_y(cell_size[0], 0, 27, "XMAXZMIN");
+  }
+
+  {   // Ridge edges in direction of X axis
+    const auto add_ref_item_list_along_x = [&](const size_t j, const size_t k, const unsigned ref_id,
+                                               const std::string& ref_name) {
+      Array<EdgeId> boundary_edges(cell_size[0]);
+      size_t l = 0;
+      for (size_t i = 0; i < cell_size[0]; ++i) {
+        const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+        const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
+
+        auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
+        Assert(i_edge != edge_to_id_map.end());
+
+        boundary_edges[l++] = i_edge->second;
+      }
+      Assert(l == cell_size[0]);
+      descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges});
+    };
+
+    add_ref_item_list_along_x(0, 0, 28, "YMINZMIN");
+    add_ref_item_list_along_x(0, cell_size[2], 29, "YMINZMAX");
+    add_ref_item_list_along_x(cell_size[1], cell_size[2], 30, "YMAXZMAX");
+    add_ref_item_list_along_x(cell_size[1], 0, 31, "YMAXZMIN");
+  }
+}
+
+template <size_t Dimension>
+void
+LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+{
+  static_assert(Dimension >= 2 and Dimension <= 3, "unexpected dimension");
+}
+
+template <>
+void
+LogicalConnectivityBuilder::_buildBoundaryFaceList(
+  const TinyVector<2, uint64_t>& cell_size,   // clazy:exclude=function-args-by-value
+  ConnectivityDescriptor& descriptor)
+{
+  const auto cell_number = [&](const TinyVector<2, uint64_t> cell_logic_id) {
+    return cell_logic_id[0] * cell_size[1] + cell_logic_id[1];
+  };
+
+  {   // xmin
+    const size_t i = 0;
+    Array<FaceId> boundary_faces(cell_size[1]);
+    for (size_t j = 0; j < cell_size[1]; ++j) {
+      constexpr size_t left_face = 3;
+
+      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+      const size_t face_id = descriptor.cell_to_face_vector[cell_id][left_face];
+
+      boundary_faces[j] = face_id;
+    }
+    descriptor.addRefItemList(RefFaceList{RefId{0, "XMIN"}, boundary_faces});
+  }
+
+  {   // xmax
+    const size_t i = cell_size[0] - 1;
+    Array<FaceId> boundary_faces(cell_size[1]);
+    for (size_t j = 0; j < cell_size[1]; ++j) {
+      constexpr size_t right_face = 1;
+
+      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+      const size_t face_id = descriptor.cell_to_face_vector[cell_id][right_face];
+
+      boundary_faces[j] = face_id;
+    }
+    descriptor.addRefItemList(RefFaceList{RefId{1, "XMAX"}, boundary_faces});
+  }
+
+  {   // ymin
+    const size_t j = 0;
+    Array<FaceId> boundary_faces(cell_size[0]);
+    for (size_t i = 0; i < cell_size[0]; ++i) {
+      constexpr size_t bottom_face = 0;
+
+      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+      const size_t face_id = descriptor.cell_to_face_vector[cell_id][bottom_face];
+
+      boundary_faces[i] = face_id;
+    }
+    descriptor.addRefItemList(RefFaceList{RefId{2, "YMIN"}, boundary_faces});
+  }
+
+  {   // ymax
+    const size_t j = cell_size[1] - 1;
+    Array<FaceId> boundary_faces(cell_size[0]);
+    for (size_t i = 0; i < cell_size[0]; ++i) {
+      constexpr size_t top_face = 2;
+
+      const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+      const size_t face_id = descriptor.cell_to_face_vector[cell_id][top_face];
+
+      boundary_faces[i] = face_id;
+    }
+    descriptor.addRefItemList(RefFaceList{RefId{3, "YMAX"}, boundary_faces});
+  }
+}
+
+template <>
+void
+LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
+                                                   ConnectivityDescriptor& descriptor)
+{
+  using Face = ConnectivityFace<3>;
+
+  const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
+    std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
+    for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
+      const auto& node_vector                                          = descriptor.face_to_node_vector[l];
+      face_to_id_map[Face(node_vector, descriptor.node_number_vector)] = l;
+    }
+    return face_to_id_map;
+  }();
+
+  const std::unordered_map<int, FaceId> face_number_id_map = [&] {
+    std::unordered_map<int, FaceId> face_number_id_map;
+    for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
+      face_number_id_map[descriptor.face_number_vector[l]] = l;
+    }
+    Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
+    return face_number_id_map;
+  }();
+
+  const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
+  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+    return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
+  };
+
+  {   // faces orthogonal to X
+    const auto add_ref_item_list_for_x = [&](const size_t i, const unsigned ref_id, const std::string& ref_name) {
+      Array<FaceId> boundary_faces(cell_size[1] * cell_size[2]);
+      size_t l = 0;
+      for (size_t j = 0; j < cell_size[1]; ++j) {
+        for (size_t k = 0; k < cell_size[2]; ++k) {
+          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k + 1});
+          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+
+          auto i_face =
+            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+          Assert(i_face != face_to_id_map.end());
+
+          boundary_faces[l++] = i_face->second;
+        }
+      }
+      Assert(l == cell_size[1] * cell_size[2]);
+      descriptor.addRefItemList(RefFaceList{RefId{ref_id, ref_name}, boundary_faces});
+    };
+
+    add_ref_item_list_for_x(0, 0, "XMIN");
+    add_ref_item_list_for_x(cell_size[0], 1, "XMAX");
+  }
+
+  {   // faces orthogonal to Y
+    const auto add_ref_item_list_for_y = [&](const size_t j, const unsigned ref_id, const std::string& ref_name) {
+      Array<FaceId> boundary_faces(cell_size[0] * cell_size[2]);
+      size_t l = 0;
+      for (size_t i = 0; i < cell_size[0]; ++i) {
+        for (size_t k = 0; k < cell_size[2]; ++k) {
+          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
+          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k + 1});
+          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+
+          auto i_face =
+            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+          Assert(i_face != face_to_id_map.end());
+
+          boundary_faces[l++] = i_face->second;
+        }
+      }
+      Assert(l == cell_size[0] * cell_size[2]);
+      descriptor.addRefItemList(RefFaceList{RefId{ref_id, ref_name}, boundary_faces});
+    };
+
+    add_ref_item_list_for_y(0, 2, "YMIN");
+    add_ref_item_list_for_y(cell_size[1], 3, "YMAX");
+  }
+
+  {   // faces orthogonal to Z
+    const auto add_ref_item_list_for_z = [&](const size_t k, const unsigned ref_id, const std::string& ref_name) {
+      Array<FaceId> boundary_faces(cell_size[0] * cell_size[1]);
+      size_t l = 0;
+      for (size_t i = 0; i < cell_size[0]; ++i) {
+        for (size_t j = 0; j < cell_size[1]; ++j) {
+          const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+          const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
+          const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k});
+          const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+
+          auto i_face =
+            face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+          Assert(i_face != face_to_id_map.end());
+
+          boundary_faces[l++] = i_face->second;
+        }
+      }
+      Assert(l == cell_size[0] * cell_size[1]);
+      descriptor.addRefItemList(RefFaceList{RefId{ref_id, ref_name}, boundary_faces});
+    };
+
+    add_ref_item_list_for_z(0, 4, "ZMIN");
+    add_ref_item_list_for_z(cell_size[2], 5, "ZMAX");
+  }
+}
+
+template <>
+void
+LogicalConnectivityBuilder::_buildConnectivity(
+  const TinyVector<1, uint64_t>& cell_size)   // clazy:exclude=function-args-by-value
+{
+  const size_t number_of_cells = cell_size[0];
+  const size_t number_of_nodes = cell_size[0] + 1;
+
+  ConnectivityDescriptor descriptor;
+  descriptor.node_number_vector.resize(number_of_nodes);
+  for (size_t i = 0; i < number_of_nodes; ++i) {
+    descriptor.node_number_vector[i] = i;
+  }
+
+  descriptor.cell_number_vector.resize(number_of_cells);
+  for (size_t i = 0; i < number_of_cells; ++i) {
+    descriptor.cell_number_vector[i] = i;
+  }
+
+  descriptor.cell_type_vector.resize(number_of_cells);
+  std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Line);
+
+  descriptor.cell_to_node_vector.resize(number_of_cells);
+  constexpr size_t nb_node_per_cell = 2;
+  for (size_t j = 0; j < number_of_cells; ++j) {
+    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
+    for (size_t r = 0; r < nb_node_per_cell; ++r) {
+      descriptor.cell_to_node_vector[j][0] = j;
+      descriptor.cell_to_node_vector[j][1] = j + 1;
+    }
+  }
+
+  this->_buildBoundaryNodeList(cell_size, descriptor);
+
+  descriptor.cell_owner_vector.resize(number_of_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 <>
+void
+LogicalConnectivityBuilder::_buildConnectivity(
+  const TinyVector<2, uint64_t>& cell_size)   // clazy:exclude=function-args-by-value
+{
+  constexpr size_t Dimension = 2;
+
+  const TinyVector<Dimension, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
+
+  const size_t number_of_cells = cell_size[0] * cell_size[1];
+  const size_t number_of_nodes = node_size[0] * node_size[1];
+
+  ConnectivityDescriptor descriptor;
+  descriptor.node_number_vector.resize(number_of_nodes);
+  for (size_t i = 0; i < number_of_nodes; ++i) {
+    descriptor.node_number_vector[i] = i;
+  }
+
+  descriptor.cell_number_vector.resize(number_of_cells);
+  for (size_t i = 0; i < number_of_cells; ++i) {
+    descriptor.cell_number_vector[i] = i;
+  }
+
+  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_number = [&](const TinyVector<Dimension, uint64_t> node_logic_id) {
+    return node_logic_id[0] * node_size[1] + node_logic_id[1];
+  };
+
+  const auto cell_logic_id = [&](size_t j) {
+    const uint64_t j0 = j / cell_size[1];
+    const uint64_t j1 = j % cell_size[1];
+    return TinyVector<Dimension, uint64_t>{j0, j1};
+  };
+
+  descriptor.cell_to_node_vector.resize(number_of_cells);
+  constexpr size_t nb_node_per_cell = 1 << Dimension;
+  for (size_t j = 0; j < number_of_cells; ++j) {
+    TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
+    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
+    for (size_t r = 0; r < nb_node_per_cell; ++r) {
+      descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0});
+      descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0});
+      descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1});
+      descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1});
+    }
+  }
+
+  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+
+  this->_buildBoundaryNodeList(cell_size, descriptor);
+  this->_buildBoundaryFaceList(cell_size, descriptor);
+
+  descriptor.cell_owner_vector.resize(number_of_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 = Connectivity<Dimension>::build(descriptor);
+}
+
+template <>
+void
+LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& cell_size)
+{
+  constexpr size_t Dimension = 3;
+
+  const TinyVector<Dimension, uint64_t> node_size = [&] {
+    TinyVector node_size{cell_size};
+    for (size_t i = 0; i < Dimension; ++i) {
+      node_size[i] += 1;
+    }
+    return node_size;
+  }();
+
+  auto count_items = [](auto&& v) {
+    using V_Type = std::decay_t<decltype(v)>;
+    static_assert(is_tiny_vector_v<V_Type>);
+
+    size_t sum = v[0];
+    for (size_t i = 1; i < Dimension; ++i) {
+      sum *= v[i];
+    }
+    return sum;
+  };
+
+  const size_t number_of_cells = count_items(cell_size);
+  const size_t number_of_nodes = count_items(node_size);
+
+  ConnectivityDescriptor descriptor;
+  descriptor.node_number_vector.resize(number_of_nodes);
+  for (size_t i = 0; i < number_of_nodes; ++i) {
+    descriptor.node_number_vector[i] = i;
+  }
+
+  descriptor.cell_number_vector.resize(number_of_cells);
+  for (size_t i = 0; i < number_of_cells; ++i) {
+    descriptor.cell_number_vector[i] = i;
+  }
+
+  descriptor.cell_type_vector.resize(number_of_cells);
+  std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Hexahedron);
+
+  const auto cell_logic_id = [&](size_t j) {
+    const size_t slice1  = cell_size[1] * cell_size[2];
+    const size_t& slice2 = cell_size[2];
+    const uint64_t j0    = j / slice1;
+    const uint64_t j1    = (j - j0 * slice1) / slice2;
+    const uint64_t j2    = j - (j0 * slice1 + j1 * slice2);
+    return TinyVector<Dimension, uint64_t>{j0, j1, j2};
+  };
+
+  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];
+  };
+
+  descriptor.cell_to_node_vector.resize(number_of_cells);
+  constexpr size_t nb_node_per_cell = 1 << Dimension;
+  for (size_t j = 0; j < number_of_cells; ++j) {
+    TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
+    descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
+    for (size_t r = 0; r < nb_node_per_cell; ++r) {
+      static_assert(Dimension == 3, "unexpected dimension");
+      descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 0});
+      descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 0});
+      descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 0});
+      descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 0});
+      descriptor.cell_to_node_vector[j][4] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 1});
+      descriptor.cell_to_node_vector[j][5] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 1});
+      descriptor.cell_to_node_vector[j][6] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 1});
+      descriptor.cell_to_node_vector[j][7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
+    }
+  }
+
+  ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+  ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
+
+  this->_buildBoundaryNodeList(cell_size, descriptor);
+  this->_buildBoundaryEdgeList(cell_size, descriptor);
+  this->_buildBoundaryFaceList(cell_size, descriptor);
+
+  descriptor.cell_owner_vector.resize(number_of_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 = Connectivity<Dimension>::build(descriptor);
+}
+
+template <size_t Dimension>
+LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size)
+{
+  if (parallel::rank() == 0) {
+    this->_buildConnectivity(size);
+  }
+}
+
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<1, uint64_t>&);
+
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<2, uint64_t>&);
+
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<3, uint64_t>&);
diff --git a/src/mesh/LogicalConnectivityBuilder.hpp b/src/mesh/LogicalConnectivityBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c63226dbfc35d97a575968a28cc7cee2bb5821ff
--- /dev/null
+++ b/src/mesh/LogicalConnectivityBuilder.hpp
@@ -0,0 +1,29 @@
+#ifndef LOGICAL_CONNECTIVITY_BUILDER_HPP
+#define LOGICAL_CONNECTIVITY_BUILDER_HPP
+
+#include <algebra/TinyVector.hpp>
+
+#include <mesh/ConnectivityBuilderBase.hpp>
+
+class LogicalConnectivityBuilder : public ConnectivityBuilderBase
+{
+ private:
+  template <size_t Dimension>
+  void _buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
+
+  template <size_t Dimension>
+  void _buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
+
+  template <size_t Dimension>
+  void _buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
+
+  template <size_t Dimension>
+  void _buildConnectivity(const TinyVector<Dimension, uint64_t>& size);
+
+ public:
+  template <size_t Dimension>
+  LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size);
+  ~LogicalConnectivityBuilder() = default;
+};
+
+#endif   // LOGICAL_CONNECTIVITY_BUILDER_HPP
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index fc9020325fe64cb2d27dd3feab97ffcf778cdb3c..9e62baf56e7915379d65a5c88c21fe7411fca253 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -41,410 +41,3 @@ MeshBuilderBase::_dispatch()
 template void MeshBuilderBase::_dispatch<1>();
 template void MeshBuilderBase::_dispatch<2>();
 template void MeshBuilderBase::_dispatch<3>();
-
-/// ============
-
-template <size_t Dimension>
-void
-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>;
-  using Face         = ConnectivityFace<Dimension>;
-
-  const auto& node_number_vector = descriptor.node_number_vector;
-  Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
-  std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
-  for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
-    const auto& cell_nodes = descriptor.cell_to_node_vector[j];
-
-    if constexpr (Dimension == 2) {
-      switch (descriptor.cell_type_vector[j]) {
-      case CellType::Triangle: {
-        cell_nb_faces[j] = 3;
-        // face 0
-        Face f0({cell_nodes[1], cell_nodes[2]}, node_number_vector);
-        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
-        // face 1
-        Face f1({cell_nodes[2], cell_nodes[0]}, node_number_vector);
-        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
-        // face 2
-        Face f2({cell_nodes[0], cell_nodes[1]}, node_number_vector);
-        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-        break;
-      }
-      case CellType::Quadrangle: {
-        cell_nb_faces[j] = 4;
-        // face 0
-        Face f0({cell_nodes[0], cell_nodes[1]}, node_number_vector);
-        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
-        // face 1
-        Face f1({cell_nodes[1], cell_nodes[2]}, node_number_vector);
-        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
-        // face 2
-        Face f2({cell_nodes[2], cell_nodes[3]}, node_number_vector);
-        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
-        // face 3
-        Face f3({cell_nodes[3], cell_nodes[0]}, node_number_vector);
-        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
-        break;
-      }
-      default: {
-        std::ostringstream error_msg;
-        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
-        throw UnexpectedError(error_msg.str());
-      }
-      }
-    } else if constexpr (Dimension == 3) {
-      switch (descriptor.cell_type_vector[j]) {
-      case CellType::Hexahedron: {
-        // face 0
-        Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
-        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
-        // face 1
-        Face f1({cell_nodes[4], cell_nodes[5], cell_nodes[6], cell_nodes[7]}, node_number_vector);
-        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
-        // face 2
-        Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
-        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
-        // face 3
-        Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
-        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
-
-        // face 4
-        Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
-        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
-
-        // face 5
-        Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
-        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
-
-        cell_nb_faces[j] = 6;
-        break;
-      }
-      case CellType::Tetrahedron: {
-        cell_nb_faces[j] = 4;
-        // face 0
-        Face f0({cell_nodes[1], cell_nodes[2], cell_nodes[3]}, node_number_vector);
-        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
-        // face 1
-        Face f1({cell_nodes[0], cell_nodes[3], cell_nodes[2]}, node_number_vector);
-        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
-        // face 2
-        Face f2({cell_nodes[0], cell_nodes[1], cell_nodes[3]}, node_number_vector);
-        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
-        // face 3
-        Face f3({cell_nodes[0], cell_nodes[2], cell_nodes[1]}, node_number_vector);
-        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
-        break;
-      }
-      case CellType::Pyramid: {
-        cell_nb_faces[j] = cell_nodes.size();
-        std::vector<unsigned int> base_nodes;
-        std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
-
-        // base face
-        {
-          Face base_face(base_nodes, node_number_vector);
-          face_cells_map[base_face].emplace_back(std::make_tuple(j, 0, base_face.reversed()));
-        }
-        // side faces
-        const auto pyramid_vertex = cell_nodes[cell_nodes.size() - 1];
-        for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
-          Face side_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], pyramid_vertex},
-                         node_number_vector);
-          face_cells_map[side_face].emplace_back(std::make_tuple(j, i_node + 1, side_face.reversed()));
-        }
-        break;
-      }
-      case CellType::Diamond: {
-        cell_nb_faces[j] = 2 * (cell_nodes.size() - 2);
-        std::vector<unsigned int> base_nodes;
-        std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
-
-        {   // top faces
-          const auto top_vertex = cell_nodes[cell_nodes.size() - 1];
-          for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
-            Face top_face({base_nodes[i_node], base_nodes[(i_node + 1) % base_nodes.size()], top_vertex},
-                          node_number_vector);
-            face_cells_map[top_face].emplace_back(std::make_tuple(j, i_node, top_face.reversed()));
-          }
-        }
-
-        {   // bottom faces
-          const auto bottom_vertex = cell_nodes[0];
-          for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
-            Face bottom_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], bottom_vertex},
-                             node_number_vector);
-            face_cells_map[bottom_face].emplace_back(
-              std::make_tuple(j, i_node + base_nodes.size(), bottom_face.reversed()));
-          }
-        }
-        break;
-      }
-      default: {
-        std::ostringstream error_msg;
-        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
-        throw UnexpectedError(error_msg.str());
-      }
-      }
-    }
-  }
-
-  {
-    descriptor.cell_to_face_vector.resize(descriptor.cell_to_node_vector.size());
-    for (CellId j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
-      descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
-    }
-    FaceId l = 0;
-    for (const auto& face_cells_vector : face_cells_map) {
-      const auto& cells_vector = face_cells_vector.second;
-      for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
-        const auto& [cell_number, cell_local_face, reversed]         = cells_vector[lj];
-        descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
-      }
-      ++l;
-    }
-  }
-
-  {
-    descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
-    for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
-      descriptor.cell_face_is_reversed_vector[j] = Array<bool>(cell_nb_faces[j]);
-    }
-    for (const auto& face_cells_vector : face_cells_map) {
-      const auto& cells_vector = face_cells_vector.second;
-      for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
-        const auto& [cell_number, cell_local_face, reversed]                  = cells_vector[lj];
-        descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
-      }
-    }
-  }
-
-  {
-    descriptor.face_to_node_vector.resize(face_cells_map.size());
-    int l = 0;
-    for (const auto& face_info : face_cells_map) {
-      const Face& face                  = face_info.first;
-      descriptor.face_to_node_vector[l] = face.nodeIdList();
-      ++l;
-    }
-  }
-
-  {
-    // Face numbers may change if numbers are provided in the file
-    descriptor.face_number_vector.resize(face_cells_map.size());
-    for (size_t l = 0; l < face_cells_map.size(); ++l) {
-      descriptor.face_number_vector[l] = l;
-    }
-  }
-}
-
-template <size_t Dimension>
-void
-ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
-{
-  static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
-  using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
-  using Edge         = ConnectivityFace<2>;
-
-  const auto& node_number_vector = descriptor.node_number_vector;
-  Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
-  std::map<Edge, std::vector<FaceEdgeInfo>> edge_faces_map;
-  for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-    const auto& face_nodes = descriptor.face_to_node_vector[l];
-
-    face_nb_edges[l] = face_nodes.size();
-    for (size_t r = 0; r < face_nodes.size() - 1; ++r) {
-      Edge e({face_nodes[r], face_nodes[r + 1]}, node_number_vector);
-      edge_faces_map[e].emplace_back(std::make_tuple(l, r, e.reversed()));
-    }
-    {
-      Edge e({face_nodes[face_nodes.size() - 1], face_nodes[0]}, node_number_vector);
-      edge_faces_map[e].emplace_back(std::make_tuple(l, face_nodes.size() - 1, e.reversed()));
-    }
-  }
-
-  std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_id_map;
-  {
-    descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
-    for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
-      descriptor.face_to_edge_vector[l].resize(face_nb_edges[l]);
-    }
-    EdgeId e = 0;
-    for (const auto& edge_faces_vector : edge_faces_map) {
-      const auto& faces_vector = edge_faces_vector.second;
-      for (unsigned short l = 0; l < faces_vector.size(); ++l) {
-        const auto& [face_number, face_local_edge, reversed]         = faces_vector[l];
-        descriptor.face_to_edge_vector[face_number][face_local_edge] = e;
-      }
-      edge_id_map[edge_faces_vector.first] = e;
-      ++e;
-    }
-  }
-
-  {
-    descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
-    for (FaceId j = 0; j < descriptor.face_edge_is_reversed_vector.size(); ++j) {
-      descriptor.face_edge_is_reversed_vector[j] = Array<bool>(face_nb_edges[j]);
-    }
-    for (const auto& edge_faces_vector : edge_faces_map) {
-      const auto& faces_vector = edge_faces_vector.second;
-      for (unsigned short lj = 0; lj < faces_vector.size(); ++lj) {
-        const auto& [face_number, face_local_edge, reversed]                  = faces_vector[lj];
-        descriptor.face_edge_is_reversed_vector[face_number][face_local_edge] = reversed;
-      }
-    }
-  }
-
-  {
-    descriptor.edge_to_node_vector.resize(edge_faces_map.size());
-    int e = 0;
-    for (const auto& edge_info : edge_faces_map) {
-      const Edge& edge                  = edge_info.first;
-      descriptor.edge_to_node_vector[e] = edge.nodeIdList();
-      ++e;
-    }
-  }
-
-  {
-    // Edge numbers may change if numbers are provided in the file
-    descriptor.edge_number_vector.resize(edge_faces_map.size());
-    for (size_t e = 0; e < edge_faces_map.size(); ++e) {
-      descriptor.edge_number_vector[e] = e;
-    }
-  }
-
-  {
-    descriptor.cell_to_edge_vector.reserve(descriptor.cell_to_node_vector.size());
-    for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
-      const auto& cell_nodes = descriptor.cell_to_node_vector[j];
-
-      switch (descriptor.cell_type_vector[j]) {
-      case CellType::Tetrahedron: {
-        constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(6);
-        for (int i_edge = 0; i_edge < 6; ++i_edge) {
-          const auto e = local_edge[i_edge];
-          Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
-          auto i = edge_id_map.find(edge);
-          if (i == edge_id_map.end()) {
-            throw NormalError("could not find this edge");
-          }
-          cell_edge_vector.push_back(i->second);
-        }
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
-      }
-      case CellType::Hexahedron: {
-        constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
-                                           {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(12);
-        for (int i_edge = 0; i_edge < 12; ++i_edge) {
-          const auto e = local_edge[i_edge];
-          Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
-          auto i = edge_id_map.find(edge);
-          if (i == edge_id_map.end()) {
-            throw NormalError("could not find this edge");
-          }
-          cell_edge_vector.push_back(i->second);
-        }
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
-      }
-      case CellType::Pyramid: {
-        const size_t number_of_edges = 2 * cell_nodes.size();
-        std::vector<unsigned int> base_nodes;
-        std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
-
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(number_of_edges);
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
-          auto i = edge_id_map.find(edge);
-          if (i == edge_id_map.end()) {
-            throw NormalError("could not find this edge");
-          }
-          cell_edge_vector.push_back(i->second);
-        }
-
-        const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
-          auto i = edge_id_map.find(edge);
-          if (i == edge_id_map.end()) {
-            throw NormalError("could not find this edge");
-          }
-          cell_edge_vector.push_back(i->second);
-        }
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-        break;
-      }
-      case CellType::Diamond: {
-        const size_t number_of_edges = 3 * cell_nodes.size();
-        std::vector<unsigned int> base_nodes;
-        std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
-
-        std::vector<unsigned int> cell_edge_vector;
-        cell_edge_vector.reserve(number_of_edges);
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
-          auto i = edge_id_map.find(edge);
-          if (i == edge_id_map.end()) {
-            throw NormalError("could not find this edge");
-          }
-          cell_edge_vector.push_back(i->second);
-        }
-
-        const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
-          auto i = edge_id_map.find(edge);
-          if (i == edge_id_map.end()) {
-            throw NormalError("could not find this edge");
-          }
-          cell_edge_vector.push_back(i->second);
-        }
-
-        const unsigned int bottom_vertex = cell_nodes[0];
-        for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
-          Edge edge{{base_nodes[i_edge], bottom_vertex}, node_number_vector};
-          auto i = edge_id_map.find(edge);
-          if (i == edge_id_map.end()) {
-            throw NormalError("could not find this edge");
-          }
-          cell_edge_vector.push_back(i->second);
-        }
-
-        descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-
-        break;
-      }
-      default: {
-        std::stringstream error_msg;
-        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
-        throw NotImplementedError(error_msg.str());
-      }
-      }
-    }
-  }
-}
-
-template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
-template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
-
-template void ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
-  ConnectivityDescriptor& descriptor);
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
index 24bd055824223a604d1a913d1e85af0a757b1965..6c0328bb8ee624e19d42fb4660ba3665ff7164de 100644
--- a/src/mesh/MeshBuilderBase.hpp
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -3,227 +3,7 @@
 
 #include <mesh/IMesh.hpp>
 
-#include <mesh/IConnectivity.hpp>
-
-#include <mesh/ItemId.hpp>
-#include <utils/PugsAssert.hpp>
-#include <utils/PugsMacros.hpp>
-
 #include <memory>
-#include <vector>
-
-class ConnectivityDescriptor;
-
-class ConnectivityBuilderBase
-{
- protected:
-  template <size_t Dimension>
-  class ConnectivityFace;
-
-  std::shared_ptr<const IConnectivity> m_connectivity;
-
-  template <size_t Dimension>
-  static void _computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
-
-  template <size_t Dimension>
-  static void _computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
-
- public:
-  std::shared_ptr<const IConnectivity>
-  connectivity() const
-  {
-    return m_connectivity;
-  }
-
-  ConnectivityBuilderBase()  = default;
-  ~ConnectivityBuilderBase() = default;
-};
-
-template <>
-class ConnectivityBuilderBase::ConnectivityFace<2>
-{
- public:
-  friend struct Hash;
-
-  struct Hash
-  {
-    size_t
-    operator()(const ConnectivityFace& f) const
-    {
-      size_t hash = 0;
-      hash ^= std::hash<unsigned int>()(f.m_node0_id);
-      hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
-      return hash;
-    }
-  };
-
- private:
-  const std::vector<int>& m_node_number_vector;
-
-  unsigned int m_node0_id;
-  unsigned int m_node1_id;
-
-  bool m_reversed;
-
- public:
-  std::vector<unsigned int>
-  nodeIdList() const
-  {
-    return {m_node0_id, m_node1_id};
-  }
-
-  bool
-  reversed() const
-  {
-    return m_reversed;
-  }
-
-  PUGS_INLINE
-  bool
-  operator==(const ConnectivityFace& f) const
-  {
-    return ((m_node0_id == f.m_node0_id) and (m_node1_id == f.m_node1_id));
-  }
-
-  PUGS_INLINE
-  bool
-  operator<(const ConnectivityFace& f) const
-  {
-    return ((m_node_number_vector[m_node0_id] < m_node_number_vector[f.m_node0_id]) or
-            ((m_node_number_vector[m_node0_id] == m_node_number_vector[f.m_node0_id]) and
-             (m_node_number_vector[m_node1_id] < m_node_number_vector[f.m_node1_id])));
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const std::vector<unsigned int>& node_id_list, const std::vector<int>& node_number_vector)
-    : m_node_number_vector(node_number_vector)
-  {
-    Assert(node_id_list.size() == 2);
-
-    if (m_node_number_vector[node_id_list[0]] < m_node_number_vector[node_id_list[1]]) {
-      m_node0_id = node_id_list[0];
-      m_node1_id = node_id_list[1];
-      m_reversed = false;
-    } else {
-      m_node0_id = node_id_list[1];
-      m_node1_id = node_id_list[0];
-      m_reversed = true;
-    }
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const ConnectivityFace&) = default;
-
-  PUGS_INLINE
-  ~ConnectivityFace() = default;
-};
-
-template <>
-class ConnectivityBuilderBase::ConnectivityFace<3>
-{
- public:
-  friend struct Hash;
-
-  struct Hash
-  {
-    size_t
-    operator()(const ConnectivityFace& f) const
-    {
-      size_t hash = 0;
-      for (size_t i = 0; i < f.m_node_id_list.size(); ++i) {
-        hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
-      }
-      return hash;
-    }
-  };
-
- private:
-  bool m_reversed;
-  std::vector<NodeId::base_type> m_node_id_list;
-  const std::vector<int>& m_node_number_vector;
-
-  PUGS_INLINE
-  std::vector<unsigned int>
-  _sort(const std::vector<unsigned int>& node_list)
-  {
-    const auto min_id = std::min_element(node_list.begin(), node_list.end());
-    const int shift   = std::distance(node_list.begin(), min_id);
-
-    std::vector<unsigned int> rotated_node_list(node_list.size());
-    if (node_list[(shift + 1) % node_list.size()] > node_list[(shift + node_list.size() - 1) % node_list.size()]) {
-      for (size_t i = 0; i < node_list.size(); ++i) {
-        rotated_node_list[i] = node_list[(shift + node_list.size() - i) % node_list.size()];
-        m_reversed           = true;
-      }
-    } else {
-      for (size_t i = 0; i < node_list.size(); ++i) {
-        rotated_node_list[i] = node_list[(shift + i) % node_list.size()];
-      }
-    }
-
-    return rotated_node_list;
-  }
-
- public:
-  PUGS_INLINE
-  const bool&
-  reversed() const
-  {
-    return m_reversed;
-  }
-
-  PUGS_INLINE
-  const std::vector<unsigned int>&
-  nodeIdList() const
-  {
-    return m_node_id_list;
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list, const std::vector<int>& node_number_vector)
-    : m_reversed(false), m_node_id_list(_sort(given_node_id_list)), m_node_number_vector(node_number_vector)
-  {
-    ;
-  }
-
- public:
-  bool
-  operator==(const ConnectivityFace& f) const
-  {
-    if (m_node_id_list.size() == f.nodeIdList().size()) {
-      for (size_t j = 0; j < m_node_id_list.size(); ++j) {
-        if (m_node_id_list[j] != f.nodeIdList()[j]) {
-          return false;
-        }
-      }
-      return true;
-    }
-    return false;
-  }
-
-  PUGS_INLINE
-  bool
-  operator<(const ConnectivityFace& f) const
-  {
-    const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
-    for (size_t i = 0; i < min_nb_nodes; ++i) {
-      if (m_node_id_list[i] < f.m_node_id_list[i])
-        return true;
-      if (m_node_id_list[i] != f.m_node_id_list[i])
-        return false;
-    }
-    return m_node_id_list.size() < f.m_node_id_list.size();
-  }
-
-  PUGS_INLINE
-  ConnectivityFace(const ConnectivityFace&) = default;
-
-  PUGS_INLINE
-  ConnectivityFace() = delete;
-
-  PUGS_INLINE
-  ~ConnectivityFace() = default;
-};
 
 class MeshBuilderBase
 {