diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index a48a72fba851266c09ede6192c6b4e658bac97fe..61c71097381797b42d7526abc1be48003ef765ca 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -101,12 +101,13 @@ add_executable (unit_tests
 
 add_executable (mpi_unit_tests
   mpi_test_main.cpp
-  test_Messenger.cpp
-  test_Partitioner.cpp
+  test_DiscreteFunctionP0.cpp
   test_ItemArray.cpp
   test_ItemArrayUtils.cpp
   test_ItemValue.cpp
   test_ItemValueUtils.cpp
+  test_Messenger.cpp
+  test_Partitioner.cpp
   test_SubItemValuePerItem.cpp
   test_SubItemArrayPerItem.cpp
   )
diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e445d86d9b30683b86a886be305cd1d3d64e5806
--- /dev/null
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -0,0 +1,1847 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionP0", "[scheme]")
+{
+  auto same_values = [](const auto& f, const auto& g) {
+    size_t number_of_cells = f.cellValues().numberOfItems();
+    for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+      if (f[cell_id] != g[cell_id]) {
+        return false;
+      }
+    }
+    return true;
+  };
+
+  SECTION("constructors")
+  {
+    SECTION("1D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();
+
+      constexpr size_t Dimension = 1;
+
+      DiscreteFunctionP0<Dimension, double> f{mesh};
+      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0);
+
+      REQUIRE(f.mesh().get() == mesh.get());
+
+      DiscreteFunctionP0 g{f};
+      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0);
+
+      CellValue<TinyVector<Dimension>> h_values{mesh->connectivity()};
+      h_values.fill(ZeroType{});
+
+      DiscreteFunctionP0 zero{mesh, [&] {
+                                CellValue<TinyVector<Dimension>> cell_value{mesh->connectivity()};
+                                cell_value.fill(ZeroType{});
+                                return cell_value;
+                              }()};
+
+      DiscreteFunctionP0 h{mesh, h_values};
+      REQUIRE(same_values(h, zero));
+      REQUIRE(same_values(h, h_values));
+
+      copy_to(MeshDataManager::instance().getMeshData(*mesh).xj(), h_values);
+
+      REQUIRE(same_values(h, h_values));
+      REQUIRE(not same_values(h, zero));
+
+      DiscreteFunctionP0 moved_h{std::move(h)};
+      REQUIRE(same_values(moved_h, h_values));
+    }
+
+    SECTION("2D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();
+
+      constexpr size_t Dimension = 2;
+
+      DiscreteFunctionP0<Dimension, double> f{mesh};
+      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0);
+
+      REQUIRE(f.mesh().get() == mesh.get());
+
+      DiscreteFunctionP0 g{f};
+      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0);
+
+      CellValue<TinyVector<Dimension>> h_values{mesh->connectivity()};
+      h_values.fill(ZeroType{});
+
+      DiscreteFunctionP0 zero{mesh, [&] {
+                                CellValue<TinyVector<Dimension>> cell_value{mesh->connectivity()};
+                                cell_value.fill(ZeroType{});
+                                return cell_value;
+                              }()};
+
+      DiscreteFunctionP0 h{mesh, h_values};
+      REQUIRE(same_values(h, zero));
+      REQUIRE(same_values(h, h_values));
+
+      copy_to(MeshDataManager::instance().getMeshData(*mesh).xj(), h_values);
+
+      REQUIRE(same_values(h, h_values));
+      REQUIRE(not same_values(h, zero));
+
+      DiscreteFunctionP0 moved_h{std::move(h)};
+      REQUIRE(same_values(moved_h, h_values));
+    }
+
+    SECTION("3D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();
+
+      constexpr size_t Dimension = 3;
+
+      DiscreteFunctionP0<Dimension, double> f{mesh};
+      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0);
+
+      REQUIRE(f.mesh().get() == mesh.get());
+
+      DiscreteFunctionP0 g{f};
+      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
+      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0);
+
+      CellValue<TinyVector<Dimension>> h_values{mesh->connectivity()};
+      h_values.fill(ZeroType{});
+
+      DiscreteFunctionP0 zero{mesh, [&] {
+                                CellValue<TinyVector<Dimension>> cell_value{mesh->connectivity()};
+                                cell_value.fill(ZeroType{});
+                                return cell_value;
+                              }()};
+
+      DiscreteFunctionP0 h{mesh, h_values};
+      REQUIRE(same_values(h, zero));
+      REQUIRE(same_values(h, h_values));
+
+      copy_to(MeshDataManager::instance().getMeshData(*mesh).xj(), h_values);
+
+      REQUIRE(same_values(h, h_values));
+      REQUIRE(not same_values(h, zero));
+
+      DiscreteFunctionP0 moved_h{std::move(h)};
+      REQUIRE(same_values(moved_h, h_values));
+    }
+  }
+
+  SECTION("binary operators")
+  {
+    SECTION("1D")
+    {
+      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")
+        {
+          DiscreteFunctionP0<Dimension, double> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              f[cell_id]     = 2 * x + 1;
+            });
+
+          DiscreteFunctionP0<Dimension, double> g{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              g[cell_id]     = std::abs((x + 1) * (x - 2)) + 1;
+            });
+
+          DiscreteFunctionP0<Dimension, const double> const_f = f;
+          DiscreteFunctionP0<Dimension, const double> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<double> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<double> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("product")
+          {
+            Array<double> product_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * g[cell_id]; });
+
+            REQUIRE(same_values(f * g, product_values));
+            REQUIRE(same_values(const_f * g, product_values));
+            REQUIRE(same_values(f * const_g, product_values));
+            REQUIRE(same_values(const_f * const_g, product_values));
+          }
+
+          SECTION("ratio")
+          {
+            Array<double> ratio_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = f[cell_id] / g[cell_id]; });
+
+            REQUIRE(same_values(f / g, ratio_values));
+            REQUIRE(same_values(const_f / g, ratio_values));
+            REQUIRE(same_values(f / const_g, ratio_values));
+            REQUIRE(same_values(const_f / const_g, ratio_values));
+          }
+        }
+
+        SECTION("vector functions")
+        {
+          constexpr std::uint64_t VectorDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyVector<VectorDimension> X{x, 2 - x};
+              f[cell_id] = 2 * X + TinyVector<2>{1, 2};
+            });
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> g{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyVector<VectorDimension> X{3 * x + 1, 2 + x};
+              g[cell_id] = X;
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_f = f;
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("matrix functions")
+        {
+          constexpr std::uint64_t MatrixDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyMatrix<MatrixDimension> A{x, 2 - x, 2 * x, x * x - 3};
+              f[cell_id] = 2 * A + TinyMatrix<2>{1, 2, 3, 4};
+            });
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> g{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyMatrix<MatrixDimension> A{3 * x + 1, 2 + x, 1 - 2 * x, 2 * x * x};
+              g[cell_id] = A;
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_f = f;
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("product")
+          {
+            Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * g[cell_id]; });
+
+            REQUIRE(same_values(f * g, product_values));
+            REQUIRE(same_values(const_f * g, product_values));
+            REQUIRE(same_values(f * const_g, product_values));
+            REQUIRE(same_values(const_f * const_g, product_values));
+          }
+        }
+      }
+
+      SECTION("external operators")
+      {
+        SECTION("scalar functions")
+        {
+          DiscreteFunctionP0<Dimension, double> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              f[cell_id]     = std::abs(2 * x) + 1;
+            });
+
+          const double a = 3;
+
+          DiscreteFunctionP0<Dimension, const double> const_f = f;
+
+          SECTION("sum")
+          {
+            {
+              Array<double> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = a + f[cell_id]; });
+
+              REQUIRE(same_values(a + f, sum_values));
+              REQUIRE(same_values(a + const_f, sum_values));
+            }
+            {
+              Array<double> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + a; });
+
+              REQUIRE(same_values(f + a, sum_values));
+              REQUIRE(same_values(const_f + a, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            {
+              Array<double> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = a - f[cell_id]; });
+              REQUIRE(same_values(a - f, difference_values));
+              REQUIRE(same_values(a - const_f, difference_values));
+            }
+
+            {
+              Array<double> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - a; });
+              REQUIRE(same_values(f - a, difference_values));
+              REQUIRE(same_values(const_f - a, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              Array<double> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+            {
+              Array<double> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * a; });
+
+              REQUIRE(same_values(f * a, product_values));
+              REQUIRE(same_values(const_f * a, product_values));
+            }
+
+            {
+              Array<TinyVector<3>> product_values{mesh->numberOfCells()};
+              const TinyVector<3> v{1, 2, 3};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * v; });
+
+              REQUIRE(same_values(f * v, product_values));
+              REQUIRE(same_values(const_f * v, product_values));
+            }
+
+            {
+              Array<TinyVector<3>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyVector<3>> v{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  v[cell_id]     = TinyVector<3>{x, 2 * x, 1 - x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * v[cell_id]; });
+
+              REQUIRE(same_values(f * v, product_values));
+              REQUIRE(same_values(const_f * v, product_values));
+            }
+
+            {
+              Array<TinyMatrix<2>> product_values{mesh->numberOfCells()};
+              const TinyMatrix<2> A{1, 2, 3, 4};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * A; });
+
+              REQUIRE(same_values(f * A, product_values));
+              REQUIRE(same_values(const_f * A, product_values));
+            }
+
+            {
+              Array<TinyMatrix<2>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyMatrix<2>> M{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  M[cell_id]     = TinyMatrix<2>{x, 2 * x, 1 - x, 2 - x * x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * M[cell_id]; });
+
+              REQUIRE(same_values(f * M, product_values));
+              REQUIRE(same_values(const_f * M, product_values));
+            }
+          }
+
+          SECTION("ratio")
+          {
+            {
+              Array<double> ratio_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = a / f[cell_id]; });
+
+              REQUIRE(same_values(a / f, ratio_values));
+              REQUIRE(same_values(a / const_f, ratio_values));
+            }
+            {
+              Array<double> ratio_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = f[cell_id] / a; });
+
+              REQUIRE(same_values(f / a, ratio_values));
+              REQUIRE(same_values(const_f / a, ratio_values));
+            }
+          }
+        }
+
+        SECTION("vector functions")
+        {
+          constexpr std::uint64_t VectorDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyVector<VectorDimension> X{x, 2 - x};
+              f[cell_id] = 2 * X + TinyVector<2>{1, 2};
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_f = f;
+
+          SECTION("sum")
+          {
+            const TinyVector<VectorDimension> v{1, 2};
+            {
+              Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = v + f[cell_id]; });
+
+              REQUIRE(same_values(v + f, sum_values));
+              REQUIRE(same_values(v + const_f, sum_values));
+            }
+            {
+              Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + v; });
+
+              REQUIRE(same_values(f + v, sum_values));
+              REQUIRE(same_values(const_f + v, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            const TinyVector<VectorDimension> v{1, 2};
+            {
+              Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = v - f[cell_id]; });
+
+              REQUIRE(same_values(v - f, difference_values));
+              REQUIRE(same_values(v - const_f, difference_values));
+            }
+            {
+              Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - v; });
+
+              REQUIRE(same_values(f - v, difference_values));
+              REQUIRE(same_values(const_f - v, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              const double a = 2.3;
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              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 * x - 1;
+                });
+
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<VectorDimension> A{1, 2, 3, 4};
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = A * f[cell_id]; });
+
+              REQUIRE(same_values(A * f, product_values));
+              REQUIRE(same_values(A * const_f, product_values));
+            }
+
+            {
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyMatrix<VectorDimension>> M{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  M[cell_id]     = TinyMatrix<2>{x, 2 * x, 1 - x, 2 - x * x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = M[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(M * f, product_values));
+              REQUIRE(same_values(M * const_f, product_values));
+            }
+          }
+        }
+
+        SECTION("matrix functions")
+        {
+          constexpr std::uint64_t MatrixDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyMatrix<MatrixDimension> X{x, 2 - x, x * x, x * 3};
+              f[cell_id] = 2 * X + TinyMatrix<2>{1, 2, 3, 4};
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_f = f;
+
+          SECTION("sum")
+          {
+            const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+            {
+              Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = A + f[cell_id]; });
+
+              REQUIRE(same_values(A + f, sum_values));
+              REQUIRE(same_values(A + const_f, sum_values));
+            }
+            {
+              Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + A; });
+
+              REQUIRE(same_values(f + A, sum_values));
+              REQUIRE(same_values(const_f + A, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+            {
+              Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = A - f[cell_id]; });
+
+              REQUIRE(same_values(A - f, difference_values));
+              REQUIRE(same_values(A - const_f, difference_values));
+            }
+            {
+              Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - A; });
+
+              REQUIRE(same_values(f - A, difference_values));
+              REQUIRE(same_values(const_f - A, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              const double a = 2.3;
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              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 * x - 1;
+                });
+
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = A * f[cell_id]; });
+
+              REQUIRE(same_values(A * f, product_values));
+              REQUIRE(same_values(A * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * A; });
+
+              REQUIRE(same_values(f * A, product_values));
+              REQUIRE(same_values(const_f * A, product_values));
+            }
+          }
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      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")
+        {
+          DiscreteFunctionP0<Dimension, double> f{mesh};
+          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]     = 2 * x + y + 1;
+            });
+
+          DiscreteFunctionP0<Dimension, double> g{mesh};
+          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]     = std::abs((x + 1) * (x - 2) + y * (1 + y)) + 1;
+            });
+
+          DiscreteFunctionP0<Dimension, const double> const_f = f;
+          DiscreteFunctionP0<Dimension, const double> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<double> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<double> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("product")
+          {
+            Array<double> product_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * g[cell_id]; });
+
+            REQUIRE(same_values(f * g, product_values));
+            REQUIRE(same_values(const_f * g, product_values));
+            REQUIRE(same_values(f * const_g, product_values));
+            REQUIRE(same_values(const_f * const_g, product_values));
+          }
+
+          SECTION("ratio")
+          {
+            Array<double> ratio_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = f[cell_id] / g[cell_id]; });
+
+            REQUIRE(same_values(f / g, ratio_values));
+            REQUIRE(same_values(const_f / g, ratio_values));
+            REQUIRE(same_values(f / const_g, ratio_values));
+            REQUIRE(same_values(const_f / const_g, ratio_values));
+          }
+        }
+
+        SECTION("vector functions")
+        {
+          constexpr std::uint64_t VectorDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyVector<VectorDimension> X{x, 2 - x};
+              f[cell_id] = 2 * X + TinyVector<2>{1, 2};
+            });
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> g{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyVector<VectorDimension> X{3 * x + 1, 2 + x};
+              g[cell_id] = X;
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_f = f;
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("matrix functions")
+        {
+          constexpr std::uint64_t MatrixDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyMatrix<MatrixDimension> A{x, 2 - x, 2 * x, x * x - 3};
+              f[cell_id] = 2 * A + TinyMatrix<2>{1, 2, 3, 4};
+            });
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> g{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyMatrix<MatrixDimension> A{3 * x + 1, 2 + x, 1 - 2 * x, 2 * x * x};
+              g[cell_id] = A;
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_f = f;
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("product")
+          {
+            Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * g[cell_id]; });
+
+            REQUIRE(same_values(f * g, product_values));
+            REQUIRE(same_values(const_f * g, product_values));
+            REQUIRE(same_values(f * const_g, product_values));
+            REQUIRE(same_values(const_f * const_g, product_values));
+          }
+        }
+      }
+
+      SECTION("external operators")
+      {
+        SECTION("scalar functions")
+        {
+          DiscreteFunctionP0<Dimension, double> f{mesh};
+          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]     = std::abs(2 * x + y) + 1;
+            });
+
+          const double a = 3;
+
+          DiscreteFunctionP0<Dimension, const double> const_f = f;
+
+          SECTION("sum")
+          {
+            {
+              Array<double> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = a + f[cell_id]; });
+
+              REQUIRE(same_values(a + f, sum_values));
+              REQUIRE(same_values(a + const_f, sum_values));
+            }
+            {
+              Array<double> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + a; });
+
+              REQUIRE(same_values(f + a, sum_values));
+              REQUIRE(same_values(const_f + a, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            {
+              Array<double> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = a - f[cell_id]; });
+              REQUIRE(same_values(a - f, difference_values));
+              REQUIRE(same_values(a - const_f, difference_values));
+            }
+
+            {
+              Array<double> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - a; });
+              REQUIRE(same_values(f - a, difference_values));
+              REQUIRE(same_values(const_f - a, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              Array<double> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+            {
+              Array<double> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * a; });
+
+              REQUIRE(same_values(f * a, product_values));
+              REQUIRE(same_values(const_f * a, product_values));
+            }
+
+            {
+              Array<TinyVector<3>> product_values{mesh->numberOfCells()};
+              const TinyVector<3> v{1, 2, 3};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * v; });
+
+              REQUIRE(same_values(f * v, product_values));
+              REQUIRE(same_values(const_f * v, product_values));
+            }
+
+            {
+              Array<TinyVector<3>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyVector<3>> v{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  v[cell_id]     = TinyVector<3>{x, 2 * x, 1 - x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * v[cell_id]; });
+
+              REQUIRE(same_values(f * v, product_values));
+              REQUIRE(same_values(const_f * v, product_values));
+            }
+
+            {
+              Array<TinyMatrix<2>> product_values{mesh->numberOfCells()};
+              const TinyMatrix<2> A{1, 2, 3, 4};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * A; });
+
+              REQUIRE(same_values(f * A, product_values));
+              REQUIRE(same_values(const_f * A, product_values));
+            }
+
+            {
+              Array<TinyMatrix<2>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyMatrix<2>> M{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  M[cell_id]     = TinyMatrix<2>{x, 2 * x, 1 - x, 2 - x * x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * M[cell_id]; });
+
+              REQUIRE(same_values(f * M, product_values));
+              REQUIRE(same_values(const_f * M, product_values));
+            }
+          }
+
+          SECTION("ratio")
+          {
+            {
+              Array<double> ratio_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = a / f[cell_id]; });
+
+              REQUIRE(same_values(a / f, ratio_values));
+              REQUIRE(same_values(a / const_f, ratio_values));
+            }
+            {
+              Array<double> ratio_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = f[cell_id] / a; });
+
+              REQUIRE(same_values(f / a, ratio_values));
+              REQUIRE(same_values(const_f / a, ratio_values));
+            }
+          }
+        }
+
+        SECTION("vector functions")
+        {
+          constexpr std::uint64_t VectorDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> f{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 TinyVector<VectorDimension> X{x + y, 2 - x * y};
+              f[cell_id] = 2 * X + TinyVector<2>{1, 2};
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_f = f;
+
+          SECTION("sum")
+          {
+            const TinyVector<VectorDimension> v{1, 2};
+            {
+              Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = v + f[cell_id]; });
+
+              REQUIRE(same_values(v + f, sum_values));
+              REQUIRE(same_values(v + const_f, sum_values));
+            }
+            {
+              Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + v; });
+
+              REQUIRE(same_values(f + v, sum_values));
+              REQUIRE(same_values(const_f + v, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            const TinyVector<VectorDimension> v{1, 2};
+            {
+              Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = v - f[cell_id]; });
+
+              REQUIRE(same_values(v - f, difference_values));
+              REQUIRE(same_values(v - const_f, difference_values));
+            }
+            {
+              Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - v; });
+
+              REQUIRE(same_values(f - v, difference_values));
+              REQUIRE(same_values(const_f - v, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              const double a = 2.3;
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              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 * x - 1;
+                });
+
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<VectorDimension> A{1, 2, 3, 4};
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = A * f[cell_id]; });
+
+              REQUIRE(same_values(A * f, product_values));
+              REQUIRE(same_values(A * const_f, product_values));
+            }
+
+            {
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyMatrix<VectorDimension>> M{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  M[cell_id]     = TinyMatrix<2>{x, 2 * x, 1 - x, 2 - x * x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = M[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(M * f, product_values));
+              REQUIRE(same_values(M * const_f, product_values));
+            }
+          }
+        }
+
+        SECTION("matrix functions")
+        {
+          constexpr std::uint64_t MatrixDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> f{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 TinyMatrix<MatrixDimension> X{x, 2 - y, x * y, y * 3};
+              f[cell_id] = 2 * X + TinyMatrix<2>{1, 2, 3, 4};
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_f = f;
+
+          SECTION("sum")
+          {
+            const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+            {
+              Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = A + f[cell_id]; });
+
+              REQUIRE(same_values(A + f, sum_values));
+              REQUIRE(same_values(A + const_f, sum_values));
+            }
+            {
+              Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + A; });
+
+              REQUIRE(same_values(f + A, sum_values));
+              REQUIRE(same_values(const_f + A, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+            {
+              Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = A - f[cell_id]; });
+
+              REQUIRE(same_values(A - f, difference_values));
+              REQUIRE(same_values(A - const_f, difference_values));
+            }
+            {
+              Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - A; });
+
+              REQUIRE(same_values(f - A, difference_values));
+              REQUIRE(same_values(const_f - A, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              const double a = 2.3;
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              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 * x - 1;
+                });
+
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = A * f[cell_id]; });
+
+              REQUIRE(same_values(A * f, product_values));
+              REQUIRE(same_values(A * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * A; });
+
+              REQUIRE(same_values(f * A, product_values));
+              REQUIRE(same_values(const_f * A, product_values));
+            }
+          }
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      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")
+        {
+          DiscreteFunctionP0<Dimension, double> f{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];
+              f[cell_id]     = 2 * x + y - z;
+            });
+
+          DiscreteFunctionP0<Dimension, double> g{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];
+              g[cell_id]     = std::abs((x + 1) * (x - 2) + y * (1 + y) + 2 * z) + 1;
+            });
+
+          DiscreteFunctionP0<Dimension, const double> const_f = f;
+          DiscreteFunctionP0<Dimension, const double> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<double> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<double> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("product")
+          {
+            Array<double> product_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * g[cell_id]; });
+
+            REQUIRE(same_values(f * g, product_values));
+            REQUIRE(same_values(const_f * g, product_values));
+            REQUIRE(same_values(f * const_g, product_values));
+            REQUIRE(same_values(const_f * const_g, product_values));
+          }
+
+          SECTION("ratio")
+          {
+            Array<double> ratio_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = f[cell_id] / g[cell_id]; });
+
+            REQUIRE(same_values(f / g, ratio_values));
+            REQUIRE(same_values(const_f / g, ratio_values));
+            REQUIRE(same_values(f / const_g, ratio_values));
+            REQUIRE(same_values(const_f / const_g, ratio_values));
+          }
+        }
+
+        SECTION("vector functions")
+        {
+          constexpr std::uint64_t VectorDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyVector<VectorDimension> X{x, 2 - x};
+              f[cell_id] = 2 * X + TinyVector<2>{1, 2};
+            });
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> g{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyVector<VectorDimension> X{3 * x + 1, 2 + x};
+              g[cell_id] = X;
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_f = f;
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("matrix functions")
+        {
+          constexpr std::uint64_t MatrixDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> f{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyMatrix<MatrixDimension> A{x, 2 - x, 2 * x, x * x - 3};
+              f[cell_id] = 2 * A + TinyMatrix<2>{1, 2, 3, 4};
+            });
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> g{mesh};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              const TinyMatrix<MatrixDimension> A{3 * x + 1, 2 + x, 1 - 2 * x, 2 * x * x};
+              g[cell_id] = A;
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_f = f;
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_g{g};
+
+          SECTION("sum")
+          {
+            Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + g[cell_id]; });
+
+            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")
+          {
+            Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - g[cell_id]; });
+
+            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("product")
+          {
+            Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(),
+              PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * g[cell_id]; });
+
+            REQUIRE(same_values(f * g, product_values));
+            REQUIRE(same_values(const_f * g, product_values));
+            REQUIRE(same_values(f * const_g, product_values));
+            REQUIRE(same_values(const_f * const_g, product_values));
+          }
+        }
+      }
+
+      SECTION("external operators")
+      {
+        SECTION("scalar functions")
+        {
+          DiscreteFunctionP0<Dimension, double> f{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];
+              f[cell_id]     = std::abs(2 * x + y * z) + 1;
+            });
+
+          const double a = 3;
+
+          DiscreteFunctionP0<Dimension, const double> const_f = f;
+
+          SECTION("sum")
+          {
+            {
+              Array<double> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = a + f[cell_id]; });
+
+              REQUIRE(same_values(a + f, sum_values));
+              REQUIRE(same_values(a + const_f, sum_values));
+            }
+            {
+              Array<double> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + a; });
+
+              REQUIRE(same_values(f + a, sum_values));
+              REQUIRE(same_values(const_f + a, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            {
+              Array<double> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = a - f[cell_id]; });
+              REQUIRE(same_values(a - f, difference_values));
+              REQUIRE(same_values(a - const_f, difference_values));
+            }
+
+            {
+              Array<double> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - a; });
+              REQUIRE(same_values(f - a, difference_values));
+              REQUIRE(same_values(const_f - a, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              Array<double> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+            {
+              Array<double> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * a; });
+
+              REQUIRE(same_values(f * a, product_values));
+              REQUIRE(same_values(const_f * a, product_values));
+            }
+
+            {
+              Array<TinyVector<3>> product_values{mesh->numberOfCells()};
+              const TinyVector<3> v{1, 2, 3};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * v; });
+
+              REQUIRE(same_values(f * v, product_values));
+              REQUIRE(same_values(const_f * v, product_values));
+            }
+
+            {
+              Array<TinyVector<3>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyVector<3>> v{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  v[cell_id]     = TinyVector<3>{x, 2 * x, 1 - x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * v[cell_id]; });
+
+              REQUIRE(same_values(f * v, product_values));
+              REQUIRE(same_values(const_f * v, product_values));
+            }
+
+            {
+              Array<TinyMatrix<2>> product_values{mesh->numberOfCells()};
+              const TinyMatrix<2> A{1, 2, 3, 4};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * A; });
+
+              REQUIRE(same_values(f * A, product_values));
+              REQUIRE(same_values(const_f * A, product_values));
+            }
+
+            {
+              Array<TinyMatrix<2>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyMatrix<2>> M{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  M[cell_id]     = TinyMatrix<2>{x, 2 * x, 1 - x, 2 - x * x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * M[cell_id]; });
+
+              REQUIRE(same_values(f * M, product_values));
+              REQUIRE(same_values(const_f * M, product_values));
+            }
+          }
+
+          SECTION("ratio")
+          {
+            {
+              Array<double> ratio_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = a / f[cell_id]; });
+
+              REQUIRE(same_values(a / f, ratio_values));
+              REQUIRE(same_values(a / const_f, ratio_values));
+            }
+            {
+              Array<double> ratio_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio_values[cell_id] = f[cell_id] / a; });
+
+              REQUIRE(same_values(f / a, ratio_values));
+              REQUIRE(same_values(const_f / a, ratio_values));
+            }
+          }
+        }
+
+        SECTION("vector functions")
+        {
+          constexpr std::uint64_t VectorDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>> f{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];
+              const TinyVector<VectorDimension> X{x + y - z, 2 - x * y};
+              f[cell_id] = 2 * X + TinyVector<2>{1, 2};
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyVector<VectorDimension>> const_f = f;
+
+          SECTION("sum")
+          {
+            const TinyVector<VectorDimension> v{1, 2};
+            {
+              Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = v + f[cell_id]; });
+
+              REQUIRE(same_values(v + f, sum_values));
+              REQUIRE(same_values(v + const_f, sum_values));
+            }
+            {
+              Array<TinyVector<VectorDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + v; });
+
+              REQUIRE(same_values(f + v, sum_values));
+              REQUIRE(same_values(const_f + v, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            const TinyVector<VectorDimension> v{1, 2};
+            {
+              Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = v - f[cell_id]; });
+
+              REQUIRE(same_values(v - f, difference_values));
+              REQUIRE(same_values(v - const_f, difference_values));
+            }
+            {
+              Array<TinyVector<VectorDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - v; });
+
+              REQUIRE(same_values(f - v, difference_values));
+              REQUIRE(same_values(const_f - v, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              const double a = 2.3;
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              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 * x - 1;
+                });
+
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<VectorDimension> A{1, 2, 3, 4};
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = A * f[cell_id]; });
+
+              REQUIRE(same_values(A * f, product_values));
+              REQUIRE(same_values(A * const_f, product_values));
+            }
+
+            {
+              Array<TinyVector<VectorDimension>> product_values{mesh->numberOfCells()};
+              DiscreteFunctionP0<Dimension, TinyMatrix<VectorDimension>> M{mesh};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                  const double x = xj[cell_id][0];
+                  M[cell_id]     = TinyMatrix<2>{x, 2 * x, 1 - x, 2 - x * x};
+                });
+
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = M[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(M * f, product_values));
+              REQUIRE(same_values(M * const_f, product_values));
+            }
+          }
+        }
+
+        SECTION("matrix functions")
+        {
+          constexpr std::uint64_t MatrixDimension = 2;
+
+          DiscreteFunctionP0<Dimension, TinyMatrix<MatrixDimension>> f{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];
+              const TinyMatrix<MatrixDimension> X{x, 2 - y, x * y, y * z + 3};
+              f[cell_id] = 2 * X + TinyMatrix<2>{1, 2, 3, 4};
+            });
+
+          DiscreteFunctionP0<Dimension, const TinyMatrix<MatrixDimension>> const_f = f;
+
+          SECTION("sum")
+          {
+            const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+            {
+              Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = A + f[cell_id]; });
+
+              REQUIRE(same_values(A + f, sum_values));
+              REQUIRE(same_values(A + const_f, sum_values));
+            }
+            {
+              Array<TinyMatrix<MatrixDimension>> sum_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum_values[cell_id] = f[cell_id] + A; });
+
+              REQUIRE(same_values(f + A, sum_values));
+              REQUIRE(same_values(const_f + A, sum_values));
+            }
+          }
+
+          SECTION("difference")
+          {
+            const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+            {
+              Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = A - f[cell_id]; });
+
+              REQUIRE(same_values(A - f, difference_values));
+              REQUIRE(same_values(A - const_f, difference_values));
+            }
+            {
+              Array<TinyMatrix<MatrixDimension>> difference_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference_values[cell_id] = f[cell_id] - A; });
+
+              REQUIRE(same_values(f - A, difference_values));
+              REQUIRE(same_values(const_f - A, difference_values));
+            }
+          }
+
+          SECTION("product")
+          {
+            {
+              const double a = 2.3;
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              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 * x - 1;
+                });
+
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(),
+                PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = a[cell_id] * f[cell_id]; });
+
+              REQUIRE(same_values(a * f, product_values));
+              REQUIRE(same_values(a * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = A * f[cell_id]; });
+
+              REQUIRE(same_values(A * f, product_values));
+              REQUIRE(same_values(A * const_f, product_values));
+            }
+
+            {
+              const TinyMatrix<MatrixDimension> A{1, 2, 3, 4};
+              Array<TinyMatrix<MatrixDimension>> product_values{mesh->numberOfCells()};
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product_values[cell_id] = f[cell_id] * A; });
+
+              REQUIRE(same_values(f * A, product_values));
+              REQUIRE(same_values(const_f * A, product_values));
+            }
+          }
+        }
+      }
+    }
+  }
+}