diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index defdc8fc9711c6be265e022e0bb0c710d078226d..f48907baac28bb49d0921d11741f463dce791642 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -49,7 +49,7 @@ class SubItemValuePerItem
 
  public:
   friend PUGS_INLINE SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
-  copy(SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
+  copy(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
     SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr> image;
 
diff --git a/src/mesh/SubItemValuePerItemUtils.hpp b/src/mesh/SubItemValuePerItemUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5ff1538764e495f85f654d70f030964f900a43d1
--- /dev/null
+++ b/src/mesh/SubItemValuePerItemUtils.hpp
@@ -0,0 +1,350 @@
+#ifndef SUB_ITEM_VALUE_PER_ITEM_UTILS_HPP
+#define SUB_ITEM_VALUE_PER_ITEM_UTILS_HPP
+
+#include <utils/Messenger.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+#include <mesh/Synchronizer.hpp>
+#include <mesh/SynchronizerManager.hpp>
+#include <utils/PugsTraits.hpp>
+
+#include <iostream>
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+min(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_value_per_item)
+{
+  using SubItemValuePerItemType = SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>;
+  constexpr ItemType item_type  = ItemOfItem::item_type;
+  using ItemIsOwnedType         = ItemValue<const bool, item_type>;
+  using data_type               = std::remove_const_t<typename SubItemValuePerItemType::data_type>;
+  using index_type              = typename SubItemValuePerItemType::index_type;
+
+  static_assert(std::is_arithmetic_v<data_type>, "min cannot be called on non-arithmetic data");
+  static_assert(not std::is_same_v<data_type, bool>, "min cannot be called on boolean data");
+
+  class SubItemValuePerItemMin
+  {
+   private:
+    const SubItemValuePerItemType& m_sub_item_value_per_item;
+    const ItemIsOwnedType m_is_owned;
+
+   public:
+    PUGS_INLINE
+    operator data_type()
+    {
+      data_type reduced_value;
+      parallel_reduce(m_sub_item_value_per_item.numberOfItems(), *this, reduced_value);
+      return reduced_value;
+    }
+
+    PUGS_INLINE
+    void
+    operator()(const index_type& item_id, data_type& data) const
+    {
+      if (m_is_owned[item_id]) {
+        Array sub_item_array = m_sub_item_value_per_item.itemValues(item_id);
+        for (size_t i = 0; i < sub_item_array.size(); ++i) {
+          if (sub_item_array[i] < data) {
+            data = sub_item_array[i];
+          }
+        }
+      }
+    }
+
+    PUGS_INLINE
+    void
+    join(volatile data_type& dst, const volatile data_type& src) const
+    {
+      if (src < dst) {
+        // cannot be reached if initial value is the min
+        dst = src;   // LCOV_EXCL_LINE
+      }
+    }
+
+    PUGS_INLINE
+    void
+    init(data_type& value) const
+    {
+      value = std::numeric_limits<data_type>::max();
+    }
+
+    PUGS_INLINE
+    SubItemValuePerItemMin(const SubItemValuePerItemType& item_value)
+      : m_sub_item_value_per_item(item_value), m_is_owned([&](const IConnectivity& connectivity) {
+          Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+                 "unexpected connectivity dimension");
+
+          switch (connectivity.dimension()) {
+          case 1: {
+            const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+            return connectivity_1d.isOwned<item_type>();
+            break;
+          }
+          case 2: {
+            const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+            return connectivity_2d.isOwned<item_type>();
+            break;
+          }
+          case 3: {
+            const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+            return connectivity_3d.isOwned<item_type>();
+            break;
+          }
+            // LCOV_EXCL_START
+          default: {
+            throw UnexpectedError("unexpected dimension");
+          }
+            // LCOV_EXCL_STOP
+          }
+        }(*item_value.connectivity_ptr()))
+    {
+      ;
+    }
+
+    PUGS_INLINE
+    ~SubItemValuePerItemMin() = default;
+  };
+
+  const DataType local_min = SubItemValuePerItemMin{sub_item_value_per_item};
+  return parallel::allReduceMin(local_min);
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+max(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_value_per_item)
+{
+  using SubItemValuePerItemType = SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>;
+  constexpr ItemType item_type  = ItemOfItem::item_type;
+  using ItemIsOwnedType         = ItemValue<const bool, item_type>;
+  using data_type               = std::remove_const_t<typename SubItemValuePerItemType::data_type>;
+  using index_type              = typename SubItemValuePerItemType::index_type;
+
+  static_assert(std::is_arithmetic_v<data_type>, "min cannot be called on non-arithmetic data");
+  static_assert(not std::is_same_v<data_type, bool>, "max cannot be called on boolean data");
+
+  class SubItemValuePerItemMax
+  {
+   private:
+    const SubItemValuePerItemType& m_sub_item_value_per_item;
+    const ItemIsOwnedType m_is_owned;
+
+   public:
+    PUGS_INLINE
+    operator data_type()
+    {
+      data_type reduced_value;
+      parallel_reduce(m_sub_item_value_per_item.numberOfItems(), *this, reduced_value);
+      return reduced_value;
+    }
+
+    PUGS_INLINE
+    void
+    operator()(const index_type& item_id, data_type& data) const
+    {
+      if (m_is_owned[item_id]) {
+        Array sub_item_array = m_sub_item_value_per_item.itemValues(item_id);
+        for (size_t i = 0; i < sub_item_array.size(); ++i) {
+          if (sub_item_array[i] > data) {
+            data = sub_item_array[i];
+          }
+        }
+      }
+    }
+
+    PUGS_INLINE
+    void
+    join(volatile data_type& dst, const volatile data_type& src) const
+    {
+      if (src > dst) {
+        // cannot be reached if initial value is the max
+        dst = src;   // LCOV_EXCL_LINE
+      }
+    }
+
+    PUGS_INLINE
+    void
+    init(data_type& value) const
+    {
+      value = std::numeric_limits<data_type>::min();
+    }
+
+    PUGS_INLINE
+    SubItemValuePerItemMax(const SubItemValuePerItemType& item_value)
+      : m_sub_item_value_per_item(item_value), m_is_owned([&](const IConnectivity& connectivity) {
+          Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+                 "unexpected connectivity dimension");
+
+          switch (connectivity.dimension()) {
+          case 1: {
+            const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+            return connectivity_1d.isOwned<item_type>();
+            break;
+          }
+          case 2: {
+            const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+            return connectivity_2d.isOwned<item_type>();
+            break;
+          }
+          case 3: {
+            const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+            return connectivity_3d.isOwned<item_type>();
+            break;
+          }
+            // LCOV_EXCL_START
+          default: {
+            throw UnexpectedError("unexpected dimension");
+          }
+            // LCOV_EXCL_STOP
+          }
+        }(*item_value.connectivity_ptr()))
+    {
+      ;
+    }
+
+    PUGS_INLINE
+    ~SubItemValuePerItemMax() = default;
+  };
+
+  const DataType local_max = SubItemValuePerItemMax{sub_item_value_per_item};
+  return parallel::allReduceMax(local_max);
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+sum(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& item_value)
+{
+  using SubItemValuePerItemType = SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>;
+  constexpr ItemType item_type  = ItemOfItem::item_type;
+  using ItemIsOwnedType         = ItemValue<const bool, item_type>;
+  using data_type               = std::remove_const_t<typename SubItemValuePerItemType::data_type>;
+  using index_type              = typename SubItemValuePerItemType::index_type;
+
+  static_assert(not std::is_same_v<data_type, bool>, "sum cannot be called on boolean data");
+
+  class SubItemValuePerItemSum
+  {
+   private:
+    const SubItemValuePerItemType& m_sub_item_value_per_item;
+    const ItemIsOwnedType m_is_owned;
+
+   public:
+    PUGS_INLINE
+    operator data_type()
+    {
+      data_type reduced_value;
+      parallel_reduce(m_sub_item_value_per_item.numberOfItems(), *this, reduced_value);
+      return reduced_value;
+    }
+
+    PUGS_INLINE
+    void
+    operator()(const index_type& item_id, data_type& data) const
+    {
+      if (m_is_owned[item_id]) {
+        Array sub_item_array = m_sub_item_value_per_item.itemValues(item_id);
+        for (size_t i = 0; i < sub_item_array.size(); ++i) {
+          data += sub_item_array[i];
+        }
+      }
+    }
+
+    PUGS_INLINE
+    void
+    join(volatile data_type& dst, const volatile data_type& src) const
+    {
+      dst += src;
+    }
+
+    PUGS_INLINE
+    void
+    init(data_type& value) const
+    {
+      if constexpr (std::is_arithmetic_v<data_type>) {
+        value = 0;
+      } else {
+        static_assert(is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>, "invalid data type");
+        value = zero;
+      }
+    }
+
+    PUGS_INLINE
+    SubItemValuePerItemSum(const SubItemValuePerItemType& item_value)
+      : m_sub_item_value_per_item(item_value), m_is_owned([&](const IConnectivity& connectivity) {
+          Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+                 "unexpected connectivity dimension");
+
+          switch (connectivity.dimension()) {
+          case 1: {
+            const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+            return connectivity_1d.isOwned<item_type>();
+            break;
+          }
+          case 2: {
+            const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+            return connectivity_2d.isOwned<item_type>();
+            break;
+          }
+          case 3: {
+            const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+            return connectivity_3d.isOwned<item_type>();
+            break;
+          }
+            // LCOV_EXCL_START
+          default: {
+            throw UnexpectedError("unexpected dimension");
+          }
+            // LCOV_EXCL_STOP
+          }
+        }(*item_value.connectivity_ptr()))
+    {
+      ;
+    }
+
+    PUGS_INLINE
+    ~SubItemValuePerItemSum() = default;
+  };
+
+  const DataType local_sum = SubItemValuePerItemSum{item_value};
+  return parallel::allReduceSum(local_sum);
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+void
+synchronize(SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_value_per_item)
+{
+  static_assert(not std::is_const_v<DataType>, "cannot synchronize SubItemValuePerItem of const data");
+  if (parallel::size() > 1) {
+    auto& manager                     = SynchronizerManager::instance();
+    const IConnectivity* connectivity = sub_item_value_per_item.connectivity_ptr().get();
+    Synchronizer& synchronizer        = manager.getConnectivitySynchronizer(connectivity);
+    synchronizer.synchronize(sub_item_value_per_item);
+  }
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+bool
+isSynchronized(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_value_per_item)
+{
+  bool is_synchronized = true;
+  if (parallel::size() > 1) {
+    SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem> sub_item_value_per_item_copy =
+      copy(sub_item_value_per_item);
+
+    synchronize(sub_item_value_per_item_copy);
+
+    for (size_t i = 0; i < sub_item_value_per_item.numberOfValues(); ++i) {
+      if (sub_item_value_per_item_copy[i] != sub_item_value_per_item[i]) {
+        is_synchronized = false;
+        break;
+      }
+    }
+
+    is_synchronized = parallel::allReduceAnd(is_synchronized);
+  }
+
+  return is_synchronized;
+}
+
+#endif   // SUB_ITEM_VALUE_PER_ITEM_UTILS_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ddfa977a6ef8ca76df9f27bbdc2bda86395cc67f..5943849313535795bef5754e678ae2e0f383fc1c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -173,6 +173,7 @@ add_executable (mpi_unit_tests
   test_Partitioner.cpp
   test_RandomEngine.cpp
   test_SubItemValuePerItem.cpp
+  test_SubItemValuePerItemUtils.cpp
   test_SubItemArrayPerItem.cpp
   test_Synchronizer.cpp
   )
diff --git a/tests/test_SubItemValuePerItemUtils.cpp b/tests/test_SubItemValuePerItemUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c93f5cb7c2c48ded7693b14c7f4066baefe18906
--- /dev/null
+++ b/tests/test_SubItemValuePerItemUtils.cpp
@@ -0,0 +1,409 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+#include <mesh/SubItemValuePerItemUtils.hpp>
+#include <utils/Messenger.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SubItemValuePerItemUtils", "[mesh]")
+{
+  SECTION("Synchronize")
+  {
+    std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_2d = named_mesh.mesh();
+
+        const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+        NodeValuePerFace<int> weak_node_value_per_face{connectivity};
+
+        for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+          auto face_array = weak_node_value_per_face.itemValues(face_id);
+          for (size_t i_node = 0; i_node < face_array.size(); ++i_node) {
+            face_array[i_node] = parallel::rank() + 2 * i_node;
+          }
+        }
+
+        NodeValuePerFace<const int> node_value_per_face{weak_node_value_per_face};
+
+        REQUIRE(node_value_per_face.connectivity_ptr() == weak_node_value_per_face.connectivity_ptr());
+
+        if (parallel::size() > 1) {
+          // before synchronization
+          auto face_owner = connectivity.faceOwner();
+
+          bool is_synchronized = true;
+          for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+            auto face_array = node_value_per_face.itemValues(face_id);
+            for (size_t i_node = 0; i_node < face_array.size(); ++i_node) {
+              if (face_array[i_node] != static_cast<int>(face_owner[face_id] + 2 * i_node)) {
+                is_synchronized = false;
+                break;
+              }
+            }
+          }
+
+          REQUIRE(not is_synchronized);
+          REQUIRE(not isSynchronized(node_value_per_face));
+        }
+
+        synchronize(weak_node_value_per_face);
+
+        {   // after synchronization
+          auto face_owner = connectivity.faceOwner();
+
+          bool is_synchronized = true;
+          for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+            auto face_array = node_value_per_face.itemValues(face_id);
+            for (size_t i_node = 0; i_node < face_array.size(); ++i_node) {
+              if (face_array[i_node] != static_cast<int>(face_owner[face_id] + 2 * i_node)) {
+                is_synchronized = false;
+                break;
+              }
+            }
+          }
+          REQUIRE(is_synchronized);
+        }
+        REQUIRE(isSynchronized(node_value_per_face));
+      }
+    }
+  }
+
+  SECTION("min")
+  {
+    SECTION("1D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+
+          const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+          NodeValuePerCell<int> node_value_per_cell{connectivity};
+          node_value_per_cell.fill(-1);
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                auto cell_array = node_value_per_cell.itemValues(cell_id);
+                for (size_t i = 0; i < cell_array.size(); ++i) {
+                  cell_array[i] = 10 + parallel::rank() + i;
+                }
+              }
+            });
+
+          REQUIRE(min(node_value_per_cell) == 10);
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+          FaceValuePerCell<int> face_value_per_cell{connectivity};
+          face_value_per_cell.fill(-1);
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                auto cell_array = face_value_per_cell.itemValues(cell_id);
+                for (size_t i = 0; i < cell_array.size(); ++i) {
+                  cell_array[i] = 10 + parallel::rank() + i;
+                }
+              }
+            });
+
+          REQUIRE(min(face_value_per_cell) == 10);
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          EdgeValuePerFace<int> edge_value_per_face{connectivity};
+          edge_value_per_face.fill(-1);
+
+          auto face_is_owned = connectivity.faceIsOwned();
+          parallel_for(
+            connectivity.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+              if (face_is_owned[face_id]) {
+                auto face_array = edge_value_per_face.itemValues(face_id);
+                for (size_t i = 0; i < face_array.size(); ++i) {
+                  face_array[i] = 10 + parallel::rank() + i;
+                }
+              }
+            });
+
+          REQUIRE(min(edge_value_per_face) == 10);
+        }
+      }
+    }
+  }
+
+  SECTION("max")
+  {
+    SECTION("1D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+
+          const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+          EdgeValuePerCell<size_t> edge_value_per_cell{connectivity};
+          edge_value_per_cell.fill(std::numeric_limits<size_t>::max());
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                auto cell_array = edge_value_per_cell.itemValues(cell_id);
+                for (size_t i = 0; i < cell_array.size(); ++i) {
+                  cell_array[i] = 10 + parallel::rank() - i;
+                }
+              }
+            });
+
+          REQUIRE(max(edge_value_per_cell) == 9 + parallel::size());
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+          EdgeValuePerCell<size_t> edge_value_per_cell{connectivity};
+          edge_value_per_cell.fill(std::numeric_limits<size_t>::max());
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+          parallel_for(
+            connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              if (cell_is_owned[cell_id]) {
+                auto cell_array = edge_value_per_cell.itemValues(cell_id);
+                for (size_t i = 0; i < cell_array.size(); ++i) {
+                  cell_array[i] = 10 + parallel::rank() - i;
+                }
+              }
+            });
+
+          REQUIRE(max(edge_value_per_cell) == 9 + parallel::size());
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          CellValue<size_t> cell_value{connectivity};
+
+          NodeValuePerEdge<size_t> node_value_per_edge{connectivity};
+          node_value_per_edge.fill(std::numeric_limits<size_t>::max());
+
+          auto edge_is_owned = connectivity.edgeIsOwned();
+          parallel_for(
+            connectivity.numberOfEdges(), PUGS_LAMBDA(EdgeId edge_id) {
+              if (edge_is_owned[edge_id]) {
+                auto edge_array = node_value_per_edge.itemValues(edge_id);
+                for (size_t i = 0; i < edge_array.size(); ++i) {
+                  edge_array[i] = 10 + parallel::rank() - i;
+                }
+              }
+            });
+
+          REQUIRE(max(node_value_per_edge) == 9 + parallel::size());
+        }
+      }
+    }
+  }
+
+  SECTION("sum")
+  {
+    SECTION("1D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+
+          const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+          NodeValuePerCell<size_t> node_value_per_cell{connectivity};
+          node_value_per_cell.fill(5);
+
+          auto cell_is_owned = connectivity.cellIsOwned();
+
+          const size_t global_number_of_nodes_per_cells = [&] {
+            size_t number_of_nodes_per_cells = 0;
+            for (CellId cell_id = 0; cell_id < cell_is_owned.numberOfItems(); ++cell_id) {
+              number_of_nodes_per_cells += cell_is_owned[cell_id] * node_value_per_cell.numberOfSubValues(cell_id);
+            }
+            return parallel::allReduceSum(number_of_nodes_per_cells);
+          }();
+
+          REQUIRE(sum(node_value_per_cell) == 5 * global_number_of_nodes_per_cells);
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name() + " for size_t data")
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+          NodeValuePerFace<size_t> node_value_per_face{connectivity};
+          node_value_per_face.fill(2);
+
+          auto face_is_owned = connectivity.faceIsOwned();
+
+          const size_t global_number_of_nodes_per_faces = [&] {
+            size_t number_of_nodes_per_faces = 0;
+            for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+              number_of_nodes_per_faces += face_is_owned[face_id] * node_value_per_face.numberOfSubValues(face_id);
+            }
+            return parallel::allReduceSum(number_of_nodes_per_faces);
+          }();
+
+          REQUIRE(sum(node_value_per_face) == 2 * global_number_of_nodes_per_faces);
+        }
+
+        SECTION(named_mesh.name() + " for N^2 data")
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+          using N2 = TinyVector<2, size_t>;
+          NodeValuePerFace<N2> node_value_per_face{connectivity};
+
+          const N2 data(3, 2);
+          node_value_per_face.fill(data);
+
+          auto face_is_owned = connectivity.faceIsOwned();
+
+          const size_t global_number_of_nodes_per_faces = [&] {
+            size_t number_of_nodes_per_faces = 0;
+            for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+              number_of_nodes_per_faces += face_is_owned[face_id] * node_value_per_face.numberOfSubValues(face_id);
+            }
+            return parallel::allReduceSum(number_of_nodes_per_faces);
+          }();
+
+          REQUIRE(sum(node_value_per_face) == global_number_of_nodes_per_faces * data);
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name() + " for size_t data")
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          EdgeValuePerFace<size_t> edge_value_per_face{connectivity};
+          edge_value_per_face.fill(3);
+
+          auto face_is_owned = connectivity.faceIsOwned();
+
+          const size_t global_number_of_edges_per_faces = [&] {
+            size_t number_of_edges_per_faces = 0;
+            for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+              number_of_edges_per_faces += face_is_owned[face_id] * edge_value_per_face.numberOfSubValues(face_id);
+            }
+            return parallel::allReduceSum(number_of_edges_per_faces);
+          }();
+
+          REQUIRE(sum(edge_value_per_face) == 3 * global_number_of_edges_per_faces);
+        }
+
+        SECTION(named_mesh.name() + " for N^2x3 data")
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          using N2x3 = TinyMatrix<2, 3, size_t>;
+          EdgeValuePerFace<N2x3> edge_value_per_face{connectivity};
+
+          const N2x3 data(1, 3, 4, 6, 2, 5);
+          edge_value_per_face.fill(data);
+
+          auto face_is_owned = connectivity.faceIsOwned();
+
+          const size_t global_number_of_edges_per_faces = [&] {
+            size_t number_of_edges_per_faces = 0;
+            for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+              number_of_edges_per_faces += face_is_owned[face_id] * edge_value_per_face.numberOfSubValues(face_id);
+            }
+            return parallel::allReduceSum(number_of_edges_per_faces);
+          }();
+
+          REQUIRE(sum(edge_value_per_face) == global_number_of_edges_per_faces * data);
+        }
+      }
+    }
+  }
+}