diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 61c71097381797b42d7526abc1be48003ef765ca..b9503526030034722a782aea8922eb5711f754b3 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -102,6 +102,7 @@ add_executable (unit_tests
 add_executable (mpi_unit_tests
   mpi_test_main.cpp
   test_DiscreteFunctionP0.cpp
+  test_DiscreteFunctionP0Vector.cpp
   test_ItemArray.cpp
   test_ItemArrayUtils.cpp
   test_ItemValue.cpp
diff --git a/tests/test_DiscreteFunctionP0Vector.cpp b/tests/test_DiscreteFunctionP0Vector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c004ba5732c41dce1af56d3074b9c038b3a1bf8c
--- /dev/null
+++ b/tests/test_DiscreteFunctionP0Vector.cpp
@@ -0,0 +1,845 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionP0Vector", "[scheme]")
+{
+  auto same_values = [](const auto& f, const auto& g) {
+    const size_t number_of_cells = f.cellArrays().numberOfItems();
+    const size_t size_of_arrays  = f.cellArrays().sizeOfArrays();
+    for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+      for (size_t i = 0; i < size_of_arrays; ++i) {
+        if (f[cell_id][i] != g[cell_id][i]) {
+          return false;
+        }
+      }
+    }
+    return true;
+  };
+
+  SECTION("constructors")
+  {
+    SECTION("1D")
+    {
+      const size_t size = 3;
+
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh1D();
+      constexpr size_t Dimension = 1;
+
+      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0Vector);
+      REQUIRE(f.size() == size);
+
+      REQUIRE(f.mesh().get() == mesh.get());
+
+      DiscreteFunctionP0Vector g{f};
+      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0Vector);
+      REQUIRE(g.size() == size);
+
+      CellArray<double> h_arrays{mesh->connectivity(), size};
+      h_arrays.fill(0);
+
+      DiscreteFunctionP0Vector zero{mesh, [&] {
+                                      CellArray<double> cell_array{mesh->connectivity(), size};
+                                      cell_array.fill(0);
+                                      return cell_array;
+                                    }()};
+
+      DiscreteFunctionP0Vector h{mesh, h_arrays};
+      REQUIRE(same_values(h, zero));
+      REQUIRE(same_values(h, h_arrays));
+
+      h_arrays.fill(1);
+
+      REQUIRE(same_values(h, h_arrays));
+      REQUIRE(not same_values(h, zero));
+
+      DiscreteFunctionP0Vector moved_h{std::move(h)};
+      REQUIRE(same_values(moved_h, h_arrays));
+    }
+
+    SECTION("2D")
+    {
+      const size_t size = 3;
+
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh2D();
+      constexpr size_t Dimension = 2;
+
+      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0Vector);
+      REQUIRE(f.size() == size);
+
+      REQUIRE(f.mesh().get() == mesh.get());
+
+      DiscreteFunctionP0Vector g{f};
+      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0Vector);
+      REQUIRE(g.size() == size);
+
+      CellArray<double> h_arrays{mesh->connectivity(), size};
+      h_arrays.fill(0);
+
+      DiscreteFunctionP0Vector zero{mesh, [&] {
+                                      CellArray<double> cell_array{mesh->connectivity(), size};
+                                      cell_array.fill(0);
+                                      return cell_array;
+                                    }()};
+
+      DiscreteFunctionP0Vector h{mesh, h_arrays};
+      REQUIRE(same_values(h, zero));
+      REQUIRE(same_values(h, h_arrays));
+
+      h_arrays.fill(1);
+
+      REQUIRE(same_values(h, h_arrays));
+      REQUIRE(not same_values(h, zero));
+
+      DiscreteFunctionP0Vector moved_h{std::move(h)};
+      REQUIRE(same_values(moved_h, h_arrays));
+    }
+
+    SECTION("3D")
+    {
+      const size_t size = 2;
+
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh3D();
+      constexpr size_t Dimension = 3;
+
+      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0Vector);
+      REQUIRE(f.size() == size);
+
+      REQUIRE(f.mesh().get() == mesh.get());
+
+      DiscreteFunctionP0Vector g{f};
+      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0Vector);
+      REQUIRE(g.size() == size);
+
+      CellArray<double> h_arrays{mesh->connectivity(), size};
+      h_arrays.fill(0);
+
+      DiscreteFunctionP0Vector zero{mesh, [&] {
+                                      CellArray<double> cell_array{mesh->connectivity(), size};
+                                      cell_array.fill(0);
+                                      return cell_array;
+                                    }()};
+
+      DiscreteFunctionP0Vector h{mesh, h_arrays};
+      REQUIRE(same_values(h, zero));
+      REQUIRE(same_values(h, h_arrays));
+
+      h_arrays.fill(1);
+
+      REQUIRE(same_values(h, h_arrays));
+      REQUIRE(not same_values(h, zero));
+
+      DiscreteFunctionP0Vector moved_h{std::move(h)};
+      REQUIRE(same_values(moved_h, h_arrays));
+    }
+  }
+
+  SECTION("fill")
+  {
+    auto all_values_equal = [](const auto& f, const auto& g) {
+      const size_t number_of_cells = f.cellArrays().numberOfItems();
+      size_t size_of_arrays        = f.cellArrays().sizeOfArrays();
+      for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+        for (size_t i = 0; i < size_of_arrays; ++i) {
+          if (f[cell_id][i] != g) {
+            return false;
+          }
+        }
+      }
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      const size_t size = 3;
+
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh1D();
+      constexpr size_t Dimension = 1;
+
+      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+      f.fill(3);
+
+      REQUIRE(all_values_equal(f, 3));
+    }
+
+    SECTION("2D")
+    {
+      const size_t size = 3;
+
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh2D();
+      constexpr size_t Dimension = 2;
+
+      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+      f.fill(2.3);
+
+      REQUIRE(all_values_equal(f, 2.3));
+    }
+
+    SECTION("3D")
+    {
+      const size_t size = 2;
+
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh3D();
+      constexpr size_t Dimension = 3;
+
+      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+      f.fill(3.2);
+
+      REQUIRE(all_values_equal(f, 3.2));
+    }
+  }
+
+  SECTION("copies")
+  {
+    auto all_values_equal = [](const auto& f, const auto& g) {
+      const size_t number_of_cells = f.cellArrays().numberOfItems();
+      const size_t size_of_arrays  = f.cellArrays().sizeOfArrays();
+
+      for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+        for (size_t i = 0; i < size_of_arrays; ++i) {
+          if (f[cell_id][i] != g) {
+            return false;
+          }
+        }
+      }
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();
+
+      constexpr size_t Dimension = 1;
+
+      const size_t size  = 3;
+      const size_t value = parallel::rank() + 1;
+      const size_t zero  = 0;
+
+      DiscreteFunctionP0Vector<Dimension, size_t> f{mesh, size};
+      f.fill(value);
+
+      REQUIRE(all_values_equal(f, value));
+
+      DiscreteFunctionP0Vector g = copy(f);
+      f.fill(zero);
+
+      REQUIRE(all_values_equal(f, zero));
+      REQUIRE(all_values_equal(g, value));
+
+      copy_to(g, f);
+      g.fill(zero);
+
+      DiscreteFunctionP0Vector<Dimension, const size_t> h = copy(f);
+
+      DiscreteFunctionP0Vector<Dimension, size_t> shallow_g{mesh, size};
+      shallow_g = g;
+
+      REQUIRE(all_values_equal(f, value));
+      REQUIRE(all_values_equal(g, zero));
+      REQUIRE(all_values_equal(shallow_g, zero));
+      REQUIRE(all_values_equal(h, value));
+
+      copy_to(h, g);
+
+      REQUIRE(all_values_equal(g, value));
+      REQUIRE(all_values_equal(shallow_g, value));
+    }
+
+    SECTION("2D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();
+
+      constexpr size_t Dimension = 2;
+
+      const size_t size  = 3;
+      const size_t value = parallel::rank() + 1;
+      const size_t zero  = 0;
+
+      DiscreteFunctionP0Vector<Dimension, size_t> f{mesh, size};
+      f.fill(value);
+
+      REQUIRE(all_values_equal(f, value));
+
+      DiscreteFunctionP0Vector g = copy(f);
+      f.fill(zero);
+
+      REQUIRE(all_values_equal(f, zero));
+      REQUIRE(all_values_equal(g, value));
+
+      copy_to(g, f);
+      g.fill(zero);
+
+      DiscreteFunctionP0Vector<Dimension, const size_t> h = copy(f);
+
+      DiscreteFunctionP0Vector<Dimension, size_t> shallow_g{mesh, size};
+      shallow_g = g;
+
+      REQUIRE(all_values_equal(f, value));
+      REQUIRE(all_values_equal(g, zero));
+      REQUIRE(all_values_equal(shallow_g, zero));
+      REQUIRE(all_values_equal(h, value));
+
+      copy_to(h, g);
+
+      REQUIRE(all_values_equal(g, value));
+      REQUIRE(all_values_equal(shallow_g, value));
+    }
+
+    SECTION("3D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();
+
+      constexpr size_t Dimension = 3;
+
+      const size_t size  = 3;
+      const size_t value = parallel::rank() + 1;
+      const size_t zero  = 0;
+
+      DiscreteFunctionP0Vector<Dimension, size_t> f{mesh, size};
+      f.fill(value);
+
+      REQUIRE(all_values_equal(f, value));
+
+      DiscreteFunctionP0Vector g = copy(f);
+      f.fill(zero);
+
+      REQUIRE(all_values_equal(f, zero));
+      REQUIRE(all_values_equal(g, value));
+
+      copy_to(g, f);
+      g.fill(zero);
+
+      DiscreteFunctionP0Vector<Dimension, const size_t> h = copy(f);
+
+      DiscreteFunctionP0Vector<Dimension, size_t> shallow_g{mesh, size};
+      shallow_g = g;
+
+      REQUIRE(all_values_equal(f, value));
+      REQUIRE(all_values_equal(g, zero));
+      REQUIRE(all_values_equal(h, value));
+
+      copy_to(h, g);
+
+      REQUIRE(all_values_equal(g, value));
+      REQUIRE(all_values_equal(shallow_g, value));
+    }
+  }
+
+  SECTION("unary operators")
+  {
+    SECTION("1D")
+    {
+      const size_t size    = 3;
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();
+
+      constexpr size_t Dimension = 1;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      SECTION("unary minus")
+      {
+        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            const double x = xj[cell_id][0];
+            for (size_t i = 0; i < size; ++i) {
+              f[cell_id][i] = 2 * x + i;
+            }
+          });
+
+        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+
+        Table<double> minus_values{mesh->numberOfCells(), size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            for (size_t i = 0; i < size; ++i) {
+              minus_values[cell_id][i] = -f[cell_id][i];
+            }
+          });
+
+        REQUIRE(same_values(-f, minus_values));
+        REQUIRE(same_values(-const_f, minus_values));
+      }
+    }
+
+    SECTION("2D")
+    {
+      const size_t size    = 3;
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();
+
+      constexpr size_t Dimension = 2;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      SECTION("unary minus")
+      {
+        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            const double x = xj[cell_id][0];
+            const double y = xj[cell_id][1];
+            for (size_t i = 0; i < size; ++i) {
+              f[cell_id][i] = 2 * x + i * y;
+            }
+          });
+
+        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+
+        Table<double> minus_values{mesh->numberOfCells(), size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            for (size_t i = 0; i < size; ++i) {
+              minus_values[cell_id][i] = -f[cell_id][i];
+            }
+          });
+
+        REQUIRE(same_values(-f, minus_values));
+        REQUIRE(same_values(-const_f, minus_values));
+      }
+    }
+
+    SECTION("3D")
+    {
+      const size_t size    = 2;
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();
+
+      constexpr size_t Dimension = 3;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      SECTION("unary minus")
+      {
+        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            const double x = xj[cell_id][0];
+            const double y = xj[cell_id][1];
+            const double z = xj[cell_id][2];
+            for (size_t i = 0; i < size; ++i) {
+              f[cell_id][i] = 2 * x + i * y - z;
+            }
+          });
+
+        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+
+        Table<double> minus_values{mesh->numberOfCells(), size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            for (size_t i = 0; i < size; ++i) {
+              minus_values[cell_id][i] = -f[cell_id][i];
+            }
+          });
+
+        REQUIRE(same_values(-f, minus_values));
+        REQUIRE(same_values(-const_f, minus_values));
+      }
+    }
+  }
+
+  SECTION("binary operators")
+  {
+    SECTION("1D")
+    {
+      const size_t size = 3;
+
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();
+
+      constexpr size_t Dimension = 1;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      SECTION("inner operators")
+      {
+        SECTION("scalar functions")
+        {
+          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              f[cell_id][0]  = 2 * x + 1;
+              f[cell_id][1]  = x * x - 1;
+              f[cell_id][2]  = 2 + x;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              g[cell_id][0]  = (x + 1) * (x - 2) + 1;
+              g[cell_id][1]  = 3 * (x + 2) - 1;
+              g[cell_id][2]  = (x + 3) * 5;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};
+
+          SECTION("sum")
+          {
+            Table<double> sum_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  sum_values[cell_id][i] = f[cell_id][i] + g[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(f + g, sum_values));
+            REQUIRE(same_values(const_f + g, sum_values));
+            REQUIRE(same_values(f + const_g, sum_values));
+            REQUIRE(same_values(const_f + const_g, sum_values));
+          }
+
+          SECTION("difference")
+          {
+            Table<double> difference_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  difference_values[cell_id][i] = f[cell_id][i] - g[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(f - g, difference_values));
+            REQUIRE(same_values(const_f - g, difference_values));
+            REQUIRE(same_values(f - const_g, difference_values));
+            REQUIRE(same_values(const_f - const_g, difference_values));
+          }
+        }
+      }
+
+      SECTION("external operators")
+      {
+        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            const double x = xj[cell_id][0];
+            for (size_t i = 0; i < size; ++i) {
+              f[cell_id][i] = std::abs(2 * x) + i;
+            }
+          });
+
+        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+
+        SECTION("product")
+        {
+          SECTION("scalar lhs")
+          {
+            const double a = 3.2;
+            Table<double> product_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  product_values[cell_id][i] = a * f[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(a * f, product_values));
+            REQUIRE(same_values(a * const_f, product_values));
+          }
+
+          SECTION("DiscreteFunctionP0 lhs")
+          {
+            DiscreteFunctionP0<Dimension, double> a{mesh};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                const double x = xj[cell_id][0];
+                a[cell_id]     = 2 * x + 1;
+              });
+
+            Table<double> product_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  product_values[cell_id][i] = a[cell_id] * f[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(a * f, product_values));
+            REQUIRE(same_values(a * const_f, product_values));
+
+            DiscreteFunctionP0<Dimension, const double> const_a = a;
+            REQUIRE(same_values(const_a * f, product_values));
+            REQUIRE(same_values(const_a * const_f, product_values));
+          }
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      const size_t size = 3;
+
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();
+
+      constexpr size_t Dimension = 2;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      SECTION("inner operators")
+      {
+        SECTION("scalar functions")
+        {
+          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const double y = xj[cell_id][1];
+              f[cell_id][0]  = 2 * x + 1;
+              f[cell_id][1]  = x * x - y;
+              f[cell_id][2]  = 2 + x * y;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const double y = xj[cell_id][1];
+              g[cell_id][0]  = (x + 1) * (y - 2) + 1;
+              g[cell_id][1]  = 3 * (x + 2) - y;
+              g[cell_id][2]  = (x + 3) + 5 * y;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};
+
+          SECTION("sum")
+          {
+            Table<double> sum_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  sum_values[cell_id][i] = f[cell_id][i] + g[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(f + g, sum_values));
+            REQUIRE(same_values(const_f + g, sum_values));
+            REQUIRE(same_values(f + const_g, sum_values));
+            REQUIRE(same_values(const_f + const_g, sum_values));
+          }
+
+          SECTION("difference")
+          {
+            Table<double> difference_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  difference_values[cell_id][i] = f[cell_id][i] - g[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(f - g, difference_values));
+            REQUIRE(same_values(const_f - g, difference_values));
+            REQUIRE(same_values(f - const_g, difference_values));
+            REQUIRE(same_values(const_f - const_g, difference_values));
+          }
+        }
+      }
+
+      SECTION("external operators")
+      {
+        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            const double x = xj[cell_id][0];
+            const double y = xj[cell_id][1];
+            for (size_t i = 0; i < size; ++i) {
+              f[cell_id][i] = std::abs(2 * x) + i * y;
+            }
+          });
+
+        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+
+        SECTION("product")
+        {
+          SECTION("scalar lhs")
+          {
+            const double a = 3.2;
+            Table<double> product_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  product_values[cell_id][i] = a * f[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(a * f, product_values));
+            REQUIRE(same_values(a * const_f, product_values));
+          }
+
+          SECTION("DiscreteFunctionP0 lhs")
+          {
+            DiscreteFunctionP0<Dimension, double> a{mesh};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                const double x = xj[cell_id][0];
+                const double y = xj[cell_id][1];
+                a[cell_id]     = 2 * x + 1 - y;
+              });
+
+            Table<double> product_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  product_values[cell_id][i] = a[cell_id] * f[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(a * f, product_values));
+            REQUIRE(same_values(a * const_f, product_values));
+
+            DiscreteFunctionP0<Dimension, const double> const_a = a;
+            REQUIRE(same_values(const_a * f, product_values));
+            REQUIRE(same_values(const_a * const_f, product_values));
+          }
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      const size_t size = 2;
+
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();
+
+      constexpr size_t Dimension = 3;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      SECTION("inner operators")
+      {
+        SECTION("scalar functions")
+        {
+          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const double y = xj[cell_id][1];
+              const double z = xj[cell_id][2];
+              f[cell_id][0]  = 2 * x * z + 1;
+              f[cell_id][1]  = x * z - y;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const double y = xj[cell_id][1];
+              const double z = xj[cell_id][2];
+              g[cell_id][0]  = (x + 1) * (y - 2) + 1 - z;
+              g[cell_id][1]  = 3 * (x + 2) - y * z;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};
+
+          SECTION("sum")
+          {
+            Table<double> sum_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  sum_values[cell_id][i] = f[cell_id][i] + g[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(f + g, sum_values));
+            REQUIRE(same_values(const_f + g, sum_values));
+            REQUIRE(same_values(f + const_g, sum_values));
+            REQUIRE(same_values(const_f + const_g, sum_values));
+          }
+
+          SECTION("difference")
+          {
+            Table<double> difference_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  difference_values[cell_id][i] = f[cell_id][i] - g[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(f - g, difference_values));
+            REQUIRE(same_values(const_f - g, difference_values));
+            REQUIRE(same_values(f - const_g, difference_values));
+            REQUIRE(same_values(const_f - const_g, difference_values));
+          }
+        }
+      }
+
+      SECTION("external operators")
+      {
+        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            const double x = xj[cell_id][0];
+            const double y = xj[cell_id][1];
+            const double z = xj[cell_id][2];
+            for (size_t i = 0; i < size; ++i) {
+              f[cell_id][i] = std::abs(2 * x) + i * y + z;
+            }
+          });
+
+        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+
+        SECTION("product")
+        {
+          SECTION("scalar lhs")
+          {
+            const double a = 3.2;
+            Table<double> product_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  product_values[cell_id][i] = a * f[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(a * f, product_values));
+            REQUIRE(same_values(a * const_f, product_values));
+          }
+
+          SECTION("DiscreteFunctionP0 lhs")
+          {
+            DiscreteFunctionP0<Dimension, double> a{mesh};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                const double x = xj[cell_id][0];
+                const double y = xj[cell_id][1];
+                const double z = xj[cell_id][2];
+                a[cell_id]     = 2 * x + 1 - y * z;
+              });
+
+            Table<double> product_values{mesh->numberOfCells(), size};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                for (size_t i = 0; i < size; ++i) {
+                  product_values[cell_id][i] = a[cell_id] * f[cell_id][i];
+                }
+              });
+
+            REQUIRE(same_values(a * f, product_values));
+            REQUIRE(same_values(a * const_f, product_values));
+
+            DiscreteFunctionP0<Dimension, const double> const_a = a;
+            REQUIRE(same_values(const_a * f, product_values));
+            REQUIRE(same_values(const_a * const_f, product_values));
+          }
+        }
+      }
+    }
+  }
+}