diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0e2991963a15afe8e4d27b5f6f0388f51e03d650..d3306db95476b2e91773271f548b3937a417b421 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -145,7 +145,7 @@ endif()
 # Search for ParMETIS
 
 find_package(ParMETIS)
-if(PARMETIS_LIBRARIES)
+if(PARMETIS_LIBRARIES AND MPI_FOUND)
   add_library(PkgConfig::ParMETIS STATIC IMPORTED)
   set_property(TARGET PkgConfig::ParMETIS PROPERTY
     IMPORTED_LOCATION "${PARMETIS_LIBRARIES}")
@@ -161,9 +161,7 @@ endif()
 
 find_package(PTScotch)
 
-set(PUGS_HAS_PTSCOTCH ${PTScotch_FOUND})
-
-if (PTScotch_FOUND)
+if (PTScotch_FOUND AND MPI_FOUND)
   add_library(PkgConfig::PTScotch STATIC IMPORTED)
   set_property(TARGET PkgConfig::PTScotch PROPERTY
     IMPORTED_LOCATION "${PTSCOTCH_LIBRARY}")
@@ -172,8 +170,10 @@ if (PTScotch_FOUND)
 
   set(PTSCOTCH_TARGET PkgConfig::PTScotch)
   include_directories(SYSTEM "${PTSCOTCH_INCLUDE_DIR}")
+  set(PUGS_HAS_PTSCOTCH TRUE)
 else()
   set(PTSCOTCH_LIBRARY "")
+  unset(PUGS_HAS_PTSCOTCH)
 endif()
 
 #------------------------------------------------------
@@ -852,9 +852,9 @@ if (MPI_FOUND)
   message(" MPI: ${MPI_CXX_LIBRARY_VERSION_STRING}")
 else()
   if (PUGS_ENABLE_MPI MATCHES "^OFF$")
-    message(" MPI: deactivated!")
-  elseif(NOT PARMETIS_LIBRARIES)
-    message(" MPI: deactivated: ParMETIS cannot be found!")
+    message(" MPI: explicitly deactivated!")
+  elseif((NOT DEFINED PUGS_HAS_PARMETIS) AND (NOT DEFINED PUGS_HAS_PTSCOTCH))
+    message(" MPI: deactivated: ParMETIS and PTScotch cannot be found!")
   else()
     if (PUGS_ENABLE_MPI MATCHES "^(AUTO|ON)$")
       message(" MPI: not found!")
@@ -864,6 +864,26 @@ else()
   endif()
 endif()
 
+if (PUGS_HAS_PARMETIS)
+  message(" ParMETIS: ${PARMETIS_LIBRARIES}")
+else()
+  if (PUGS_HAS_MPI)
+  message(" ParMETIS: not found!")
+  else()
+  message(" ParMETIS: deactivated (MPI not found)")
+  endif()
+endif()
+
+if (PUGS_HAS_PTSCOTCH)
+  message(" PTScotch: ${PTSCOTCH_LIBRARIES}")
+else()
+  if (PUGS_HAS_MPI)
+  message(" PTScotch: not found!")
+  else()
+  message(" PTScotch: deactivated (MPI not found)")
+  endif()
+endif()
+
 if (PETSC_FOUND)
   message(" PETSc: ${PETSC_VERSION}")
 else()
diff --git a/README.md b/README.md
index 0fb43ab24b26b2aab63744031199dbd4553c6fca..059c2960300263220474050aef84f1e81c68f961 100644
--- a/README.md
+++ b/README.md
@@ -54,12 +54,12 @@ functionalities or development tools.
 
 `pugs` integrates MPI support and should compile with any MPI-3
 implementation. However in order to dispatch meshes on MPI processes
-`pugs` requires `ParMETIS`. Without, MPI support will be automatically
-deactivated.
+`pugs` requires `ParMETIS` or `PTScotch`. Without one of them, the MPI
+support will be automatically deactivated.
 
 The easiest way to enable MPI support on Debian-like systems is to run
 ```shell
-apt install libparmetis-dev
+apt install libparmetis-dev libscotch-dev
 ```
 
 #### `PETSc`
diff --git a/doc/userdoc.org b/doc/userdoc.org
index 243fbde6af57b702f741d211bc9cb734ebb73364..cca6363ed3a545ea0515d43f03ee716810e8ef3d 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -2689,9 +2689,9 @@ actually stored in the ~checkpoint.h5~ file, the script is not a command
 line argument.
 
 By now, it is not possible to *simply* modify the script while
-resuming. This is done in purpose since it remains unclear if such a
-dangerous functionality should be make easy. Indeed allowing
-script modifications may lead to undefined behaviors.
+resuming. This is done on purpose since it remains unclear if such a
+dangerous functionality should be make easy. Indeed allowing script
+modifications may lead to undefined behaviors.
 #+END_note
 
 #+BEGIN_warning
@@ -3431,6 +3431,30 @@ In this example, we set three arrays defined at all nodes, all the
 - ~face_values~ interpolates ~f~ at the centers of the faces of ~m~,
 - ~cell_values~ interpolates ~f~ at the centers of the cells of ~m~.
 
+***** ~load_balance:mesh -> mesh~
+
+Creates a new partition of a mesh in parallel (~MPI~).
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+
+  let m1:mesh, m1 = cartesianMesh([0,0], [1,1], (10,10));
+  let m2:mesh, m2 = load_balance(m1);
+#+END_SRC
+In this example the mesh ~m2~ is a new parallel partition of the mesh
+~m1~. This implies that the two meshes, while being identical from the
+mathematical point of view, are now different (and are not defined on
+a same connectivity: the connectivity is actually balanced).
+
+Discrete functions defined on ~m1~ are *not* implicitly transferred on the
+mesh ~m2~! The function ~load_balance:(Vh)->(Vh)~ must be used to perform
+mesh and data load balancing.
+
+#+BEGIN_note
+To simplify development of large calculations. A sequential call of
+this function has the same effect: the mesh and its connectivity are
+copied!
+#+END_note
+
 ***** ~transform: mesh*function -> mesh~
 
 This function allows to compute a new mesh as the transformation of a
@@ -4471,6 +4495,109 @@ to different meshes.
 
 A good example of the use of this kind of function is mass fractions.
 
+***** ~get_mesh:Vh -> mesh~
+
+This function allows to retrieve the mesh that is used to define a
+discrete function.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+
+  let m1:mesh, m1 = cartesianMesh(0, [1,1], (10,10));
+
+  let f: R^2 -> R, x -> 2*x[0]+x[1];
+  let fh:Vh, fh = interpolate(m1, P0(), f);
+
+  let m2: mesh, m2 = get_mesh(fh);
+#+END_SRC
+Observe that in this example the two variables ~m1~ and ~m2~ refer to the
+same exact mesh. There is no duplication.
+
+***** ~load_balance: (Vh) -> (Vh)~
+
+This function performs a parallel load balancing of a list of ~Vh~
+variables. All the input variables must be defined on a same mesh. The
+return list consists of the balanced variables sorted using the input
+ordering. The new (balanced) mesh is carried out by the new variables
+and can be retrieved using the ~get_mesh: Vh -> mesh~ function.
+
+Here is an example of use.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+
+  let m:mesh, m = cartesianMesh(0, [1,1], (10,10));
+
+  let f: R^2 -> R, x -> 2*x[0]+x[1];
+  let g: R^2 -> R, x -> 3*x[1]-x[0];
+
+  let fh:Vh, fh = interpolate(m, P0(), f);
+  let gh:Vh, gh = interpolate(m, P0(), g);
+
+  (fh,gh) = load_balance((fh, gh));
+  let sum:Vh, sum = fh+gh;
+
+  m = get_mesh(sum);
+#+END_SRC
+
+#+BEGIN_note
+One notices the parenthesis enclosing the list ~(fh,gh)~ when invoking
+the ~load_balance~ function.
+#+END_note
+
+One should also observe that the in the previous example, balancing ~fh~
+and ~gh~ separately defines two new different meshes so that the resulting
+discrete functions are no more defined on the same one. The following
+example enlightens this behavior.
+
+#+NAME: separate-load-balance
+#+BEGIN_SRC pugs-error :exports both :results output
+  import mesh;
+  import scheme;
+
+  let m:mesh, m = cartesianMesh(0, [1,1], (10,10));
+
+  let f: R^2 -> R, x -> 2*x[0]+x[1];
+  let g: R^2 -> R, x -> 3*x[1]-x[0];
+
+  let fh:Vh, fh = interpolate(m, P0(), f);
+  let gh:Vh, gh = interpolate(m, P0(), g);
+
+  fh = load_balance(fh);
+  gh = load_balance(gh);
+  let sum:Vh, sum = fh+gh;
+#+END_SRC
+#+RESULTS: separate-load-balance
+
+#+BEGIN_warning
+Using this function in a sequential calculation has the same behavior:
+the discrete functions, the mesh and its connectivity are copied. This
+is intended to simplify development of large calculations on simpler
+tests, and for consistency.
+#+END_warning
+
+The following example illustrates this remark: even for sequential
+calculations a new mesh (and connectivity) is built *on purpose* when
+invoking ~load_balance~
+#+NAME: serial-effect-load-balance
+#+BEGIN_SRC pugs-error :exports both :results output
+  import mesh;
+  import scheme;
+
+  let m:mesh, m = cartesianMesh(0, [1,1], (10,10));
+
+  let f: R^2 -> R, x -> 2*x[0]+x[1];
+
+  let fh:Vh, fh = interpolate(m, P0(), f);
+  let gh:Vh, gh = load_balance(fh);
+
+  let sum:Vh, sum = fh+gh;
+#+END_SRC
+#+RESULTS: serial-effect-load-balance
+Even if the functions are similar (mathematically ~fh~ = ~gh~) they do not
+live on the same mesh from ~pugs~ point of view, thus one cannot perform
+simple algebra.
+
 ***** Numerical methods
 
 We describe rapidly two functions that embed numerical methods. These
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 9b698b9e095ff6c227c58a2e6fe6bae5f7f72ad4..b1faa718bb00272dcf546c91b9021449a04c40a6 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -24,6 +24,7 @@
 #include <mesh/ItemArrayVariant.hpp>
 #include <mesh/ItemValueVariant.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshBalancer.hpp>
 #include <mesh/MeshRelaxer.hpp>
 #include <mesh/MeshTransformer.hpp>
 #include <mesh/MeshUtils.hpp>
@@ -349,6 +350,16 @@ MeshModule::MeshModule()
                                             }
 
                                             ));
+
+  this->_addBuiltinFunction("load_balance", std::function(
+
+                                              [](const std::shared_ptr<const MeshVariant>& mesh_v)
+                                                -> std::shared_ptr<const MeshVariant> {
+                                                MeshBalancer mesh_balancer(mesh_v);
+                                                return mesh_balancer.mesh();
+                                              }
+
+                                              ));
 }
 
 void
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 40cc703cc6dcc0029da81366fb95134fc1329599..b9b789cfa8e3d7bc5328062c6085217dbe142120 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -46,6 +46,7 @@
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/InflowBoundaryConditionDescriptor.hpp>
 #include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/LoadBalancer.hpp>
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -731,6 +732,30 @@ SchemeModule::SchemeModule()
 
                                                    ));
 
+  this->_addBuiltinFunction("get_mesh",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_variant)
+                                -> std::shared_ptr<const MeshVariant> {
+                                return std::visit(   //
+                                  [](auto&& discrete_function) {
+                                    return discrete_function.meshVariant();   //
+                                  },
+                                  discrete_function_variant->discreteFunction());
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("load_balance", std::function(
+
+                                              [](const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
+                                                   discrete_function_variant_list)
+                                                -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
+                                                return LoadBalancer{}.balance(discrete_function_variant_list);
+                                              }
+
+                                              ));
+
   MathFunctionRegisterForVh{*this};
 }
 
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 76f5f96d0b98944d3bcd65e42e36a705541886f9..bd45f495eaf08c8b9e000649a4500ddb146aac95 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -20,6 +20,7 @@ add_library(
   MedianDualConnectivityBuilder.cpp
   MedianDualMeshBuilder.cpp
   Mesh.cpp
+  MeshBalancer.cpp
   MeshBuilderBase.cpp
   MeshCellZone.cpp
   MeshData.cpp
diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
index 4ae7410ce7b896e7cf1f8147a0b64365cebc9352..d0b34bb128b07d5c47ac1e958ab3bf20dafb9cc1 100644
--- a/src/mesh/CartesianMeshBuilder.cpp
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -144,7 +144,7 @@ CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<Dimension>& a,
 
     this->_buildCartesianMesh(corner0, corner1, size);
   }
-  this->_dispatch<Dimension>();
+  this->_dispatch<Dimension>(DispatchType::initial);
 }
 
 template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<1>&,
diff --git a/src/mesh/CartesianMeshBuilder.hpp b/src/mesh/CartesianMeshBuilder.hpp
index 7b32b55b64d021c655ff9c35245ece4b31b2ed76..c798ae6233e2aa4b699568e9ba58c48075409ad3 100644
--- a/src/mesh/CartesianMeshBuilder.hpp
+++ b/src/mesh/CartesianMeshBuilder.hpp
@@ -5,8 +5,6 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshBuilderBase.hpp>
 
-#include <memory>
-
 class CartesianMeshBuilder : public MeshBuilderBase
 {
  private:
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 65a8804950abad22fe2205569217cd8ff162b043..44c07b62877451940e975fbc929b6da1530a36f4 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -46,14 +46,6 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
   m_cell_number = WeakCellValue<const int>(*this, descriptor.cellNumberVector());
   m_node_number = WeakNodeValue<const int>(*this, descriptor.nodeNumberVector());
 
-  {
-    WeakCellValue<int> cell_global_index(*this);
-    int first_index = 0;
-    parallel_for(
-      this->numberOfCells(), PUGS_LAMBDA(CellId j) { cell_global_index[j] = first_index + j; });
-    m_cell_global_index = cell_global_index;
-  }
-
   m_cell_owner = WeakCellValue<const int>(*this, descriptor.cellOwnerVector());
 
   {
@@ -310,6 +302,80 @@ Connectivity<Dimension>::_write(std::ostream& os) const
   return os;
 }
 
+template <size_t Dimension>
+CRSGraph
+Connectivity<Dimension>::ownCellToCellGraph() const
+{
+  std::vector<std::vector<int>> cell_cells(this->numberOfCells());
+  const auto& face_to_cell_matrix = this->faceToCellMatrix();
+
+  for (FaceId l = 0; l < this->numberOfFaces(); ++l) {
+    const auto& face_to_cell = face_to_cell_matrix[l];
+    if (face_to_cell.size() > 1) {
+      const CellId cell_0 = face_to_cell[0];
+      const CellId cell_1 = face_to_cell[1];
+
+      if (m_cell_is_owned[cell_0]) {
+        cell_cells[cell_0].push_back(cell_1);
+      }
+      if (m_cell_is_owned[cell_1]) {
+        cell_cells[cell_1].push_back(cell_0);
+      }
+    }
+  }
+
+  std::vector<std::vector<int>> owned_cell_cells;
+  owned_cell_cells.reserve(this->numberOfCells());
+
+  size_t number_of_owned_cells = 0;
+  for (CellId cell_id = 0; cell_id < this->numberOfCells(); ++cell_id) {
+    if (m_cell_is_owned[cell_id]) {
+      ++number_of_owned_cells;
+    }
+  }
+
+  for (size_t i_cell = 0; i_cell < this->numberOfCells(); ++i_cell) {
+    if (cell_cells[i_cell].size() > 0) {
+      owned_cell_cells.emplace_back(std::move(cell_cells[i_cell]));
+    }
+  }
+
+  Array number_of_owned_cells_by_rank = parallel::allGather(number_of_owned_cells);
+
+  // index is defined locally. First we compute the first index on the
+  // processor
+  size_t global_index = 0;
+  for (size_t i = 0; i < parallel::rank(); ++i) {
+    global_index += number_of_owned_cells_by_rank[i];
+  }
+
+  CellValue<int> global_cell_index{*this};
+  for (CellId cell_id = 0; cell_id < this->numberOfCells(); ++cell_id) {
+    if (m_cell_is_owned[cell_id]) {
+      global_cell_index[cell_id] = global_index++;
+    }
+  }
+  synchronize(global_cell_index);
+
+  Array<int> entries(owned_cell_cells.size() + 1);
+  entries[0] = 0;
+  for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) {
+    entries[i_cell + 1] = entries[i_cell] + owned_cell_cells[i_cell].size();
+  }
+  Array<int> neighbors(entries[owned_cell_cells.size()]);
+  {
+    size_t i = 0;
+    for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) {
+      for (CellId cell_id : owned_cell_cells[i_cell]) {
+        neighbors[i] = global_cell_index[cell_id];
+        ++i;
+      }
+    }
+  }
+
+  return CRSGraph(entries, neighbors);
+}
+
 template std::ostream& Connectivity<1>::_write(std::ostream&) const;
 template std::ostream& Connectivity<2>::_write(std::ostream&) const;
 template std::ostream& Connectivity<3>::_write(std::ostream&) const;
@@ -333,3 +399,7 @@ template Connectivity<3>::Connectivity();
 template void Connectivity<1>::_buildFrom(const ConnectivityDescriptor&);
 template void Connectivity<2>::_buildFrom(const ConnectivityDescriptor&);
 template void Connectivity<3>::_buildFrom(const ConnectivityDescriptor&);
+
+template CRSGraph Connectivity<1>::ownCellToCellGraph() const;
+template CRSGraph Connectivity<2>::ownCellToCellGraph() const;
+template CRSGraph Connectivity<3>::ownCellToCellGraph() const;
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index aa1a471d6e84b554c1715ca7eecc7a2361744b07..6d89a820fc25f2ed662f0e4f277a7f44fb136830 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -59,8 +59,6 @@ class Connectivity final : public IConnectivity
   ConnectivityMatrix m_item_to_item_matrix[Dimension + 1][Dimension + 1];
   WeakCellValue<const CellType> m_cell_type;
 
-  WeakCellValue<const int> m_cell_global_index;
-
   WeakCellValue<const int> m_cell_number;
   WeakFaceValue<const int> m_face_number;
   WeakEdgeValue<const int> m_edge_number;
@@ -655,41 +653,7 @@ class Connectivity final : public IConnectivity
     }
   }
 
-  PUGS_INLINE
-  CRSGraph
-  cellToCellGraph() const
-  {
-    std::vector<std::set<int>> cell_cells(this->numberOfCells());
-    const auto& face_to_cell_matrix = this->faceToCellMatrix();
-
-    for (FaceId l = 0; l < this->numberOfFaces(); ++l) {
-      const auto& face_to_cell = face_to_cell_matrix[l];
-      if (face_to_cell.size() > 1) {
-        const CellId cell_0 = face_to_cell[0];
-        const CellId cell_1 = face_to_cell[1];
-
-        cell_cells[cell_0].insert(cell_1);
-        cell_cells[cell_1].insert(cell_0);
-      }
-    }
-
-    Array<int> entries(this->numberOfCells() + 1);
-    entries[0] = 0;
-    for (size_t j = 0; j < this->numberOfCells(); ++j) {
-      entries[j + 1] = entries[j] + cell_cells[j].size();
-    }
-    Array<int> neighbors(entries[this->numberOfCells()]);
-    {
-      size_t k = 0;
-      for (size_t j = 0; j < this->numberOfCells(); ++j) {
-        for (CellId cell_id : cell_cells[j]) {
-          neighbors[k] = m_cell_global_index[cell_id];
-          ++k;
-        }
-      }
-    }
-    return CRSGraph(entries, neighbors);
-  }
+  CRSGraph ownCellToCellGraph() const;
 
   PUGS_INLINE
   size_t
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index 65326b56bd49c1ce155dbf07fe32da8622c2e66c..1270dd78f39183d198e7155198e904799203398a 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -5,18 +5,29 @@
 #include <utils/Partitioner.hpp>
 
 #include <iostream>
-#include <unordered_map>
+#include <unordered_set>
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildNewOwner()
 {
   if constexpr (item_type == ItemType::cell) {
-    CRSGraph connectivity_graph = m_connectivity.cellToCellGraph();
+    CRSGraph connectivity_graph = m_connectivity.ownCellToCellGraph();
     Partitioner P;
+    Array new_owner_array = P.partition(connectivity_graph);
 
-    CellValue<int> cell_new_owner(m_connectivity, P.partition(connectivity_graph));
+    const auto& cell_is_owned = m_connectivity.cellIsOwned();
+
+    CellValue<int> cell_new_owner(m_connectivity);
+
+    size_t i = 0;
+    for (CellId cell_id = 0; cell_id < m_connectivity.numberOfCells(); ++cell_id) {
+      if (cell_is_owned[cell_id]) {
+        cell_new_owner[cell_id] = new_owner_array[i++];
+      }
+    }
+    synchronize(cell_new_owner);
 
     this->_dispatchedInfo<ItemType::cell>().m_new_owner = cell_new_owner;
   } else {
@@ -42,11 +53,12 @@ ConnectivityDispatcher<Dimension>::_buildNewOwner()
       });
 
     synchronize(item_new_owner);
+
     this->_dispatchedInfo<item_type>().m_new_owner = item_new_owner;
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildItemToExchangeLists()
@@ -59,7 +71,7 @@ ConnectivityDispatcher<Dimension>::_buildItemToExchangeLists()
   this->_buildRecvItemIdCorrespondanceByProc<item_type>();
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildItemListToSend()
@@ -67,27 +79,31 @@ ConnectivityDispatcher<Dimension>::_buildItemListToSend()
   if constexpr (item_type == ItemType::cell) {
     const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
     const auto& cell_to_node_matrix = m_connectivity.cellToNodeMatrix();
+    const auto& cell_is_owned       = m_connectivity.cellIsOwned();
 
     const auto& cell_new_owner = this->_dispatchedInfo<ItemType::cell>().m_new_owner;
 
     std::vector<std::vector<CellId>> cell_vector_to_send_by_proc(parallel::size());
     Array<bool> send_to_rank(parallel::size());
-    for (CellId j = 0; j < m_connectivity.numberOfCells(); ++j) {
+
+    for (CellId cell_id = 0; cell_id < m_connectivity.numberOfCells(); ++cell_id) {
       send_to_rank.fill(false);
-      const auto& cell_to_node = cell_to_node_matrix[j];
-
-      for (size_t R = 0; R < cell_to_node.size(); ++R) {
-        const NodeId& r          = cell_to_node[R];
-        const auto& node_to_cell = node_to_cell_matrix[r];
-        for (size_t K = 0; K < node_to_cell.size(); ++K) {
-          const CellId& k                 = node_to_cell[K];
-          send_to_rank[cell_new_owner[k]] = true;
+      const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+      if (cell_is_owned[cell_id]) {
+        for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+          const NodeId& node_id    = cell_to_node[i_node];
+          const auto& node_to_cell = node_to_cell_matrix[node_id];
+          for (size_t i_node_cell = 0; i_node_cell < node_to_cell.size(); ++i_node_cell) {
+            const CellId& node_cell                 = node_to_cell[i_node_cell];
+            send_to_rank[cell_new_owner[node_cell]] = true;
+          }
         }
-      }
 
-      for (size_t k = 0; k < send_to_rank.size(); ++k) {
-        if (send_to_rank[k]) {
-          cell_vector_to_send_by_proc[k].push_back(j);
+        for (size_t i_rank = 0; i_rank < send_to_rank.size(); ++i_rank) {
+          if (send_to_rank[i_rank]) {
+            cell_vector_to_send_by_proc[i_rank].push_back(cell_id);
+          }
         }
       }
     }
@@ -100,24 +116,141 @@ ConnectivityDispatcher<Dimension>::_buildItemListToSend()
   } else {
     const auto& cell_list_to_send_by_proc = this->_dispatchedInfo<ItemType::cell>().m_list_to_send_by_proc;
 
+    const auto& item_is_owned = m_connectivity.template isOwned<item_type>();
+    const auto& item_owner    = m_connectivity.template owner<item_type>();
+    const auto& item_number   = m_connectivity.template number<item_type>();
+
     using ItemId                        = ItemIdT<item_type>;
     const auto& cell_to_sub_item_matrix = m_connectivity.template getItemToItemMatrix<ItemType::cell, item_type>();
 
     auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
     item_list_to_send_by_proc.resize(parallel::size());
 
+    Array<bool> tag(m_connectivity.template numberOf<item_type>());
+
+    struct ToNumber
+    {
+      int to;
+      int number;
+      bool
+      operator==(const ToNumber& x) const
+      {
+        return to == x.to and number == x.number;
+      }
+      ToNumber() = default;
+      ToNumber(int _to, int _number) : to{_to}, number{_number} {}
+    };
+
+    struct hashFunction
+    {
+      size_t
+      operator()(const ToNumber& x) const
+      {
+        return x.to ^ x.number;
+      }
+    };
+
+    std::vector<std::unordered_set<ToNumber, hashFunction>> request_not_owned_to_send_by_proc(parallel::size());
+
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      tag.fill(false);
+      for (size_t i_cell = 0; i_cell < cell_list_to_send_by_proc[i_rank].size(); ++i_cell) {
+        const CellId& cell_id          = cell_list_to_send_by_proc[i_rank][i_cell];
+        const auto& cell_sub_item_list = cell_to_sub_item_matrix[cell_id];
+        for (size_t i_item = 0; i_item < cell_sub_item_list.size(); ++i_item) {
+          const ItemId& sub_item_id = cell_sub_item_list[i_item];
+          if (not item_is_owned[sub_item_id] and (not tag[sub_item_id])) {
+            request_not_owned_to_send_by_proc[item_owner[sub_item_id]].insert(
+              ToNumber{static_cast<int>(i_rank), item_number[sub_item_id]});
+            tag[sub_item_id] = true;
+          }
+        }
+      }
+    }
+
+    std::vector<Array<int>> not_owned_info_to_send_by_proc(parallel::size());
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      Array<int> info_to_send_by_proc(2 * request_not_owned_to_send_by_proc[i_rank].size());
+      size_t i = 0;
+      for (auto&& to_send_info : request_not_owned_to_send_by_proc[i_rank]) {
+        info_to_send_by_proc[i]     = to_send_info.to;
+        info_to_send_by_proc[i + 1] = to_send_info.number;
+        i += 2;
+      }
+      not_owned_info_to_send_by_proc[i_rank] = info_to_send_by_proc;
+    }
+
+    std::vector<Array<int>> number_of_not_owned_to_send_by_proc(parallel::size());
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      number_of_not_owned_to_send_by_proc[i_rank]    = Array<int>(1);
+      number_of_not_owned_to_send_by_proc[i_rank][0] = request_not_owned_to_send_by_proc[i_rank].size();
+    }
+
+    std::vector<Array<int>> number_of_possibly_item_to_send_by_proc(parallel::size());
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      number_of_possibly_item_to_send_by_proc[i_rank] = Array<int>(1);
+    }
+    parallel::exchange(number_of_not_owned_to_send_by_proc, number_of_possibly_item_to_send_by_proc);
+
+    std::vector<Array<int>> possible_info_to_send_demanded_by_rank(parallel::size());
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      possible_info_to_send_demanded_by_rank[i_rank] =
+        Array<int>(2 * number_of_possibly_item_to_send_by_proc[i_rank][0]);
+    }
+
+    parallel::exchange(not_owned_info_to_send_by_proc, possible_info_to_send_demanded_by_rank);
+
+    std::vector<std::vector<int>> possible_item_to_send_to_rank(parallel::size());
+    for (size_t i_demanding_rank = 0; i_demanding_rank < parallel::size(); ++i_demanding_rank) {
+      const auto& possible_info_to_send_to = possible_info_to_send_demanded_by_rank[i_demanding_rank];
+      for (size_t i = 0; i < possible_info_to_send_to.size(); i += 2) {
+        const int i_rank = possible_info_to_send_to[i];
+        const int number = possible_info_to_send_to[i + 1];
+        possible_item_to_send_to_rank[i_rank].push_back(number);
+      }
+    }
+
+    const bool has_possible_item_to_send = [&] {
+      for (auto&& item_to_send_list : possible_item_to_send_to_rank) {
+        if (item_to_send_list.size() > 0) {
+          return true;
+        }
+      }
+      return false;
+    }();
+
+    std::unordered_map<int, ItemId> owned_number_to_item_id_map;
+    if (has_possible_item_to_send) {
+      for (ItemId item_id = 0; item_id < m_connectivity.template numberOf<item_type>(); ++item_id) {
+        if (item_is_owned[item_id]) {
+          owned_number_to_item_id_map[item_number[item_id]] = item_id;
+        }
+      }
+    }
+
     for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-      Array<bool> tag(m_connectivity.template numberOf<item_type>());
       tag.fill(false);
       std::vector<ItemId> item_id_vector;
-      for (size_t j = 0; j < cell_list_to_send_by_proc[i_rank].size(); ++j) {
-        const CellId& cell_id          = cell_list_to_send_by_proc[i_rank][j];
+      const auto& possible_item_to_send = possible_item_to_send_to_rank[i_rank];
+
+      for (auto&& possible_item_number : possible_item_to_send) {
+        const auto& searched_item_id = owned_number_to_item_id_map.find(possible_item_number);
+        Assert(searched_item_id != owned_number_to_item_id_map.end());
+        const ItemId item_id = searched_item_id->second;
+        if (not tag[item_id]) {
+          item_id_vector.push_back(item_id);
+          tag[item_id] = true;
+        }
+      }
+
+      for (size_t i_cell = 0; i_cell < cell_list_to_send_by_proc[i_rank].size(); ++i_cell) {
+        const CellId& cell_id          = cell_list_to_send_by_proc[i_rank][i_cell];
         const auto& cell_sub_item_list = cell_to_sub_item_matrix[cell_id];
-        for (size_t r = 0; r < cell_sub_item_list.size(); ++r) {
-          const ItemId& item_id = cell_sub_item_list[r];
-          if (not tag[item_id]) {
-            item_id_vector.push_back(item_id);
-            tag[item_id] = true;
+        for (size_t i_item = 0; i_item < cell_sub_item_list.size(); ++i_item) {
+          const ItemId& sub_item_id = cell_sub_item_list[i_item];
+          if (item_is_owned[sub_item_id] and (not tag[sub_item_id])) {
+            item_id_vector.push_back(sub_item_id);
+            tag[sub_item_id] = true;
           }
         }
       }
@@ -126,7 +259,7 @@ ConnectivityDispatcher<Dimension>::_buildItemListToSend()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildNumberOfItemToExchange()
@@ -141,7 +274,7 @@ ConnectivityDispatcher<Dimension>::_buildNumberOfItemToExchange()
   this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc = parallel::allToAll(nb_item_to_send_by_proc);
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
 ConnectivityDispatcher<Dimension>::_gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
@@ -150,6 +283,7 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(const ItemValue<DataType, item_ty
   std::vector<Array<const DataType>> recv_item_data_by_proc = this->exchange(data_to_gather);
 
   const auto& recv_id_correspondance_by_proc = this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
+
   Assert(recv_id_correspondance_by_proc.size() == parallel::size());
 
   gathered_array = Array<std::remove_const_t<DataType>>(this->_dispatchedInfo<item_type>().m_number_to_id_map.size());
@@ -162,7 +296,7 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(const ItemValue<DataType, item_ty
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
 void
 ConnectivityDispatcher<Dimension>::_gatherFrom(
@@ -221,14 +355,14 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 void
 ConnectivityDispatcher<Dimension>::_buildCellNumberIdMap()
 {
   const auto recv_cell_number_by_proc = this->exchange(m_connectivity.template number<ItemType::cell>());
   auto& cell_number_id_map            = this->_dispatchedInfo<ItemType::cell>().m_number_to_id_map;
+  CellId cell_id                      = 0;
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-    CellId cell_id = 0;
     for (size_t i = 0; i < recv_cell_number_by_proc[i_rank].size(); ++i) {
       const int cell_number     = recv_cell_number_by_proc[i_rank][i];
       auto [iterator, inserted] = cell_number_id_map.insert(std::make_pair(cell_number, cell_id));
@@ -238,7 +372,7 @@ ConnectivityDispatcher<Dimension>::_buildCellNumberIdMap()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename ItemOfItemT>
 void
 ConnectivityDispatcher<Dimension>::_buildSubItemNumberToIdMap()
@@ -250,18 +384,20 @@ ConnectivityDispatcher<Dimension>::_buildSubItemNumberToIdMap()
     this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
 
   auto& sub_item_number_id_map = this->_dispatchedInfo<ItemOfItemT::sub_item_type>().m_number_to_id_map;
+  int sub_item_id              = 0;
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-    int sub_item_id = 0;
     for (size_t i = 0; i < cell_sub_item_number_to_recv_by_proc[i_rank].size(); ++i) {
-      int sub_item_number       = cell_sub_item_number_to_recv_by_proc[i_rank][i];
+      int sub_item_number = cell_sub_item_number_to_recv_by_proc[i_rank][i];
+
       auto [iterator, inserted] = sub_item_number_id_map.insert(std::make_pair(sub_item_number, sub_item_id));
-      if (inserted)
-        sub_item_id++;
+      if (inserted) {
+        ++sub_item_id;
+      }
     }
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename SubItemOfItemT>
 void
 ConnectivityDispatcher<Dimension>::_buildNumberOfSubItemPerItemToRecvByProc()
@@ -280,7 +416,7 @@ ConnectivityDispatcher<Dimension>::_buildNumberOfSubItemPerItemToRecvByProc()
     this->exchange(number_of_sub_item_per_item);
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename SubItemOfItemT>
 void
 ConnectivityDispatcher<Dimension>::_buildSubItemNumbersToRecvByProc()
@@ -330,7 +466,7 @@ ConnectivityDispatcher<Dimension>::_buildSubItemNumbersToRecvByProc()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename ItemOfItemT>
 void
 ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
@@ -340,6 +476,8 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
 
   const auto& item_list_to_recv_size_by_proc = this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
 
+  const auto& recv_id_correspondance_by_proc = this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
+
   const auto& number_of_sub_item_per_item_to_recv_by_proc =
     this->_dispatchedInfo<ItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc;
 
@@ -349,7 +487,14 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
     this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
 
   std::vector<std::vector<unsigned int>> item_to_subitem_legacy;
-  size_t number_of_node_by_cell = 0;
+  item_to_subitem_legacy.resize([&] {
+    size_t size = 0;
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      size += item_list_to_recv_size_by_proc[i_rank];
+    }
+    return size;
+  }());
+  size_t number_of_subitem_by_item = 0;
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
     int l = 0;
     for (size_t i = 0; i < item_list_to_recv_size_by_proc[i_rank]; ++i) {
@@ -359,16 +504,14 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
         Assert(searched_sub_item_id != sub_item_number_id_map.end());
         sub_item_vector.push_back(searched_sub_item_id->second);
       }
-      number_of_node_by_cell += sub_item_vector.size();
+      number_of_subitem_by_item += sub_item_vector.size();
 
-      item_to_subitem_legacy.emplace_back(sub_item_vector);
+      item_to_subitem_legacy[recv_id_correspondance_by_proc[i_rank][i]] = std::move(sub_item_vector);
     }
   }
-  Array<unsigned int> item_to_subitem_row_map(item_to_subitem_legacy.size() + 1);
-  Array<unsigned int> item_to_subitem_list(number_of_node_by_cell);
 
-  item_to_subitem_row_map.fill(10000000);
-  item_to_subitem_list.fill(10000000);
+  Array<unsigned int> item_to_subitem_row_map(item_to_subitem_legacy.size() + 1);
+  Array<unsigned int> item_to_subitem_list(number_of_subitem_by_item);
 
   item_to_subitem_row_map[0] = 0;
   for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
@@ -385,7 +528,7 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
   m_new_descriptor.itemOfItemVector<ItemOfItemT>() = ConnectivityMatrix(item_to_subitem_row_map, item_to_subitem_list);
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildRecvItemIdCorrespondanceByProc()
@@ -413,6 +556,7 @@ ConnectivityDispatcher<Dimension>::_buildRecvItemIdCorrespondanceByProc()
   parallel::exchange(send_item_number_by_proc, recv_item_number_by_proc);
 
   const auto& item_number_to_id_map = this->_dispatchedInfo<item_type>().m_number_to_id_map;
+
   for (size_t i_rank = 0; i_rank < item_list_to_recv_size_by_proc.size(); ++i_rank) {
     Array<ItemId> item_id_correspondance(item_list_to_recv_size_by_proc[i_rank]);
     for (size_t l = 0; l < item_list_to_recv_size_by_proc[i_rank]; ++l) {
@@ -426,7 +570,7 @@ ConnectivityDispatcher<Dimension>::_buildRecvItemIdCorrespondanceByProc()
   this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc = recv_item_id_correspondance_by_proc;
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
@@ -446,162 +590,153 @@ ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
   }();
 
   if (number_of_item_list_sender > 0) {
-    if (number_of_item_list_sender > 1) {
-      std::cerr << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red
-                << "need to check that knowing procs know the same item_ref_lists!" << rang::fg::reset << '\n';
-    }
-
-    if (number_of_item_list_sender < parallel::size()) {
-      const size_t sender_rank = [&]() {
-        size_t i_rank = 0;
-        for (; i_rank < parallel::size(); ++i_rank) {
-          if (number_of_item_ref_list_per_proc[i_rank] > 0) {
-            break;
-          }
-        }
-        return i_rank;
-      }();
-
-      Assert(number_of_item_list_sender < parallel::size());
-
-      // sending is boundary property
-      Array<RefItemListBase::Type> ref_item_list_type{number_of_item_ref_list_per_proc[sender_rank]};
-      if (parallel::rank() == sender_rank) {
-        for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
-             ++i_item_ref_list) {
-          auto item_ref_list                  = m_connectivity.template refItemList<item_type>(i_item_ref_list);
-          ref_item_list_type[i_item_ref_list] = item_ref_list.type();
+    const size_t sender_rank = [&]() {
+      size_t i_rank = 0;
+      for (; i_rank < parallel::size(); ++i_rank) {
+        if (number_of_item_ref_list_per_proc[i_rank] > 0) {
+          break;
         }
       }
-      parallel::broadcast(ref_item_list_type, sender_rank);
-
-      // sending references tags
-      Array<RefId::TagNumberType> ref_tag_list{number_of_item_ref_list_per_proc[sender_rank]};
-      if (parallel::rank() == sender_rank) {
-        for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
-             ++i_item_ref_list) {
-          auto item_ref_list            = m_connectivity.template refItemList<item_type>(i_item_ref_list);
-          ref_tag_list[i_item_ref_list] = item_ref_list.refId().tagNumber();
-        }
+      return i_rank;
+    }();
+
+    // sending is boundary property
+    Array<RefItemListBase::Type> ref_item_list_type{number_of_item_ref_list_per_proc[sender_rank]};
+    if (parallel::rank() == sender_rank) {
+      for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
+           ++i_item_ref_list) {
+        auto item_ref_list                  = m_connectivity.template refItemList<item_type>(i_item_ref_list);
+        ref_item_list_type[i_item_ref_list] = item_ref_list.type();
       }
-      parallel::broadcast(ref_tag_list, sender_rank);
-
-      // sending references name size
-      Array<size_t> ref_name_size_list{number_of_item_ref_list_per_proc[sender_rank]};
-      if (parallel::rank() == sender_rank) {
-        for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
-             ++i_item_ref_list) {
-          auto item_ref_list                  = m_connectivity.template refItemList<item_type>(i_item_ref_list);
-          ref_name_size_list[i_item_ref_list] = item_ref_list.refId().tagName().size();
-        }
+    }
+    parallel::broadcast(ref_item_list_type, sender_rank);
+
+    // sending references tags
+    Array<RefId::TagNumberType> ref_tag_list{number_of_item_ref_list_per_proc[sender_rank]};
+    if (parallel::rank() == sender_rank) {
+      for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
+           ++i_item_ref_list) {
+        auto item_ref_list            = m_connectivity.template refItemList<item_type>(i_item_ref_list);
+        ref_tag_list[i_item_ref_list] = item_ref_list.refId().tagNumber();
       }
-      parallel::broadcast(ref_name_size_list, sender_rank);
-
-      // sending references name size
-      Array<RefId::TagNameType::value_type> ref_name_cat{sum(ref_name_size_list)};
-      if (parallel::rank() == sender_rank) {
-        size_t i_char = 0;
-        for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
-             ++i_item_ref_list) {
-          auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
-          for (auto c : item_ref_list.refId().tagName()) {
-            ref_name_cat[i_char++] = c;
-          }
-        }
+    }
+    parallel::broadcast(ref_tag_list, sender_rank);
+
+    // sending references name size
+    Array<size_t> ref_name_size_list{number_of_item_ref_list_per_proc[sender_rank]};
+    if (parallel::rank() == sender_rank) {
+      for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
+           ++i_item_ref_list) {
+        auto item_ref_list                  = m_connectivity.template refItemList<item_type>(i_item_ref_list);
+        ref_name_size_list[i_item_ref_list] = item_ref_list.refId().tagName().size();
       }
-      parallel::broadcast(ref_name_cat, sender_rank);
-
-      std::vector<RefId> ref_id_list = [&]() {
-        std::vector<RefId> mutable_ref_id_list;
-        mutable_ref_id_list.reserve(ref_name_size_list.size());
-        size_t begining = 0;
-        for (size_t i_ref = 0; i_ref < ref_name_size_list.size(); ++i_ref) {
-          const size_t size = ref_name_size_list[i_ref];
-          mutable_ref_id_list.emplace_back(ref_tag_list[i_ref], std::string{&(ref_name_cat[begining]), size});
-          begining += size;
+    }
+    parallel::broadcast(ref_name_size_list, sender_rank);
+
+    // sending references name size
+    Array<RefId::TagNameType::value_type> ref_name_cat{sum(ref_name_size_list)};
+    if (parallel::rank() == sender_rank) {
+      size_t i_char = 0;
+      for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
+           ++i_item_ref_list) {
+        auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
+        for (auto c : item_ref_list.refId().tagName()) {
+          ref_name_cat[i_char++] = c;
         }
-        return mutable_ref_id_list;
-      }();
-
-      using block_type            = int32_t;
-      constexpr size_t block_size = sizeof(block_type);
-      const size_t nb_block       = ref_id_list.size() / block_size + (ref_id_list.size() % block_size != 0);
-      for (size_t i_block = 0; i_block < nb_block; ++i_block) {
-        ItemValue<block_type, item_type> item_references(m_connectivity);
-        item_references.fill(0);
-
-        if (m_connectivity.template numberOfRefItemList<item_type>() > 0) {
-          const size_t max_i_ref = std::min(ref_id_list.size(), block_size * (i_block + 1));
-          for (size_t i_ref = block_size * i_block, i = 0; i_ref < max_i_ref; ++i_ref, ++i) {
-            block_type ref_bit{1 << i};
-            auto item_ref_list = m_connectivity.template refItemList<item_type>(i_ref);
-
-            const auto& item_list = item_ref_list.list();
-            for (size_t i_item = 0; i_item < item_list.size(); ++i_item) {
-              const ItemId& item_id = item_list[i_item];
-              item_references[item_id] |= ref_bit;
-            }
+      }
+    }
+    parallel::broadcast(ref_name_cat, sender_rank);
+
+    std::vector<RefId> ref_id_list = [&]() {
+      std::vector<RefId> mutable_ref_id_list;
+      mutable_ref_id_list.reserve(ref_name_size_list.size());
+      size_t begining = 0;
+      for (size_t i_ref = 0; i_ref < ref_name_size_list.size(); ++i_ref) {
+        const size_t size = ref_name_size_list[i_ref];
+        mutable_ref_id_list.emplace_back(ref_tag_list[i_ref], std::string{&(ref_name_cat[begining]), size});
+        begining += size;
+      }
+      return mutable_ref_id_list;
+    }();
+
+    using block_type            = int32_t;
+    constexpr size_t block_size = sizeof(block_type);
+    const size_t nb_block       = ref_id_list.size() / block_size + (ref_id_list.size() % block_size != 0);
+    for (size_t i_block = 0; i_block < nb_block; ++i_block) {
+      ItemValue<block_type, item_type> item_references(m_connectivity);
+      item_references.fill(0);
+
+      if (m_connectivity.template numberOfRefItemList<item_type>() > 0) {
+        const size_t max_i_ref = std::min(ref_id_list.size(), block_size * (i_block + 1));
+        for (size_t i_ref = block_size * i_block, i = 0; i_ref < max_i_ref; ++i_ref, ++i) {
+          block_type ref_bit{1 << i};
+          auto item_ref_list = m_connectivity.template refItemList<item_type>(i_ref);
+
+          const auto& item_list = item_ref_list.list();
+          for (size_t i_item = 0; i_item < item_list.size(); ++i_item) {
+            const ItemId& item_id = item_list[i_item];
+            item_references[item_id] |= ref_bit;
           }
         }
+      }
 
-        const auto& nb_item_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_size_by_proc;
+      const auto& nb_item_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_size_by_proc;
 
-        const auto& send_item_id_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
+      const auto& send_item_id_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
 
-        std::vector<Array<const block_type>> send_item_refs_by_proc(parallel::size());
+      std::vector<Array<const block_type>> send_item_refs_by_proc(parallel::size());
 
-        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-          Array<block_type> send_item_refs(nb_item_to_send_by_proc[i_rank]);
-          const Array<const ItemId> send_item_id = send_item_id_by_proc[i_rank];
-          parallel_for(
-            send_item_id.size(), PUGS_LAMBDA(size_t l) {
-              const ItemId& item_id = send_item_id[l];
-              send_item_refs[l]     = item_references[item_id];
-            });
-          send_item_refs_by_proc[i_rank] = send_item_refs;
-        }
+      for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+        Array<block_type> send_item_refs(nb_item_to_send_by_proc[i_rank]);
+        const Array<const ItemId> send_item_id = send_item_id_by_proc[i_rank];
+        parallel_for(
+          send_item_id.size(), PUGS_LAMBDA(size_t l) {
+            const ItemId& item_id = send_item_id[l];
+            send_item_refs[l]     = item_references[item_id];
+          });
+        send_item_refs_by_proc[i_rank] = send_item_refs;
+      }
 
-        std::vector<Array<block_type>> recv_item_refs_by_proc(parallel::size());
-        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-          recv_item_refs_by_proc[i_rank] =
-            Array<block_type>(this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc[i_rank]);
-        }
-        parallel::exchange(send_item_refs_by_proc, recv_item_refs_by_proc);
-
-        const auto& recv_item_id_correspondance_by_proc =
-          this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
-        std::vector<block_type> item_refs(m_new_descriptor.template itemNumberVector<item_type>().size());
-        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-          for (size_t r = 0; r < recv_item_refs_by_proc[i_rank].size(); ++r) {
-            const ItemId& item_id = recv_item_id_correspondance_by_proc[i_rank][r];
-            item_refs[item_id]    = recv_item_refs_by_proc[i_rank][r];
-          }
+      std::vector<Array<block_type>> recv_item_refs_by_proc(parallel::size());
+      for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+        recv_item_refs_by_proc[i_rank] =
+          Array<block_type>(this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc[i_rank]);
+      }
+      parallel::exchange(send_item_refs_by_proc, recv_item_refs_by_proc);
+
+      const auto& recv_item_id_correspondance_by_proc =
+        this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
+      std::vector<block_type> item_refs(m_new_descriptor.template itemNumberVector<item_type>().size());
+      for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+        for (size_t r = 0; r < recv_item_refs_by_proc[i_rank].size(); ++r) {
+          const ItemId& item_id = recv_item_id_correspondance_by_proc[i_rank][r];
+          item_refs[item_id]    = recv_item_refs_by_proc[i_rank][r];
         }
+      }
 
-        const size_t max_i_ref = std::min(ref_id_list.size(), block_size * (i_block + 1));
-        for (size_t i_ref = block_size * i_block, i = 0; i_ref < max_i_ref; ++i_ref, ++i) {
-          block_type ref_bit{1 << i};
+      const size_t max_i_ref = std::min(ref_id_list.size(), block_size * (i_block + 1));
+      for (size_t i_ref = block_size * i_block, i = 0; i_ref < max_i_ref; ++i_ref, ++i) {
+        block_type ref_bit{1 << i};
 
-          std::vector<ItemId> item_id_vector;
+        std::vector<ItemId> item_id_vector;
 
-          for (uint32_t i_item = 0; i_item < item_refs.size(); ++i_item) {
-            const ItemId item_id{i_item};
-            if (item_refs[item_id] & ref_bit) {
-              item_id_vector.push_back(item_id);
-            }
+        for (uint32_t i_item = 0; i_item < item_refs.size(); ++i_item) {
+          const ItemId item_id{i_item};
+          if (item_refs[item_id] & ref_bit) {
+            item_id_vector.push_back(item_id);
           }
+        }
 
-          Array<const ItemId> item_id_array = convert_to_array(item_id_vector);
+        Array<const ItemId> item_id_array = convert_to_array(item_id_vector);
 
-          RefItemListBase::Type type = ref_item_list_type[i_ref];
-          m_new_descriptor.addRefItemList(RefItemList<item_type>(ref_id_list[i_ref], item_id_array, type));
-        }
+        RefItemListBase::Type type = ref_item_list_type[i_ref];
+        m_new_descriptor.addRefItemList(RefItemList<item_type>(ref_id_list[i_ref], item_id_array, type));
       }
     }
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 void
 ConnectivityDispatcher<Dimension>::_dispatchEdges()
 {
@@ -611,6 +746,10 @@ ConnectivityDispatcher<Dimension>::_dispatchEdges()
     this->_buildSubItemNumberToIdMap<EdgeOfCell>();
     this->_buildItemToExchangeLists<ItemType::edge>();
 
+    this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfEdge>();
+    this->_buildSubItemNumbersToRecvByProc<NodeOfEdge>();
+    this->_buildItemToSubItemDescriptor<NodeOfEdge>();
+
     m_new_descriptor.setEdgeNumberVector([&] {
       Array<int> edge_number_vector;
       this->_gatherFrom(m_connectivity.template number<ItemType::edge>(), edge_number_vector);
@@ -619,10 +758,6 @@ ConnectivityDispatcher<Dimension>::_dispatchEdges()
 
     this->_buildItemToSubItemDescriptor<EdgeOfCell>();
 
-    this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfEdge>();
-    this->_buildSubItemNumbersToRecvByProc<NodeOfEdge>();
-    this->_buildItemToSubItemDescriptor<NodeOfEdge>();
-
     this->_buildNumberOfSubItemPerItemToRecvByProc<EdgeOfFace>();
     this->_buildSubItemNumbersToRecvByProc<EdgeOfFace>();
     this->_buildItemToSubItemDescriptor<EdgeOfFace>();
@@ -643,7 +778,7 @@ ConnectivityDispatcher<Dimension>::_dispatchEdges()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 void
 ConnectivityDispatcher<Dimension>::_dispatchFaces()
 {
@@ -681,7 +816,7 @@ ConnectivityDispatcher<Dimension>::_dispatchFaces()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType& connectivity)
   : m_connectivity(connectivity)
 {
@@ -693,12 +828,15 @@ ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType
   }
 
   this->_buildNewOwner<ItemType::cell>();
+
   if constexpr (Dimension > 1) {
     this->_buildNewOwner<ItemType::face>();
   }
+
   if constexpr (Dimension > 2) {
     this->_buildNewOwner<ItemType::edge>();
   }
+
   this->_buildNewOwner<ItemType::node>();
 
   this->_buildItemToExchangeLists<ItemType::cell>();
@@ -750,7 +888,6 @@ ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType
   this->_dispatchEdges();
 
   this->_buildItemReferenceList<ItemType::node>();
-
   m_dispatched_connectivity = ConnectivityType::build(m_new_descriptor);
 
   {
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index 9ff5913b60799148f99bb6fdf4cff2cedd56973c..46932004c5d585214f71f4395d6719ee2a529835 100644
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -8,7 +8,7 @@
 
 #include <unordered_map>
 
-template <int Dimension>
+template <size_t Dimension>
 class ConnectivityDispatcher
 {
  public:
@@ -253,19 +253,19 @@ class ConnectivityDispatcher
 
     using MutableDataType = std::remove_const_t<DataType>;
     std::vector<Array<DataType>> item_value_to_send_by_proc(parallel::size());
-    for (size_t i = 0; i < parallel::size(); ++i) {
-      const Array<const ItemId>& item_list = item_list_to_send_by_proc[i];
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      const Array<const ItemId>& item_list = item_list_to_send_by_proc[i_rank];
       Array<MutableDataType> item_value_list(item_list.size());
       parallel_for(
         item_list.size(),
         PUGS_LAMBDA(const ItemId& item_id) { item_value_list[item_id] = item_value[item_list[item_id]]; });
-      item_value_to_send_by_proc[i] = item_value_list;
+      item_value_to_send_by_proc[i_rank] = item_value_list;
     }
 
     const auto& item_list_to_recv_size_by_proc = this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
     std::vector<Array<MutableDataType>> recv_item_value_by_proc(parallel::size());
-    for (size_t i = 0; i < parallel::size(); ++i) {
-      recv_item_value_by_proc[i] = Array<MutableDataType>(item_list_to_recv_size_by_proc[i]);
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      recv_item_value_by_proc[i_rank] = Array<MutableDataType>(item_list_to_recv_size_by_proc[i_rank]);
     }
 
     parallel::exchange(item_value_to_send_by_proc, recv_item_value_by_proc);
@@ -277,14 +277,68 @@ class ConnectivityDispatcher
       const auto& recv_item_id_correspondance = recv_item_id_correspondance_by_proc[i_rank];
       const auto& recv_item_value             = recv_item_value_by_proc[i_rank];
       parallel_for(
-        recv_item_value.size(), PUGS_LAMBDA(size_t r) {
-          const ItemId& item_id   = recv_item_id_correspondance[r];
-          new_item_value[item_id] = recv_item_value[r];
+        recv_item_value.size(), PUGS_LAMBDA(size_t index) {
+          const ItemId item_id    = recv_item_id_correspondance[index];
+          new_item_value[item_id] = recv_item_value[index];
         });
     }
     return new_item_value;
   }
 
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
+  dispatch(ItemArray<DataType, item_type, ConnectivityPtr> item_array) const
+  {
+    using ItemId = ItemIdT<item_type>;
+
+    Assert(m_dispatched_connectivity.use_count() > 0, "cannot dispatch quantity before connectivity");
+
+    const auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
+
+    const size_t size_of_arrays = item_array.sizeOfArrays();
+
+    using MutableDataType = std::remove_const_t<DataType>;
+    std::vector<Array<DataType>> item_array_to_send_by_proc(parallel::size());
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      const Array<const ItemId>& item_list = item_list_to_send_by_proc[i_rank];
+      Array<MutableDataType> item_array_list(item_list.size() * item_array.sizeOfArrays());
+      parallel_for(
+        item_list.size(), PUGS_LAMBDA(const ItemId& i) {
+          const size_t j   = i * size_of_arrays;
+          const auto array = item_array[item_list[i]];
+          for (size_t k = 0; k < size_of_arrays; ++k) {
+            item_array_list[j + k] = array[k];
+          }
+        });
+      item_array_to_send_by_proc[i_rank] = item_array_list;
+    }
+
+    const auto& item_list_to_recv_size_by_proc = this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
+    std::vector<Array<MutableDataType>> recv_item_array_by_proc(parallel::size());
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      recv_item_array_by_proc[i_rank] = Array<MutableDataType>(item_list_to_recv_size_by_proc[i_rank] * size_of_arrays);
+    }
+
+    parallel::exchange(item_array_to_send_by_proc, recv_item_array_by_proc);
+
+    const auto& recv_item_id_correspondance_by_proc =
+      this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
+    ItemArray<MutableDataType, item_type> new_item_array(*m_dispatched_connectivity, size_of_arrays);
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      const auto& recv_item_id_correspondance = recv_item_id_correspondance_by_proc[i_rank];
+      const auto& recv_item_array             = recv_item_array_by_proc[i_rank];
+      parallel_for(
+        recv_item_id_correspondance.size(), PUGS_LAMBDA(size_t i) {
+          const size_t j = i * size_of_arrays;
+          auto array     = new_item_array[recv_item_id_correspondance[i]];
+          for (size_t k = 0; k < size_of_arrays; ++k) {
+            array[k] = recv_item_array[j + k];
+          }
+        });
+    }
+    return new_item_array;
+  }
+
   ConnectivityDispatcher(const ConnectivityType& mesh);
   ConnectivityDispatcher(const ConnectivityDispatcher&) = delete;
   ~ConnectivityDispatcher()                             = default;
diff --git a/src/mesh/ConnectivityDispatcherVariant.hpp b/src/mesh/ConnectivityDispatcherVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..70ff61f2aa0ba20d738e16cd67f2bdaf7172a03e
--- /dev/null
+++ b/src/mesh/ConnectivityDispatcherVariant.hpp
@@ -0,0 +1,71 @@
+#ifndef CONNECTIVITY_DISPATCHER_VARIANT_HPP
+#define CONNECTIVITY_DISPATCHER_VARIANT_HPP
+
+#include <utils/Demangle.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <rang.hpp>
+
+#include <memory>
+#include <sstream>
+#include <variant>
+
+template <size_t Dimension>
+class ConnectivityDispatcher;
+
+class ConnectivityDispatcherVariant
+{
+ private:
+  using Variant = std::variant<std::shared_ptr<const ConnectivityDispatcher<1>>,   //
+                               std::shared_ptr<const ConnectivityDispatcher<2>>,   //
+                               std::shared_ptr<const ConnectivityDispatcher<3>>>;
+
+  Variant m_p_connectivity_dispatcher;
+
+ public:
+  template <size_t Dimension>
+  PUGS_INLINE std::shared_ptr<const ConnectivityDispatcher<Dimension>>
+  get() const
+  {
+    if (not std::holds_alternative<std::shared_ptr<const ConnectivityDispatcher<Dimension>>>(
+          this->m_p_connectivity_dispatcher)) {
+      std::ostringstream error_msg;
+      error_msg << "invalid connectivity dispatcher type type\n";
+      error_msg << "- required " << rang::fgB::red << demangle<ConnectivityDispatcher<Dimension>>() << rang::fg::reset
+                << '\n';
+      error_msg << "- contains " << rang::fgB::yellow
+                << std::visit(
+                     [](auto&& p_connectivity_dispatcher) -> std::string {
+                       using CDType = std::decay_t<decltype(*p_connectivity_dispatcher)>;
+                       return demangle<std::decay_t<CDType>>();
+                     },
+                     this->m_p_connectivity_dispatcher)
+                << rang::fg::reset;
+      throw NormalError(error_msg.str());
+    }
+    return std::get<std::shared_ptr<const ConnectivityDispatcher<Dimension>>>(m_p_connectivity_dispatcher);
+  }
+
+  PUGS_INLINE
+  Variant
+  variant() const
+  {
+    return m_p_connectivity_dispatcher;
+  }
+
+  ConnectivityDispatcherVariant() = delete;
+
+  template <size_t Dimension>
+  ConnectivityDispatcherVariant(
+    const std::shared_ptr<const ConnectivityDispatcher<Dimension>>& p_connectivity_dispatcher)
+    : m_p_connectivity_dispatcher{p_connectivity_dispatcher}
+  {}
+
+  ConnectivityDispatcherVariant(const ConnectivityDispatcherVariant&) = default;
+  ConnectivityDispatcherVariant(ConnectivityDispatcherVariant&&)      = default;
+
+  ~ConnectivityDispatcherVariant() = default;
+};
+
+#endif   // CONNECTIVITY_DISPATCHER_VARIANT_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index bb409dd28c7e80b12537beb584b421e81045893a..a2d147da4d747e294ac7b0757df875c732654d2d 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -1086,15 +1086,15 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
 
     switch (mesh_dimension) {
     case 1: {
-      this->_dispatch<1>();
+      this->_dispatch<1>(DispatchType::initial);
       break;
     }
     case 2: {
-      this->_dispatch<2>();
+      this->_dispatch<2>(DispatchType::initial);
       break;
     }
     case 3: {
-      this->_dispatch<3>();
+      this->_dispatch<3>(DispatchType::initial);
       break;
     }
     default: {
diff --git a/src/mesh/MeshBalancer.cpp b/src/mesh/MeshBalancer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..23e4759b09bb4b9b440f9488964445a09c9f02b2
--- /dev/null
+++ b/src/mesh/MeshBalancer.cpp
@@ -0,0 +1,17 @@
+#include <mesh/MeshBalancer.hpp>
+
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+
+MeshBalancer::MeshBalancer(const std::shared_ptr<const MeshVariant>& initial_mesh)
+{
+  m_mesh = initial_mesh;
+  std::visit(
+    [this](auto&& mesh) {
+      using MeshType = mesh_type_t<decltype(mesh)>;
+
+      this->_dispatch<MeshType::Dimension>(DispatchType::balance);
+    },
+    initial_mesh->variant());
+}
diff --git a/src/mesh/MeshBalancer.hpp b/src/mesh/MeshBalancer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..08abbd08190e7962e7ac02ffe6a6a3c810b293f4
--- /dev/null
+++ b/src/mesh/MeshBalancer.hpp
@@ -0,0 +1,13 @@
+#ifndef MESH_BALANCER_HPP
+#define MESH_BALANCER_HPP
+
+#include <mesh/MeshBuilderBase.hpp>
+
+class MeshBalancer : public MeshBuilderBase
+{
+ public:
+  MeshBalancer(const std::shared_ptr<const MeshVariant>& initial_mesh);
+  ~MeshBalancer() = default;
+};
+
+#endif   // MESH_BALANCER_HPP
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index 4537f83b2177210f2bb0bf7e8e55b3219b57e6f8..b8efb24e95ad3abcb0a46581e469d5bf86afdf31 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -3,6 +3,7 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/ConnectivityDescriptor.hpp>
 #include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ConnectivityDispatcherVariant.hpp>
 #include <mesh/ItemId.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
@@ -15,35 +16,45 @@
 
 template <size_t Dimension>
 void
-MeshBuilderBase::_dispatch()
+MeshBuilderBase::_dispatch(const DispatchType dispatch_type)
 {
-  if (parallel::size() == 1) {
-    return;
-  }
-
   using ConnectivityType = Connectivity<Dimension>;
   using Rd               = TinyVector<Dimension>;
   using MeshType         = Mesh<Dimension>;
 
-  if (not m_mesh) {
-    ConnectivityDescriptor descriptor;
-    std::shared_ptr connectivity = ConnectivityType::build(descriptor);
-    NodeValue<Rd> xr;
-    m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(connectivity, xr));
-  }
+  if ((parallel::size() == 1) and (dispatch_type == DispatchType::initial)) {
+    const MeshType& mesh = *(m_mesh->get<const MeshType>());
+
+    // force "creation" of a new mesh to avoid different
+    // parallel/sequential behaviors: if the mesh should change in
+    // parallel, is also changes in sequential.
+    m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(mesh.shared_connectivity(), mesh.xr()));
+  } else {
+    Assert((parallel::size() >= 1) or (dispatch_type == DispatchType::balance));
+
+    if (not m_mesh) {
+      ConnectivityDescriptor descriptor;
+      std::shared_ptr connectivity = ConnectivityType::build(descriptor);
+      NodeValue<Rd> xr;
+      m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(connectivity, xr));
+    }
 
-  const MeshType& mesh = *(m_mesh->get<const MeshType>());
-  ConnectivityDispatcher<Dimension> dispatcher(mesh.connectivity());
+    const MeshType& mesh = *(m_mesh->get<const MeshType>());
+
+    auto p_dispatcher = std::make_shared<const ConnectivityDispatcher<Dimension>>(mesh.connectivity());
 
-  std::shared_ptr dispatched_connectivity = dispatcher.dispatchedConnectivity();
-  NodeValue<Rd> dispatched_xr             = dispatcher.dispatch(mesh.xr());
+    m_connectivity_dispatcher = std::make_shared<ConnectivityDispatcherVariant>(p_dispatcher);
 
-  m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(dispatched_connectivity, dispatched_xr));
+    std::shared_ptr dispatched_connectivity = p_dispatcher->dispatchedConnectivity();
+    NodeValue<Rd> dispatched_xr             = p_dispatcher->dispatch(mesh.xr());
+
+    m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(dispatched_connectivity, dispatched_xr));
+  }
 }
 
-template void MeshBuilderBase::_dispatch<1>();
-template void MeshBuilderBase::_dispatch<2>();
-template void MeshBuilderBase::_dispatch<3>();
+template void MeshBuilderBase::_dispatch<1>(const DispatchType);
+template void MeshBuilderBase::_dispatch<2>(const DispatchType);
+template void MeshBuilderBase::_dispatch<3>(const DispatchType);
 
 template <size_t Dimension>
 void
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
index 665f0a88cc5a7f766ceb792594c6f5ace23529d6..a64105f3a5933630ee2a557a9484c05b753b05d8 100644
--- a/src/mesh/MeshBuilderBase.hpp
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -2,29 +2,44 @@
 #define MESH_BUILDER_BASE_HPP
 
 class MeshVariant;
+class ConnectivityDispatcherVariant;
 
 #include <memory>
 
 class MeshBuilderBase
 {
+ public:
+  enum class DispatchType : uint8_t
+  {
+    initial,
+    balance
+  };
+
  protected:
   std::shared_ptr<const MeshVariant> m_mesh;
+  std::shared_ptr<const ConnectivityDispatcherVariant> m_connectivity_dispatcher;
 
   template <size_t Dimension>
-  void _dispatch();
+  void _dispatch(const DispatchType balance);
 
   template <size_t Dimension>
   void _checkMesh() const;
 
  public:
+  std::shared_ptr<const ConnectivityDispatcherVariant>
+  connectivityDispatcher() const
+  {
+    return m_connectivity_dispatcher;
+  }
+
   std::shared_ptr<const MeshVariant>
   mesh() const
   {
     return m_mesh;
   }
 
-  MeshBuilderBase()  = default;
-  ~MeshBuilderBase() = default;
+  MeshBuilderBase()          = default;
+  virtual ~MeshBuilderBase() = default;
 };
 
 #endif   // MESH_BUILDER_BASE_HPP
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 64d35434462e5980a11b69ea2257ce779b3ec581..03e6337f572e494e94d1a31b48991c9444c881bd 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -45,7 +45,9 @@ class CharArrayEmbedder : public ICharArrayEmbedder
   void
   write(std::ostream& os) const final
   {
-    os.write(&(m_char_cast_array[0]), m_char_cast_array.size());
+    if (m_char_cast_array.size() > 0) {
+      os.write(&(m_char_cast_array[0]), m_char_cast_array.size());
+    }
   }
 
   CharArrayEmbedder(Array<InputDataT> array) : m_char_cast_array{array} {}
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index dadbe73bb24ad84e491248b89016baf85f0dbc31..a4642f38ffbc8a07e6e341ada1ce63e54a2137cf 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -3,13 +3,14 @@
 add_library(
   PugsScheme
   AcousticSolver.cpp
-  HyperelasticSolver.cpp
   DiscreteFunctionIntegrator.cpp
   DiscreteFunctionInterpoler.cpp
   DiscreteFunctionUtils.cpp
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   FluxingAdvectionSolver.cpp
+  HyperelasticSolver.cpp
+  LoadBalancer.cpp
 )
 
 target_link_libraries(
diff --git a/src/scheme/LoadBalancer.cpp b/src/scheme/LoadBalancer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ae49d645316d20331f16852894c58d10f3b350ff
--- /dev/null
+++ b/src/scheme/LoadBalancer.cpp
@@ -0,0 +1,52 @@
+#include <scheme/LoadBalancer.hpp>
+
+#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ConnectivityDispatcherVariant.hpp>
+#include <mesh/MeshBalancer.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <utils/Messenger.hpp>
+
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
+LoadBalancer::balance(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_list)
+{
+  if (not hasSameMesh(discrete_function_list)) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+  auto mesh_v = getCommonMesh(discrete_function_list);
+  Assert(mesh_v.use_count() > 0);
+
+  std::vector<std::shared_ptr<const DiscreteFunctionVariant>> balanced_discrete_function_list;
+
+  MeshBalancer mesh_balancer(mesh_v);
+
+  auto p_balanced_mesh_v = mesh_balancer.mesh();
+
+  std::visit(
+    [&mesh_balancer, &balanced_discrete_function_list, &discrete_function_list](auto&& p_balanced_mesh) {
+      using MeshType             = mesh_type_t<decltype(p_balanced_mesh)>;
+      constexpr size_t Dimension = MeshType::Dimension;
+      const auto& dispatcher     = mesh_balancer.connectivityDispatcher()->get<Dimension>();
+
+      for (size_t i_discrete_function = 0; i_discrete_function < discrete_function_list.size(); ++i_discrete_function) {
+        std::visit(
+          [&balanced_discrete_function_list, &dispatcher, &p_balanced_mesh](auto&& discrete_function) {
+            using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+            if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+              const auto& dispatched_cell_value = dispatcher->dispatch(discrete_function.cellValues());
+              balanced_discrete_function_list.push_back(
+                std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionT{p_balanced_mesh, dispatched_cell_value}));
+            } else {
+              const auto& dispatched_cell_array = dispatcher->dispatch(discrete_function.cellArrays());
+              balanced_discrete_function_list.push_back(
+                std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionT{p_balanced_mesh, dispatched_cell_array}));
+            }
+          },
+          discrete_function_list[i_discrete_function]->discreteFunction());
+      }
+    },
+    p_balanced_mesh_v->variant());
+
+  return balanced_discrete_function_list;
+}
diff --git a/src/scheme/LoadBalancer.hpp b/src/scheme/LoadBalancer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..465db1bb659b09745aae61c477635b7b5a28ebdb
--- /dev/null
+++ b/src/scheme/LoadBalancer.hpp
@@ -0,0 +1,21 @@
+#ifndef LOAD_BALANCER_HPP
+#define LOAD_BALANCER_HPP
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+class MeshVariant;
+class DiscreteFunctionVariant;
+
+class LoadBalancer
+{
+ public:
+  std::vector<std::shared_ptr<const DiscreteFunctionVariant>> balance(
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_list);
+
+  LoadBalancer()  = default;
+  ~LoadBalancer() = default;
+};
+
+#endif   // LOAD_BALANCER_HPP
diff --git a/src/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 2413bdbaeb5712e3b7eaf96fd9d09df5e8f0a6e0..8cdb5fc83a61e4c18ed2df1a24e3ec5d52fe7bdf 100644
--- a/src/utils/BuildInfo.cpp
+++ b/src/utils/BuildInfo.cpp
@@ -9,6 +9,14 @@
 #include <mpi.h>
 #endif   //  PUGS_HAS_MPI
 
+#ifdef PUGS_HAS_PARMETIS
+#include <parmetis.h>
+#endif   //  PUGS_HAS_PARMETIS
+
+#ifdef PUGS_HAS_PTSCOTCH
+#include <ptscotch.h>
+#endif   //  PUGS_HAS_PTSCOTCH
+
 #ifdef PUGS_HAS_PETSC
 #include <petsc.h>
 #endif   // PUGS_HAS_PETSC
@@ -64,6 +72,27 @@ BuildInfo::mpiLibrary()
 #endif   // PUGS_HAS_MPI
 }
 
+std::string
+BuildInfo::parmetisLibrary()
+{
+#ifdef PUGS_HAS_PARMETIS
+  return stringify(PARMETIS_MAJOR_VERSION) + "." + stringify(PARMETIS_MINOR_VERSION) + "." +
+         stringify(PARMETIS_SUBMINOR_VERSION);
+#else    // PUGS_HAS_PARMETIS
+  return "none";
+#endif   // PUGS_HAS_PARMETIS
+}
+
+std::string
+BuildInfo::ptscotchLibrary()
+{
+#ifdef PUGS_HAS_PTSCOTCH
+  return stringify(SCOTCH_VERSION) + "." + stringify(SCOTCH_RELEASE) + "." + stringify(SCOTCH_PATCHLEVEL);
+#else    // PUGS_HAS_PTSCOTCH
+  return "none";
+#endif   // PUGS_HAS_PTSCOTCH
+}
+
 std::string
 BuildInfo::petscLibrary()
 {
diff --git a/src/utils/BuildInfo.hpp b/src/utils/BuildInfo.hpp
index 86f7c04f035827a936a00c9b1f42e304e33e7b26..0142b5d088e6323858116c7949cb7b0b681e2acc 100644
--- a/src/utils/BuildInfo.hpp
+++ b/src/utils/BuildInfo.hpp
@@ -9,6 +9,8 @@ struct BuildInfo
   static std::string compiler();
   static std::string kokkosDevices();
   static std::string mpiLibrary();
+  static std::string parmetisLibrary();
+  static std::string ptscotchLibrary();
   static std::string petscLibrary();
   static std::string slepcLibrary();
   static std::string eigen3Library();
diff --git a/src/utils/ParMETISPartitioner.cpp b/src/utils/ParMETISPartitioner.cpp
index f1b17c3d070a71e0290c8c6434b0e13d2e442bad..5a7bdb6a27d1ce459d19a4cfac38c2cca2166d64 100644
--- a/src/utils/ParMETISPartitioner.cpp
+++ b/src/utils/ParMETISPartitioner.cpp
@@ -53,8 +53,23 @@ ParMETISPartitioner::partition(const CRSGraph& graph)
 
     int local_number_of_nodes = graph.numberOfNodes();
 
+    int group_size;
+    MPI_Group_size(mesh_group, &group_size);
+
+    std::vector<int> local_number_of_nodes_list;
+    local_number_of_nodes_list.resize(group_size);
+    MPI_Datatype mpi_int = parallel::Messenger::helper::mpiType<int>();
+    MPI_Allgather(&local_number_of_nodes, 1, mpi_int,             //
+                  &(local_number_of_nodes_list[0]), 1, mpi_int,   //
+                  partitioning_comm);
+
+    std::vector<int> vtxdist;
+    vtxdist.push_back(0);
+    for (size_t i = 0; i < local_number_of_nodes_list.size(); ++i) {
+      vtxdist.push_back(vtxdist[i] + local_number_of_nodes_list[i]);
+    }
+
     partition = Array<int>(local_number_of_nodes);
-    std::vector<int> vtxdist{0, local_number_of_nodes};
 
     const Array<const int>& entries   = graph.entries();
     const Array<const int>& neighbors = graph.neighbors();
diff --git a/src/utils/Partitioner.cpp b/src/utils/Partitioner.cpp
index 4f53d322f5e52c130122ccc3966270e7220d34d7..9ef628245f73b3be4f6ecd875ae5dd255821027b 100644
--- a/src/utils/Partitioner.cpp
+++ b/src/utils/Partitioner.cpp
@@ -15,8 +15,10 @@ Partitioner::partition(const CRSGraph& graph)
   case PartitionerLibrary::ptscotch: {
     return PTScotchPartitioner::partition(graph);
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid partition library");
   }
+    // LCOV_EXCL_STOP
   }
 }
diff --git a/src/utils/PartitionerOptions.hpp b/src/utils/PartitionerOptions.hpp
index e4157f7659cac94ef8b4c7cd02a54d289156d4a8..7adce0db4ae35ace3a5d75509efea5742e2957fb 100644
--- a/src/utils/PartitionerOptions.hpp
+++ b/src/utils/PartitionerOptions.hpp
@@ -29,7 +29,7 @@ name(const PartitionerLibrary library)
   case PartitionerLibrary::PT__end: {
   }
   }
-  throw UnexpectedError("Linear system library name is not defined!");
+  throw UnexpectedError("Partitioner library name is not defined!");
 }
 
 template <typename PartitionerEnumType>
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index ac64a50cfe0b026fe0185d79e30d4f6d648317ac..356869699a7d76f04648dce9a95d4d091a619df2 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -62,6 +62,8 @@ pugsBuildInfo()
   os << "compiler: " << rang::style::bold << BuildInfo::compiler() << rang::style::reset << '\n';
   os << "kokkos:   " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
   os << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
+  os << "ParMETIS: " << rang::style::bold << BuildInfo::parmetisLibrary() << rang::style::reset << '\n';
+  os << "PTScotch: " << rang::style::bold << BuildInfo::ptscotchLibrary() << rang::style::reset << '\n';
   os << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
   os << "SLEPc:    " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
   os << "Eigen3:   " << rang::style::bold << BuildInfo::eigen3Library() << rang::style::reset << '\n';
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 1dd13e6e3ce1f23c31e046db7597a1ed2ffb947d..169d375538eae053ff85a70919e11c24620d8c8a 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -134,6 +134,7 @@ add_executable (unit_tests
   test_LinearSolverOptions.cpp
   test_LineTransformation.cpp
   test_ListAffectationProcessor.cpp
+  test_LoadBalancer_serial.cpp
   test_MathModule.cpp
   test_MedianDualConnectivityBuilder.cpp
   test_MedianDualMeshBuilder.cpp
@@ -142,6 +143,7 @@ add_executable (unit_tests
   test_OStream.cpp
   test_ParallelChecker_write.cpp
   test_ParseError.cpp
+  test_PartitionerOptions.cpp
   test_PETScUtils.cpp
   test_PrimalToDiamondDualConnectivityDataMapper.cpp
   test_PrimalToDual1DConnectivityDataMapper.cpp
@@ -254,6 +256,7 @@ add_executable (mpi_unit_tests
   test_ItemValueUtils.cpp
   test_ItemValueVariant.cpp
   test_ItemValueVariantFunctionInterpoler.cpp
+  test_LoadBalancer_parallel.cpp
   test_MeshEdgeBoundary.cpp
   test_MeshEdgeInterface.cpp
   test_MeshFaceBoundary.cpp
diff --git a/tests/test_BuildInfo.cpp b/tests/test_BuildInfo.cpp
index 9414902f4d63b498e68bffd4d84b3cae32a165bd..02cf31adc012a1d900d5540838a946e828ac301c 100644
--- a/tests/test_BuildInfo.cpp
+++ b/tests/test_BuildInfo.cpp
@@ -12,6 +12,14 @@
 #include <mpi.h>
 #endif   //  PUGS_HAS_MPI
 
+#ifdef PUGS_HAS_PARMETIS
+#include <parmetis.h>
+#endif   //  PUGS_HAS_PARMETIS
+
+#ifdef PUGS_HAS_PTSCOTCH
+#include <ptscotch.h>
+#endif   //  PUGS_HAS_PTSCOTCH
+
 #ifdef PUGS_HAS_PETSC
 #include <petsc.h>
 #endif   // PUGS_HAS_PETSC
@@ -57,6 +65,30 @@ TEST_CASE("BuildInfo", "[utils]")
 #endif   // PUGS_HAS_MPI
   }
 
+  SECTION("ParMETIS")
+  {
+#ifdef PUGS_HAS_PARMETIS
+    const std::string parmetis_library = stringify(PARMETIS_MAJOR_VERSION) + "." + stringify(PARMETIS_MINOR_VERSION) +
+                                         "." + stringify(PARMETIS_SUBMINOR_VERSION);
+
+    REQUIRE(BuildInfo::parmetisLibrary() == parmetis_library);
+#else
+    REQUIRE(BuildInfo::parmetisLibrary() == "none");
+#endif   // PUGS_HAS_PARMETIS
+  }
+
+  SECTION("PTScotch")
+  {
+#ifdef PUGS_HAS_PTSCOTCH
+    const std::string ptscotch_library =
+      stringify(SCOTCH_VERSION) + "." + stringify(SCOTCH_RELEASE) + "." + stringify(SCOTCH_PATCHLEVEL);
+
+    REQUIRE(BuildInfo::ptscotchLibrary() == ptscotch_library);
+#else
+    REQUIRE(BuildInfo::ptscotchLibrary() == "none");
+#endif   // PUGS_HAS_PTSCOTCH
+  }
+
   SECTION("petsc")
   {
 #ifdef PUGS_HAS_PETSC
diff --git a/tests/test_LoadBalancer_parallel.cpp b/tests/test_LoadBalancer_parallel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a252c74e0705a8cc2d1b34e54d0cae167c1ea674
--- /dev/null
+++ b/tests/test_LoadBalancer_parallel.cpp
@@ -0,0 +1,271 @@
+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <scheme/LoadBalancer.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <utils/PartitionerOptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LoadBalancer parallel", "[scheme]")
+{
+  auto pugs_has_partitioner = [&](const PartitionerLibrary& partitioner_library) {
+    switch (partitioner_library) {
+    case PartitionerLibrary::parmetis: {
+#ifdef PUGS_HAS_PARMETIS
+      return true;
+#else    // PUGS_HAS_PARMETIS
+      return false;
+#endif   // PUGS_HAS_PARMETIS
+    }
+    case PartitionerLibrary::ptscotch: {
+#ifdef PUGS_HAS_PARMETIS
+      return true;
+#else    // PUGS_HAS_PARMETIS
+      return false;
+#endif   // PUGS_HAS_PARMETIS
+    }
+    default: {
+      return false;
+    }
+    }
+  };
+
+  constexpr std::array partitioner_libs = {PartitionerLibrary::parmetis, PartitionerLibrary::ptscotch};
+
+  const auto initial_library = PartitionerOptions::default_options.library();
+
+  for (auto lib : partitioner_libs) {
+    if (pugs_has_partitioner(lib)) {
+      SECTION(name(lib))
+      {
+        SECTION("1D")
+        {
+          for (auto&& named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+            std::shared_ptr<const Mesh<1>> initial_mesh_1d = named_mesh.mesh()->get<Mesh<1>>();
+            SECTION(named_mesh.name())
+            {
+              DiscreteFunctionP0<double> initial_fh{initial_mesh_1d};
+              DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_1d, 3};
+
+              {
+                auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_1d).xj();
+
+                for (CellId cell_id = 0; cell_id < initial_mesh_1d->numberOfCells(); ++cell_id) {
+                  initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+                  initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+                  initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0]);
+                  initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] + 4;
+                }
+              }
+
+              LoadBalancer load_balancer;
+              auto balanced_discrete_function_list =
+                load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                       std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+              REQUIRE(balanced_discrete_function_list.size() == 2);
+
+              DiscreteFunctionP0<const double> final_fh =
+                balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+              DiscreteFunctionP0Vector<const double> final_vh =
+                balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+              REQUIRE(final_vh.size() == 3);
+              {
+                auto final_mesh_1d_v = final_fh.meshVariant();
+                auto final_mesh_1d   = final_mesh_1d_v->get<Mesh<1>>();
+
+                REQUIRE(final_mesh_1d->id() > initial_mesh_1d->id());
+
+                REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                           std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                                    "error: discrete functions are not defined on the same mesh");
+
+                {
+                  auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_1d).xj();
+
+                  bool final_fh_is_same = true;
+                  bool final_vh_is_same = true;
+                  for (CellId cell_id = 0; cell_id < final_mesh_1d->numberOfCells(); ++cell_id) {
+                    if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                      final_fh_is_same = false;
+                    }
+
+                    if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0])) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] + 4) {
+                      final_vh_is_same = false;
+                    }
+                  }
+
+                  REQUIRE(parallel::allReduceAnd(final_fh_is_same));
+                  REQUIRE(parallel::allReduceAnd(final_vh_is_same));
+                }
+              }
+            }
+          }
+        }
+
+        SECTION("2D")
+        {
+          for (auto&& named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+            std::shared_ptr<const Mesh<2>> initial_mesh_2d = named_mesh.mesh()->get<Mesh<2>>();
+            SECTION(named_mesh.name())
+            {
+              DiscreteFunctionP0<double> initial_fh{initial_mesh_2d};
+              DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_2d, 3};
+
+              {
+                auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_2d).xj();
+
+                for (CellId cell_id = 0; cell_id < initial_mesh_2d->numberOfCells(); ++cell_id) {
+                  initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+                  initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+                  initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0] + initial_xj[cell_id][1]);
+                  initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+                }
+              }
+
+              LoadBalancer load_balancer;
+              auto balanced_discrete_function_list =
+                load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                       std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+              REQUIRE(balanced_discrete_function_list.size() == 2);
+
+              DiscreteFunctionP0<const double> final_fh =
+                balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+              DiscreteFunctionP0Vector<const double> final_vh =
+                balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+              REQUIRE(final_vh.size() == 3);
+              {
+                auto final_mesh_2d_v = final_fh.meshVariant();
+                auto final_mesh_2d   = final_mesh_2d_v->get<Mesh<2>>();
+
+                REQUIRE(final_mesh_2d->id() > initial_mesh_2d->id());
+
+                REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                           std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                                    "error: discrete functions are not defined on the same mesh");
+
+                {
+                  auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_2d).xj();
+
+                  bool final_fh_is_same = true;
+                  bool final_vh_is_same = true;
+                  for (CellId cell_id = 0; cell_id < final_mesh_2d->numberOfCells(); ++cell_id) {
+                    if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                      final_fh_is_same = false;
+                    }
+
+                    if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0] + final_xj[cell_id][1])) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                      final_vh_is_same = false;
+                    }
+                  }
+
+                  REQUIRE(parallel::allReduceAnd(final_fh_is_same));
+                  REQUIRE(parallel::allReduceAnd(final_vh_is_same));
+                }
+              }
+            }
+          }
+        }
+
+        SECTION("3D")
+        {
+          for (auto&& named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+            std::shared_ptr<const Mesh<3>> initial_mesh_3d = named_mesh.mesh()->get<Mesh<3>>();
+            SECTION(named_mesh.name())
+            {
+              DiscreteFunctionP0<double> initial_fh{initial_mesh_3d};
+              DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_3d, 3};
+
+              {
+                auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_3d).xj();
+
+                for (CellId cell_id = 0; cell_id < initial_mesh_3d->numberOfCells(); ++cell_id) {
+                  initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+                  initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+                  initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][1] + initial_xj[cell_id][2]);
+                  initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+                }
+              }
+
+              LoadBalancer load_balancer;
+              auto balanced_discrete_function_list =
+                load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                       std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+              REQUIRE(balanced_discrete_function_list.size() == 2);
+
+              DiscreteFunctionP0<const double> final_fh =
+                balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+              DiscreteFunctionP0Vector<const double> final_vh =
+                balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+              REQUIRE(final_vh.size() == 3);
+              {
+                auto final_mesh_3d_v = final_fh.meshVariant();
+                auto final_mesh_3d   = final_mesh_3d_v->get<Mesh<3>>();
+
+                REQUIRE(final_mesh_3d->id() > initial_mesh_3d->id());
+
+                REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                           std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                                    "error: discrete functions are not defined on the same mesh");
+
+                {
+                  auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_3d).xj();
+
+                  bool final_fh_is_same = true;
+                  bool final_vh_is_same = true;
+                  for (CellId cell_id = 0; cell_id < final_mesh_3d->numberOfCells(); ++cell_id) {
+                    if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                      final_fh_is_same = false;
+                    }
+
+                    if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][1] + final_xj[cell_id][2])) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                      final_vh_is_same = false;
+                    }
+                  }
+
+                  REQUIRE(parallel::allReduceAnd(final_fh_is_same));
+                  REQUIRE(parallel::allReduceAnd(final_vh_is_same));
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  PartitionerOptions::default_options.library() = initial_library;
+}
diff --git a/tests/test_LoadBalancer_serial.cpp b/tests/test_LoadBalancer_serial.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..63f6271770983c7a5dd8b6da82f3c39be8cfd2c5
--- /dev/null
+++ b/tests/test_LoadBalancer_serial.cpp
@@ -0,0 +1,235 @@
+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <scheme/LoadBalancer.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LoadBalancer serial", "[scheme]")
+{
+  SECTION("1D")
+  {
+    for (auto&& named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+      std::shared_ptr<const Mesh<1>> initial_mesh_1d = named_mesh.mesh()->get<Mesh<1>>();
+      SECTION(named_mesh.name())
+      {
+        DiscreteFunctionP0<double> initial_fh{initial_mesh_1d};
+        DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_1d, 3};
+
+        {
+          auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_1d).xj();
+
+          for (CellId cell_id = 0; cell_id < initial_mesh_1d->numberOfCells(); ++cell_id) {
+            initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+            initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+            initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0]);
+            initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] + 4;
+          }
+        }
+
+        LoadBalancer load_balancer;
+        auto balanced_discrete_function_list =
+          load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                 std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+        REQUIRE(balanced_discrete_function_list.size() == 2);
+
+        DiscreteFunctionP0<const double> final_fh =
+          balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+        DiscreteFunctionP0Vector<const double> final_vh =
+          balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+        REQUIRE(final_vh.size() == 3);
+        {
+          auto final_mesh_1d_v = final_fh.meshVariant();
+          auto final_mesh_1d   = final_mesh_1d_v->get<Mesh<1>>();
+
+          REQUIRE(final_mesh_1d->id() > initial_mesh_1d->id());
+
+          REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                     std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                              "error: discrete functions are not defined on the same mesh");
+
+          {
+            auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_1d).xj();
+
+            bool final_fh_is_same = true;
+            bool final_vh_is_same = true;
+            for (CellId cell_id = 0; cell_id < final_mesh_1d->numberOfCells(); ++cell_id) {
+              if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                final_fh_is_same = false;
+              }
+
+              if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0])) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] + 4) {
+                final_vh_is_same = false;
+              }
+            }
+
+            REQUIRE(final_fh_is_same);
+            REQUIRE(final_vh_is_same);
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    for (auto&& named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+      std::shared_ptr<const Mesh<2>> initial_mesh_2d = named_mesh.mesh()->get<Mesh<2>>();
+      SECTION(named_mesh.name())
+      {
+        DiscreteFunctionP0<double> initial_fh{initial_mesh_2d};
+        DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_2d, 3};
+
+        {
+          auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_2d).xj();
+
+          for (CellId cell_id = 0; cell_id < initial_mesh_2d->numberOfCells(); ++cell_id) {
+            initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+            initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+            initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0] + initial_xj[cell_id][1]);
+            initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+          }
+        }
+
+        LoadBalancer load_balancer;
+        auto balanced_discrete_function_list =
+          load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                 std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+        REQUIRE(balanced_discrete_function_list.size() == 2);
+
+        DiscreteFunctionP0<const double> final_fh =
+          balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+        DiscreteFunctionP0Vector<const double> final_vh =
+          balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+        REQUIRE(final_vh.size() == 3);
+        {
+          auto final_mesh_2d_v = final_fh.meshVariant();
+          auto final_mesh_2d   = final_mesh_2d_v->get<Mesh<2>>();
+
+          REQUIRE(final_mesh_2d->id() > initial_mesh_2d->id());
+
+          REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                     std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                              "error: discrete functions are not defined on the same mesh");
+
+          {
+            auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_2d).xj();
+
+            bool final_fh_is_same = true;
+            bool final_vh_is_same = true;
+            for (CellId cell_id = 0; cell_id < final_mesh_2d->numberOfCells(); ++cell_id) {
+              if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                final_fh_is_same = false;
+              }
+
+              if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0] + final_xj[cell_id][1])) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                final_vh_is_same = false;
+              }
+            }
+
+            REQUIRE(final_fh_is_same);
+            REQUIRE(final_vh_is_same);
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    for (auto&& named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+      std::shared_ptr<const Mesh<3>> initial_mesh_3d = named_mesh.mesh()->get<Mesh<3>>();
+      SECTION(named_mesh.name())
+      {
+        DiscreteFunctionP0<double> initial_fh{initial_mesh_3d};
+        DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_3d, 3};
+
+        {
+          auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_3d).xj();
+
+          for (CellId cell_id = 0; cell_id < initial_mesh_3d->numberOfCells(); ++cell_id) {
+            initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+            initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+            initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][1] + initial_xj[cell_id][2]);
+            initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+          }
+        }
+
+        LoadBalancer load_balancer;
+        auto balanced_discrete_function_list =
+          load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                 std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+        REQUIRE(balanced_discrete_function_list.size() == 2);
+
+        DiscreteFunctionP0<const double> final_fh =
+          balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+        DiscreteFunctionP0Vector<const double> final_vh =
+          balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+        REQUIRE(final_vh.size() == 3);
+        {
+          auto final_mesh_3d_v = final_fh.meshVariant();
+          auto final_mesh_3d   = final_mesh_3d_v->get<Mesh<3>>();
+
+          REQUIRE(final_mesh_3d->id() > initial_mesh_3d->id());
+
+          REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                     std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                              "error: discrete functions are not defined on the same mesh");
+
+          {
+            auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_3d).xj();
+
+            bool final_fh_is_same = true;
+            bool final_vh_is_same = true;
+            for (CellId cell_id = 0; cell_id < final_mesh_3d->numberOfCells(); ++cell_id) {
+              if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                final_fh_is_same = false;
+              }
+
+              if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][1] + final_xj[cell_id][2])) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                final_vh_is_same = false;
+              }
+            }
+
+            REQUIRE(final_fh_is_same);
+            REQUIRE(final_vh_is_same);
+          }
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_Partitioner.cpp b/tests/test_Partitioner.cpp
index 78a427f6d31fd35f27519498a0382e681f3421c4..67ddbb11161d5ab95a065c5f1f7b6bb0dbfb1ad3 100644
--- a/tests/test_Partitioner.cpp
+++ b/tests/test_Partitioner.cpp
@@ -11,8 +11,6 @@ TEST_CASE("Partitioner", "[utils]")
 {
   SECTION("one graph split to all")
   {
-    Partitioner partitioner;
-
     std::vector<int> entries_vector;
     std::vector<int> neighbors_vector;
 
@@ -70,17 +68,54 @@ TEST_CASE("Partitioner", "[utils]")
 
     CRSGraph graph{entries, neighbors};
 
-    Array<int> partitioned = partitioner.partition(graph);
+    SECTION("ParMETIS")
+    {
+      PartitionerOptions partitioner_options;
+      partitioner_options.library() = PartitionerLibrary::parmetis;
+      Partitioner partitioner{partitioner_options};
 
-    REQUIRE((partitioned.size() + 1) == entries.size());
+      Array<int> partitioned = partitioner.partition(graph);
 
-    if (parallel::rank() == 0) {
-      std::set<int> assigned_ranks;
-      for (size_t i = 0; i < partitioned.size(); ++i) {
-        assigned_ranks.insert(partitioned[i]);
+      REQUIRE((partitioned.size() + 1) == entries.size());
+
+#ifdef PUGS_HAS_PARMETIS
+      if (parallel::rank() == 0) {
+        std::set<int> assigned_ranks;
+        for (size_t i = 0; i < partitioned.size(); ++i) {
+          assigned_ranks.insert(partitioned[i]);
+        }
+
+        REQUIRE(assigned_ranks.size() == parallel::size());
       }
+#else    // PUGS_HAS_PARMETIS
+      REQUIRE(min(partitioned) == 0);
+      REQUIRE(max(partitioned) == 0);
+#endif   // PUGS_HAS_PARMETIS
+    }
 
-      REQUIRE(assigned_ranks.size() == parallel::size());
+    SECTION("PTScotch")
+    {
+      PartitionerOptions partitioner_options;
+      partitioner_options.library() = PartitionerLibrary::ptscotch;
+      Partitioner partitioner{partitioner_options};
+
+      Array<int> partitioned = partitioner.partition(graph);
+
+      REQUIRE((partitioned.size() + 1) == entries.size());
+
+#ifdef PUGS_HAS_PTSCOTCH
+      if (parallel::rank() == 0) {
+        std::set<int> assigned_ranks;
+        for (size_t i = 0; i < partitioned.size(); ++i) {
+          assigned_ranks.insert(partitioned[i]);
+        }
+
+        REQUIRE(assigned_ranks.size() == parallel::size());
+      }
+#else    // PUGS_HAS_PTSCOTCH
+      REQUIRE(min(partitioned) == 0);
+      REQUIRE(max(partitioned) == 0);
+#endif   // PUGS_HAS_PTSCOTCH
     }
   }
 }
diff --git a/tests/test_PartitionerOptions.cpp b/tests/test_PartitionerOptions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d1de1dde0e010fab8f54126783631f4a7cf9c499
--- /dev/null
+++ b/tests/test_PartitionerOptions.cpp
@@ -0,0 +1,45 @@
+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <utils/PartitionerOptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PartitionerOptions", "[utils]")
+{
+  SECTION("name")
+  {
+    REQUIRE(name(PartitionerLibrary::parmetis) == "ParMETIS");
+    REQUIRE(name(PartitionerLibrary::ptscotch) == "PTScotch");
+    REQUIRE_THROWS_WITH(name(PartitionerLibrary::PT__end),
+                        "unexpected error: Partitioner library name is not defined!");
+  }
+
+  SECTION("type from name")
+  {
+    REQUIRE(getPartitionerEnumFromName<PartitionerLibrary>("ParMETIS") == PartitionerLibrary::parmetis);
+    REQUIRE(getPartitionerEnumFromName<PartitionerLibrary>("PTScotch") == PartitionerLibrary::ptscotch);
+
+    REQUIRE_THROWS_WITH(getPartitionerEnumFromName<PartitionerLibrary>("foobar"),
+                        "error: could not find 'foobar' associate type!");
+  }
+
+  SECTION("output")
+  {
+    {
+      PartitionerOptions options;
+      options.library() = PartitionerLibrary::parmetis;
+      std::stringstream os;
+      os << options;
+      REQUIRE(os.str() == "  library: ParMETIS\n");
+    }
+
+    {
+      PartitionerOptions options;
+      options.library() = PartitionerLibrary::ptscotch;
+      std::stringstream os;
+      os << options;
+      REQUIRE(os.str() == "  library: PTScotch\n");
+    }
+  }
+}
diff --git a/tests/test_PugsUtils.cpp b/tests/test_PugsUtils.cpp
index 5ba8817fcaa34ddf31f0955ed9fd118f841bd87d..175ed929467e17bc1fd29e2279f60d3a04b1a1c9 100644
--- a/tests/test_PugsUtils.cpp
+++ b/tests/test_PugsUtils.cpp
@@ -49,6 +49,8 @@ TEST_CASE("PugsUtils", "[utils]")
       os << "compiler: " << rang::style::bold << BuildInfo::compiler() << rang::style::reset << '\n';
       os << "kokkos:   " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
       os << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
+      os << "ParMETIS: " << rang::style::bold << BuildInfo::parmetisLibrary() << rang::style::reset << '\n';
+      os << "PTScotch: " << rang::style::bold << BuildInfo::ptscotchLibrary() << rang::style::reset << '\n';
       os << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
       os << "SLEPc:    " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
       os << "Eigen3:   " << rang::style::bold << BuildInfo::eigen3Library() << rang::style::reset << '\n';
diff --git a/tests/test_checkpointing_Checkpoint.cpp b/tests/test_checkpointing_Checkpoint.cpp
index 317683b26dfca317b4b6f0059b6a8e55c8d5e4f0..9436797a726b303698156fc0f5991c29d0f146e9 100644
--- a/tests/test_checkpointing_Checkpoint.cpp
+++ b/tests/test_checkpointing_Checkpoint.cpp
@@ -152,13 +152,13 @@ for(let i:N, i=0; i<3; ++i) {
 
     HighFive::Group m = embedded1.getGroup("m");
     REQUIRE(m.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+    REQUIRE(m.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
     HighFive::Group singleton        = checkpoint.getGroup("singleton");
     HighFive::Group global_variables = singleton.getGroup("global_variables");
     REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() ==
             initial_connectivity_id + 1 + (parallel::size() > 1));
-    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 1 + (parallel::size() > 1));
+    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 2);
     HighFive::Group execution_info = singleton.getGroup("execution_info");
     REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 1);
 
@@ -170,10 +170,10 @@ for(let i:N, i=0; i<3; ++i) {
     REQUIRE(connectivity0.getAttribute("type").read<std::string>() == "unstructured");
 
     HighFive::Group mesh  = checkpoint.getGroup("mesh");
-    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + (parallel::size() > 1)));
+    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
     REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + (parallel::size() > 1));
     REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
     HighFive::Group functions = checkpoint.getGroup("functions");
diff --git a/tests/test_checkpointing_Checkpoint_sequential.cpp b/tests/test_checkpointing_Checkpoint_sequential.cpp
index 160404f8835a7fcf6a0c94700ae79d023391b2ca..825138fabde43964790c900eb2ce0ad4338e73e5 100644
--- a/tests/test_checkpointing_Checkpoint_sequential.cpp
+++ b/tests/test_checkpointing_Checkpoint_sequential.cpp
@@ -158,19 +158,19 @@ for(let i:N, i=0; i<3; ++i) {
     REQUIRE(duals.getAttribute("type").read<std::string>() == "(mesh)");
     HighFive::Group duals_0 = duals.getGroup("0");
     REQUIRE(duals_0.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
+    REQUIRE(duals_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
     HighFive::Group duals_1 = duals.getGroup("1");
     REQUIRE(duals_1.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(duals_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
 
     HighFive::Group m = embedded1.getGroup("m");
     REQUIRE(m.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(m.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
     HighFive::Group singleton        = checkpoint.getGroup("singleton");
     HighFive::Group global_variables = singleton.getGroup("global_variables");
     REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() == initial_connectivity_id + 3);
-    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 4);
     HighFive::Group execution_info = singleton.getGroup("execution_info");
     REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 1);
 
@@ -193,21 +193,21 @@ for(let i:N, i=0; i<3; ++i) {
     REQUIRE(connectivity2.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Median);
 
     HighFive::Group mesh  = checkpoint.getGroup("mesh");
-    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id));
+    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
     REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id);
     REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
-    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
-    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
+    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh1.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh1.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
-    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
-    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
-    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
+    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh2.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh2.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Median);
 
diff --git a/tests/test_checkpointing_Connectivity.cpp b/tests/test_checkpointing_Connectivity.cpp
index dcb6bb14665bebec43b08eaa710e67ab3fc5f76d..01c4a6b70ac8defc28e4e0aa1bfb84e986a67201 100644
--- a/tests/test_checkpointing_Connectivity.cpp
+++ b/tests/test_checkpointing_Connectivity.cpp
@@ -57,7 +57,7 @@ TEST_CASE("checkpointing_Connectivity", "[utils/checkpointing]")
         auto new_connectivity_1d = test_only::duplicateConnectivity(mesh_1d->connectivity());
         checkpointing::writeConnectivity(*new_connectivity_1d, file, checkpoint_group_0);
         checkpointing::writeConnectivity(*new_connectivity_1d, file, checkpoint_group_1);
-        connectivity_id_map[mesh_1d->id()] = new_connectivity_1d->id();
+        connectivity_id_map[mesh_1d->connectivity().id()] = new_connectivity_1d->id();
 
         auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
 
@@ -69,7 +69,7 @@ TEST_CASE("checkpointing_Connectivity", "[utils/checkpointing]")
 
         checkpointing::writeConnectivity(*new_connectivity_2d, file, checkpoint_group_0);
         checkpointing::writeConnectivity(*new_connectivity_2d, file, checkpoint_group_1);
-        connectivity_id_map[mesh_2d->id()] = new_connectivity_2d->id();
+        connectivity_id_map[mesh_2d->connectivity().id()] = new_connectivity_2d->id();
 
         HighFive::Group global_variables_group_0 = checkpoint_group_0.createGroup("singleton/global_variables");
         global_variables_group_0.createAttribute("connectivity_id",
@@ -83,7 +83,7 @@ TEST_CASE("checkpointing_Connectivity", "[utils/checkpointing]")
 
         auto new_connectivity_3d = test_only::duplicateConnectivity(mesh_3d->connectivity());
         checkpointing::writeConnectivity(*new_connectivity_3d, file, checkpoint_group_1);
-        connectivity_id_map[mesh_3d->id()] = new_connectivity_3d->id();
+        connectivity_id_map[mesh_3d->connectivity().id()] = new_connectivity_3d->id();
 
         // creates artificially holes in numbering
         test_only::duplicateConnectivity(mesh_3d->connectivity());
diff --git a/tests/test_checkpointing_Resume.cpp b/tests/test_checkpointing_Resume.cpp
index 82ea9202b75958e7d77a33c7af47d9a4516860bc..13f1b341ec80e8ad654d8f03bb7305e2a9420f44 100644
--- a/tests/test_checkpointing_Resume.cpp
+++ b/tests/test_checkpointing_Resume.cpp
@@ -229,22 +229,21 @@ for(let i:N, i=0; i<3; ++i) {
 
       HighFive::Group m1d = embedded1.getGroup("m1d");
       REQUIRE(m1d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (1 + 2 * (parallel::size() > 1)));
+      REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
 
       HighFive::Group m2d = embedded1.getGroup("m2d");
       REQUIRE(m2d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+      REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
       HighFive::Group m3d = embedded1.getGroup("m3d");
       REQUIRE(m3d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (2 + 3 * (parallel::size() > 1)));
+      REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
 
       HighFive::Group singleton        = checkpoint.getGroup("singleton");
       HighFive::Group global_variables = singleton.getGroup("global_variables");
       REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() ==
               initial_connectivity_id + 3 * (1 + (parallel::size() > 1)));
-      REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() ==
-              initial_mesh_id + 3 * (1 + (parallel::size() > 1)));
+      REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 6);
       HighFive::Group execution_info = singleton.getGroup("execution_info");
       REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 1);
 
@@ -270,24 +269,24 @@ for(let i:N, i=0; i<3; ++i) {
       REQUIRE(connectivity2.getAttribute("type").read<std::string>() == "unstructured");
 
       HighFive::Group mesh  = checkpoint.getGroup("mesh");
-      HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + (parallel::size() > 1)));
+      HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
       REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + (parallel::size() > 1));
       REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-      REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+      REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
       REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
-      HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + (1 + 2 * (parallel::size() > 1))));
+      HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
       REQUIRE(mesh1.getAttribute("connectivity").read<uint64_t>() ==
               initial_connectivity_id + (1 + 2 * (parallel::size() > 1)));
       REQUIRE(mesh1.getAttribute("dimension").read<uint64_t>() == 1);
-      REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + (1 + 2 * (parallel::size() > 1)));
+      REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
       REQUIRE(mesh1.getAttribute("type").read<std::string>() == "polygonal");
 
-      HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + (2 + 3 * (parallel::size() > 1))));
+      HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
       REQUIRE(mesh2.getAttribute("connectivity").read<uint64_t>() ==
               initial_connectivity_id + (2 + 3 * (parallel::size() > 1)));
       REQUIRE(mesh2.getAttribute("dimension").read<uint64_t>() == 3);
-      REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + (2 + 3 * (parallel::size() > 1)));
+      REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
       REQUIRE(mesh2.getAttribute("type").read<std::string>() == "polygonal");
 
       HighFive::Group functions = checkpoint.getGroup("functions");
@@ -357,22 +356,21 @@ for(let i:N, i=0; i<3; ++i) {
 
       HighFive::Group m1d = embedded1.getGroup("m1d");
       REQUIRE(m1d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (1 + 2 * (parallel::size() > 1)));
+      REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
 
       HighFive::Group m2d = embedded1.getGroup("m2d");
       REQUIRE(m2d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+      REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
       HighFive::Group m3d = embedded1.getGroup("m3d");
       REQUIRE(m3d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (2 + 3 * (parallel::size() > 1)));
+      REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
 
       HighFive::Group singleton        = checkpoint.getGroup("singleton");
       HighFive::Group global_variables = singleton.getGroup("global_variables");
       REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() ==
               initial_connectivity_id + 3 * (1 + (parallel::size() > 1)));
-      REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() ==
-              initial_mesh_id + 3 * (1 + (parallel::size() > 1)));
+      REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 6);
       HighFive::Group execution_info = singleton.getGroup("execution_info");
       REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 1);
 
@@ -398,24 +396,24 @@ for(let i:N, i=0; i<3; ++i) {
       REQUIRE(connectivity2.getAttribute("type").read<std::string>() == "unstructured");
 
       HighFive::Group mesh  = checkpoint.getGroup("mesh");
-      HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + (parallel::size() > 1)));
+      HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
       REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + (parallel::size() > 1));
       REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-      REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+      REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
       REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
-      HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + (1 + 2 * (parallel::size() > 1))));
+      HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
       REQUIRE(mesh1.getAttribute("connectivity").read<uint64_t>() ==
               initial_connectivity_id + (1 + 2 * (parallel::size() > 1)));
       REQUIRE(mesh1.getAttribute("dimension").read<uint64_t>() == 1);
-      REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + (1 + 2 * (parallel::size() > 1)));
+      REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
       REQUIRE(mesh1.getAttribute("type").read<std::string>() == "polygonal");
 
-      HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + (2 + 3 * (parallel::size() > 1))));
+      HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
       REQUIRE(mesh2.getAttribute("connectivity").read<uint64_t>() ==
               initial_connectivity_id + (2 + 3 * (parallel::size() > 1)));
       REQUIRE(mesh2.getAttribute("dimension").read<uint64_t>() == 3);
-      REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + (2 + 3 * (parallel::size() > 1)));
+      REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
       REQUIRE(mesh2.getAttribute("type").read<std::string>() == "polygonal");
 
       HighFive::Group functions = checkpoint.getGroup("functions");
@@ -488,22 +486,21 @@ for(let i:N, i=0; i<3; ++i) {
 
       HighFive::Group m1d = embedded1.getGroup("m1d");
       REQUIRE(m1d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (1 + 2 * (parallel::size() > 1)));
+      REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
 
       HighFive::Group m2d = embedded1.getGroup("m2d");
       REQUIRE(m2d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+      REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
       HighFive::Group m3d = embedded1.getGroup("m3d");
       REQUIRE(m3d.getAttribute("type").read<std::string>() == "mesh");
-      REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + (2 + 3 * (parallel::size() > 1)));
+      REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
 
       HighFive::Group singleton        = checkpoint.getGroup("singleton");
       HighFive::Group global_variables = singleton.getGroup("global_variables");
       REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() ==
               initial_connectivity_id + 3 * (1 + (parallel::size() > 1)));
-      REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() ==
-              initial_mesh_id + 3 * (1 + (parallel::size() > 1)));
+      REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 6);
       HighFive::Group execution_info = singleton.getGroup("execution_info");
       REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 2);
 
@@ -529,24 +526,24 @@ for(let i:N, i=0; i<3; ++i) {
       REQUIRE(connectivity2.getAttribute("type").read<std::string>() == "unstructured");
 
       HighFive::Group mesh  = checkpoint.getGroup("mesh");
-      HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + (parallel::size() > 1)));
+      HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
       REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + (parallel::size() > 1));
       REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-      REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + (parallel::size() > 1));
+      REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
       REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
-      HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + (1 + 2 * (parallel::size() > 1))));
+      HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
       REQUIRE(mesh1.getAttribute("connectivity").read<uint64_t>() ==
               initial_connectivity_id + (1 + 2 * (parallel::size() > 1)));
       REQUIRE(mesh1.getAttribute("dimension").read<uint64_t>() == 1);
-      REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + (1 + 2 * (parallel::size() > 1)));
+      REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
       REQUIRE(mesh1.getAttribute("type").read<std::string>() == "polygonal");
 
-      HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + (2 + 3 * (parallel::size() > 1))));
+      HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
       REQUIRE(mesh2.getAttribute("connectivity").read<uint64_t>() ==
               initial_connectivity_id + (2 + 3 * (parallel::size() > 1)));
       REQUIRE(mesh2.getAttribute("dimension").read<uint64_t>() == 3);
-      REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + (2 + 3 * (parallel::size() > 1)));
+      REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
       REQUIRE(mesh2.getAttribute("type").read<std::string>() == "polygonal");
 
       HighFive::Group functions = checkpoint.getGroup("functions");
diff --git a/tests/test_checkpointing_Resume_sequential.cpp b/tests/test_checkpointing_Resume_sequential.cpp
index e8a3aa1038a3089c5188bfb2dc19f4dc6a4abdbc..f2ea03bbee24e6568db15b10001ec20c4f49c25a 100644
--- a/tests/test_checkpointing_Resume_sequential.cpp
+++ b/tests/test_checkpointing_Resume_sequential.cpp
@@ -233,37 +233,37 @@ for(let i:N, i=0; i<3; ++i) {
 
     HighFive::Group m1d = embedded1.getGroup("m1d");
     REQUIRE(m1d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(m1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
 
     HighFive::Group dual_1d = embedded1.getGroup("dual_1d");
     REQUIRE(dual_1d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(dual_1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 4);
+    REQUIRE(dual_1d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
 
     HighFive::Group m2d = embedded1.getGroup("m2d");
     REQUIRE(m2d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
     HighFive::Group duals_2d = embedded1.getGroup("duals_2d");
     REQUIRE(duals_2d.getAttribute("type").read<std::string>() == "(mesh)");
     HighFive::Group duals_2d_0 = duals_2d.getGroup("0");
     REQUIRE(duals_2d_0.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_2d_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
+    REQUIRE(duals_2d_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
     HighFive::Group duals_2d_1 = duals_2d.getGroup("1");
     REQUIRE(duals_2d_1.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_2d_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(duals_2d_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
 
     HighFive::Group m3d = embedded1.getGroup("m3d");
     REQUIRE(m3d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
+    REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 8);
 
     HighFive::Group dual_3d = embedded1.getGroup("dual_3d");
     REQUIRE(dual_3d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(dual_3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
+    REQUIRE(dual_3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 9);
 
     HighFive::Group singleton        = checkpoint.getGroup("singleton");
     HighFive::Group global_variables = singleton.getGroup("global_variables");
     REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() == initial_connectivity_id + 7);
-    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 7);
+    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 10);
     HighFive::Group execution_info = singleton.getGroup("execution_info");
     REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 1);
 
@@ -308,45 +308,45 @@ for(let i:N, i=0; i<3; ++i) {
     REQUIRE(connectivity6.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
     HighFive::Group mesh  = checkpoint.getGroup("mesh");
-    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id));
+    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
     REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id);
     REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
-    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
-    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
+    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh1.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh1.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
-    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
-    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
-    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
+    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh2.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh2.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Median);
 
-    HighFive::Group mesh3 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
+    HighFive::Group mesh3 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
     REQUIRE(mesh3.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + 3);
     REQUIRE(mesh3.getAttribute("dimension").read<uint64_t>() == 1);
-    REQUIRE(mesh3.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(mesh3.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
     REQUIRE(mesh3.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh4 = mesh.getGroup(std::to_string(initial_mesh_id + 4));
-    REQUIRE(mesh4.getAttribute("id").read<uint64_t>() == initial_mesh_id + 4);
-    REQUIRE(mesh4.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 3);
+    HighFive::Group mesh4 = mesh.getGroup(std::to_string(initial_mesh_id + 6));
+    REQUIRE(mesh4.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
+    REQUIRE(mesh4.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 5);
     REQUIRE(mesh4.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh4.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Dual1D);
 
-    HighFive::Group mesh5 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
+    HighFive::Group mesh5 = mesh.getGroup(std::to_string(initial_mesh_id + 8));
     REQUIRE(mesh5.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + 5);
     REQUIRE(mesh5.getAttribute("dimension").read<uint64_t>() == 3);
-    REQUIRE(mesh5.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
+    REQUIRE(mesh5.getAttribute("id").read<uint64_t>() == initial_mesh_id + 8);
     REQUIRE(mesh5.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh6 = mesh.getGroup(std::to_string(initial_mesh_id + 6));
-    REQUIRE(mesh6.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
-    REQUIRE(mesh6.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 5);
+    HighFive::Group mesh6 = mesh.getGroup(std::to_string(initial_mesh_id + 9));
+    REQUIRE(mesh6.getAttribute("id").read<uint64_t>() == initial_mesh_id + 9);
+    REQUIRE(mesh6.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 8);
     REQUIRE(mesh6.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh6.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
@@ -417,29 +417,29 @@ for(let i:N, i=0; i<3; ++i) {
 
     HighFive::Group m2d = embedded1.getGroup("m2d");
     REQUIRE(m2d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
     HighFive::Group duals_2d = embedded1.getGroup("duals_2d");
     REQUIRE(duals_2d.getAttribute("type").read<std::string>() == "(mesh)");
     HighFive::Group duals_2d_0 = duals_2d.getGroup("0");
     REQUIRE(duals_2d_0.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_2d_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
+    REQUIRE(duals_2d_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
     HighFive::Group duals_2d_1 = duals_2d.getGroup("1");
     REQUIRE(duals_2d_1.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_2d_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(duals_2d_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
 
     HighFive::Group m3d = embedded1.getGroup("m3d");
     REQUIRE(m3d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
+    REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 8);
 
     HighFive::Group dual_3d = embedded1.getGroup("dual_3d");
     REQUIRE(dual_3d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(dual_3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
+    REQUIRE(dual_3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 9);
 
     HighFive::Group singleton        = checkpoint.getGroup("singleton");
     HighFive::Group global_variables = singleton.getGroup("global_variables");
     REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() == initial_connectivity_id + 7);
-    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 7);
+    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 10);
     HighFive::Group execution_info = singleton.getGroup("execution_info");
     REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 1);
 
@@ -484,45 +484,45 @@ for(let i:N, i=0; i<3; ++i) {
     REQUIRE(connectivity6.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
     HighFive::Group mesh  = checkpoint.getGroup("mesh");
-    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id));
+    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
     REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id);
     REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
-    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
-    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
+    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh1.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh1.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
-    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
-    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
-    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
+    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh2.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh2.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Median);
 
-    HighFive::Group mesh3 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
+    HighFive::Group mesh3 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
     REQUIRE(mesh3.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + 3);
     REQUIRE(mesh3.getAttribute("dimension").read<uint64_t>() == 1);
-    REQUIRE(mesh3.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(mesh3.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
     REQUIRE(mesh3.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh4 = mesh.getGroup(std::to_string(initial_mesh_id + 4));
-    REQUIRE(mesh4.getAttribute("id").read<uint64_t>() == initial_mesh_id + 4);
-    REQUIRE(mesh4.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 3);
+    HighFive::Group mesh4 = mesh.getGroup(std::to_string(initial_mesh_id + 6));
+    REQUIRE(mesh4.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
+    REQUIRE(mesh4.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 5);
     REQUIRE(mesh4.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh4.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Dual1D);
 
-    HighFive::Group mesh5 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
+    HighFive::Group mesh5 = mesh.getGroup(std::to_string(initial_mesh_id + 8));
     REQUIRE(mesh5.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + 5);
     REQUIRE(mesh5.getAttribute("dimension").read<uint64_t>() == 3);
-    REQUIRE(mesh5.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
+    REQUIRE(mesh5.getAttribute("id").read<uint64_t>() == initial_mesh_id + 8);
     REQUIRE(mesh5.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh6 = mesh.getGroup(std::to_string(initial_mesh_id + 6));
-    REQUIRE(mesh6.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
-    REQUIRE(mesh6.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 5);
+    HighFive::Group mesh6 = mesh.getGroup(std::to_string(initial_mesh_id + 9));
+    REQUIRE(mesh6.getAttribute("id").read<uint64_t>() == initial_mesh_id + 9);
+    REQUIRE(mesh6.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 8);
     REQUIRE(mesh6.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh6.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
@@ -596,29 +596,29 @@ for(let i:N, i=0; i<3; ++i) {
 
     HighFive::Group m2d = embedded1.getGroup("m2d");
     REQUIRE(m2d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(m2d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
 
     HighFive::Group duals_2d = embedded1.getGroup("duals_2d");
     REQUIRE(duals_2d.getAttribute("type").read<std::string>() == "(mesh)");
     HighFive::Group duals_2d_0 = duals_2d.getGroup("0");
     REQUIRE(duals_2d_0.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_2d_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
+    REQUIRE(duals_2d_0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
     HighFive::Group duals_2d_1 = duals_2d.getGroup("1");
     REQUIRE(duals_2d_1.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(duals_2d_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(duals_2d_1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
 
     HighFive::Group m3d = embedded1.getGroup("m3d");
     REQUIRE(m3d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
+    REQUIRE(m3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 8);
 
     HighFive::Group dual_3d = embedded1.getGroup("dual_3d");
     REQUIRE(dual_3d.getAttribute("type").read<std::string>() == "mesh");
-    REQUIRE(dual_3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
+    REQUIRE(dual_3d.getAttribute("id").read<uint64_t>() == initial_mesh_id + 9);
 
     HighFive::Group singleton        = checkpoint.getGroup("singleton");
     HighFive::Group global_variables = singleton.getGroup("global_variables");
     REQUIRE(global_variables.getAttribute("connectivity_id").read<uint64_t>() == initial_connectivity_id + 7);
-    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 7);
+    REQUIRE(global_variables.getAttribute("mesh_id").read<uint64_t>() == initial_mesh_id + 10);
     HighFive::Group execution_info = singleton.getGroup("execution_info");
     REQUIRE(execution_info.getAttribute("run_number").read<uint64_t>() == 2);
 
@@ -663,45 +663,45 @@ for(let i:N, i=0; i<3; ++i) {
     REQUIRE(connectivity6.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
     HighFive::Group mesh  = checkpoint.getGroup("mesh");
-    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id));
+    HighFive::Group mesh0 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
     REQUIRE(mesh0.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id);
     REQUIRE(mesh0.getAttribute("dimension").read<uint64_t>() == 2);
-    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id);
+    REQUIRE(mesh0.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh0.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 1));
-    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 1);
-    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh1 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
+    REQUIRE(mesh1.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
+    REQUIRE(mesh1.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh1.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh1.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);
 
-    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 2));
-    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 2);
-    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id);
+    HighFive::Group mesh2 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
+    REQUIRE(mesh2.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(mesh2.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 1);
     REQUIRE(mesh2.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh2.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Median);
 
-    HighFive::Group mesh3 = mesh.getGroup(std::to_string(initial_mesh_id + 3));
+    HighFive::Group mesh3 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
     REQUIRE(mesh3.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + 3);
     REQUIRE(mesh3.getAttribute("dimension").read<uint64_t>() == 1);
-    REQUIRE(mesh3.getAttribute("id").read<uint64_t>() == initial_mesh_id + 3);
+    REQUIRE(mesh3.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
     REQUIRE(mesh3.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh4 = mesh.getGroup(std::to_string(initial_mesh_id + 4));
-    REQUIRE(mesh4.getAttribute("id").read<uint64_t>() == initial_mesh_id + 4);
-    REQUIRE(mesh4.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 3);
+    HighFive::Group mesh4 = mesh.getGroup(std::to_string(initial_mesh_id + 6));
+    REQUIRE(mesh4.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
+    REQUIRE(mesh4.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 5);
     REQUIRE(mesh4.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh4.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Dual1D);
 
-    HighFive::Group mesh5 = mesh.getGroup(std::to_string(initial_mesh_id + 5));
+    HighFive::Group mesh5 = mesh.getGroup(std::to_string(initial_mesh_id + 8));
     REQUIRE(mesh5.getAttribute("connectivity").read<uint64_t>() == initial_connectivity_id + 5);
     REQUIRE(mesh5.getAttribute("dimension").read<uint64_t>() == 3);
-    REQUIRE(mesh5.getAttribute("id").read<uint64_t>() == initial_mesh_id + 5);
+    REQUIRE(mesh5.getAttribute("id").read<uint64_t>() == initial_mesh_id + 8);
     REQUIRE(mesh5.getAttribute("type").read<std::string>() == "polygonal");
 
-    HighFive::Group mesh6 = mesh.getGroup(std::to_string(initial_mesh_id + 6));
-    REQUIRE(mesh6.getAttribute("id").read<uint64_t>() == initial_mesh_id + 6);
-    REQUIRE(mesh6.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 5);
+    HighFive::Group mesh6 = mesh.getGroup(std::to_string(initial_mesh_id + 9));
+    REQUIRE(mesh6.getAttribute("id").read<uint64_t>() == initial_mesh_id + 9);
+    REQUIRE(mesh6.getAttribute("primal_mesh_id").read<uint64_t>() == initial_mesh_id + 8);
     REQUIRE(mesh6.getAttribute("type").read<std::string>() == "dual_mesh");
     REQUIRE(mesh6.getAttribute("type_of_dual").read<DualMeshType>() == DualMeshType::Diamond);