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/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.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.hpp b/src/mesh/Connectivity.hpp
index aa1a471d6e84b554c1715ca7eecc7a2361744b07..faee4c28c34a879117c6418f2fe13e0f753a1e22 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -659,7 +659,7 @@ class Connectivity final : public IConnectivity
   CRSGraph
   cellToCellGraph() const
   {
-    std::vector<std::set<int>> cell_cells(this->numberOfCells());
+    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) {
@@ -668,23 +668,36 @@ class Connectivity final : public IConnectivity
         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);
+        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());
+
+    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<int> entries(this->numberOfCells() + 1);
+    Array<int> entries(owned_cell_cells.size() + 1);
     entries[0] = 0;
-    for (size_t j = 0; j < this->numberOfCells(); ++j) {
-      entries[j + 1] = entries[j] + cell_cells[j].size();
+    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[this->numberOfCells()]);
+    Array<int> neighbors(entries[owned_cell_cells.size()]);
     {
-      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;
+      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] = m_cell_global_index[cell_id];
+          ++i;
         }
       }
     }
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index 65326b56bd49c1ce155dbf07fe32da8622c2e66c..65133fda261190c008b7169c4d32090488379237 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -7,7 +7,7 @@
 #include <iostream>
 #include <unordered_map>
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildNewOwner()
@@ -15,8 +15,18 @@ ConnectivityDispatcher<Dimension>::_buildNewOwner()
   if constexpr (item_type == ItemType::cell) {
     CRSGraph connectivity_graph = m_connectivity.cellToCellGraph();
     Partitioner P;
+    Array new_owner_array   = P.partition(connectivity_graph);
+    CellValue cell_is_owned = m_connectivity.cellIsOwned();
 
-    CellValue<int> cell_new_owner(m_connectivity, P.partition(connectivity_graph));
+    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 {
@@ -46,7 +56,7 @@ ConnectivityDispatcher<Dimension>::_buildNewOwner()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildItemToExchangeLists()
@@ -59,7 +69,7 @@ ConnectivityDispatcher<Dimension>::_buildItemToExchangeLists()
   this->_buildRecvItemIdCorrespondanceByProc<item_type>();
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildItemListToSend()
@@ -126,7 +136,7 @@ ConnectivityDispatcher<Dimension>::_buildItemListToSend()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <ItemType item_type>
 void
 ConnectivityDispatcher<Dimension>::_buildNumberOfItemToExchange()
@@ -141,7 +151,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,
@@ -162,7 +172,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,7 +231,7 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 void
 ConnectivityDispatcher<Dimension>::_buildCellNumberIdMap()
 {
@@ -238,7 +248,7 @@ ConnectivityDispatcher<Dimension>::_buildCellNumberIdMap()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename ItemOfItemT>
 void
 ConnectivityDispatcher<Dimension>::_buildSubItemNumberToIdMap()
@@ -261,7 +271,7 @@ ConnectivityDispatcher<Dimension>::_buildSubItemNumberToIdMap()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename SubItemOfItemT>
 void
 ConnectivityDispatcher<Dimension>::_buildNumberOfSubItemPerItemToRecvByProc()
@@ -280,7 +290,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 +340,7 @@ ConnectivityDispatcher<Dimension>::_buildSubItemNumbersToRecvByProc()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 template <typename ItemOfItemT>
 void
 ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
@@ -385,7 +395,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()
@@ -426,7 +436,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()
@@ -601,7 +611,7 @@ ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 void
 ConnectivityDispatcher<Dimension>::_dispatchEdges()
 {
@@ -643,7 +653,7 @@ ConnectivityDispatcher<Dimension>::_dispatchEdges()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 void
 ConnectivityDispatcher<Dimension>::_dispatchFaces()
 {
@@ -681,7 +691,7 @@ ConnectivityDispatcher<Dimension>::_dispatchFaces()
   }
 }
 
-template <int Dimension>
+template <size_t Dimension>
 ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType& connectivity)
   : m_connectivity(connectivity)
 {
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index 9ff5913b60799148f99bb6fdf4cff2cedd56973c..9d9d6a14fa95f10a8b78553c0cd69732b8804393 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:
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/MeshBalancer.cpp b/src/mesh/MeshBalancer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d076c7d48b47d188c83e30b37bb53a877e4be555
--- /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>();
+    },
+    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..861050ed7e3cfd070ab13d6eaa7eed489f4352a4 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>
@@ -33,10 +34,13 @@ MeshBuilderBase::_dispatch()
   }
 
   const MeshType& mesh = *(m_mesh->get<const MeshType>());
-  ConnectivityDispatcher<Dimension> dispatcher(mesh.connectivity());
 
-  std::shared_ptr dispatched_connectivity = dispatcher.dispatchedConnectivity();
-  NodeValue<Rd> dispatched_xr             = dispatcher.dispatch(mesh.xr());
+  auto p_dispatcher = std::make_shared<const ConnectivityDispatcher<Dimension>>(mesh.connectivity());
+
+  m_connectivity_dispatcher = std::make_shared<ConnectivityDispatcherVariant>(p_dispatcher);
+
+  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));
 }
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
index 665f0a88cc5a7f766ceb792594c6f5ace23529d6..fb811acf499a7b4bedd3f124495fa7bd24eb92a7 100644
--- a/src/mesh/MeshBuilderBase.hpp
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -2,6 +2,7 @@
 #define MESH_BUILDER_BASE_HPP
 
 class MeshVariant;
+class ConnectivityDispatcherVariant;
 
 #include <memory>
 
@@ -9,6 +10,7 @@ class MeshBuilderBase
 {
  protected:
   std::shared_ptr<const MeshVariant> m_mesh;
+  std::shared_ptr<const ConnectivityDispatcherVariant> m_connectivity_dispatcher;
 
   template <size_t Dimension>
   void _dispatch();
@@ -23,8 +25,8 @@ class MeshBuilderBase
     return m_mesh;
   }
 
-  MeshBuilderBase()  = default;
-  ~MeshBuilderBase() = default;
+  MeshBuilderBase()          = default;
+  virtual ~MeshBuilderBase() = default;
 };
 
 #endif   // MESH_BUILDER_BASE_HPP
diff --git a/src/utils/ParMETISPartitioner.cpp b/src/utils/ParMETISPartitioner.cpp
index 47864de1fd3337bc6e193ad3453ad8d6c75eb797..7c7d81ca712fe285ec967555493fc153faada5b8 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();
@@ -62,6 +77,12 @@ ParMETISPartitioner::partition(const CRSGraph& graph)
     int* entries_ptr   = const_cast<int*>(&(entries[0]));
     int* neighbors_ptr = const_cast<int*>(&(neighbors[0]));
 
+    if (group_size > 1) {
+#warning remove
+      std::cout << parallel::rank() << ": entries " << graph.entries() << "\n";
+      std::cout << parallel::rank() << ": neighbors " << graph.neighbors() << "\n";
+    }
+
     int result =
       ParMETIS_V3_PartKway(&(vtxdist[0]), entries_ptr, neighbors_ptr, NULL, NULL, &wgtflag, &numflag, &ncon, &npart,
                            &(tpwgts[0]), &(ubvec[0]), &(options[0]), &edgecut, &(partition[0]), &partitioning_comm);
@@ -69,6 +90,12 @@ ParMETISPartitioner::partition(const CRSGraph& graph)
     if (result == METIS_ERROR) {
       throw UnexpectedError("Metis Error");
     }
+
+    if (group_size > 1) {
+#warning remove
+      std::cout << parallel::rank() << ": partition " << partition << "\n";
+    }
+
     // LCOV_EXCL_STOP
   }