diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index 8827de489048f78e7c672e79b5189d5303b66f56..acfa7a0064d5368c251af31ff3cb8b9a6a9a749a 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -2329,9 +2329,10 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
     DiscreteFunctionP0 lhs    = lhs_expression;                                \
     DiscreteFunctionP0 rhs    = rhs_expression;                                \
     DiscreteFunctionP0 result = FCT(lhs, rhs);                                 \
-    bool is_same              = true;                                          \
+    using namespace std;                                                       \
+    bool is_same = true;                                                       \
     parallel_for(lhs.cellValues().numberOfItems(), [&](const CellId cell_id) { \
-      if (result[cell_id] != std::FCT(lhs[cell_id], rhs[cell_id])) {           \
+      if (result[cell_id] != FCT(lhs[cell_id], rhs[cell_id])) {                \
         is_same = false;                                                       \
       }                                                                        \
     });                                                                        \
@@ -2343,8 +2344,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
     DiscreteFunctionP0 rhs    = rhs_expression;                                 \
     DiscreteFunctionP0 result = FCT(lhs, rhs);                                  \
     bool is_same              = true;                                           \
+    using namespace std;                                                        \
     parallel_for(rhs.cellValues().numberOfItems(), [&](const CellId cell_id) {  \
-      if (result[cell_id] != std::FCT(lhs, rhs[cell_id])) {                     \
+      if (result[cell_id] != FCT(lhs, rhs[cell_id])) {                          \
         is_same = false;                                                        \
       }                                                                         \
     });                                                                         \
@@ -2356,8 +2358,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
     DiscreteFunctionP0 lhs    = lhs_expression;                                 \
     DiscreteFunctionP0 result = FCT(lhs, rhs);                                  \
     bool is_same              = true;                                           \
+    using namespace std;                                                        \
     parallel_for(lhs.cellValues().numberOfItems(), [&](const CellId cell_id) {  \
-      if (result[cell_id] != std::FCT(lhs[cell_id], rhs)) {                     \
+      if (result[cell_id] != FCT(lhs[cell_id], rhs)) {                          \
         is_same = false;                                                        \
       }                                                                         \
     });                                                                         \
@@ -2563,6 +2566,813 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
       {
         CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(sin(positive_function), 0.1, max);
       }
+
+      SECTION("dot(uh,hv)")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        DiscreteFunctionP0<Dimension, TinyVector<2>> vh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            vh[cell_id]    = TinyVector<2>{2.3 * x, 1 - x};
+          });
+
+        CHECK_STD_BINARY_MATH_FUNCTION(uh, vh, dot);
+      }
+
+      SECTION("dot(uh,v)")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        const TinyVector<2> v{1, 2};
+
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(uh, v, dot);
+      }
+
+      SECTION("dot(u,hv)")
+      {
+        const TinyVector<2> u{3, -2};
+
+        DiscreteFunctionP0<Dimension, TinyVector<2>> vh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            vh[cell_id]    = TinyVector<2>{2.3 * x, 1 - x};
+          });
+
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, dot);
+      }
+
+      SECTION("scalar sum")
+      {
+        const CellValue<const double> cell_value = positive_function.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(positive_function));
+      }
+
+      SECTION("vector sum")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+        const CellValue<const TinyVector<2>> cell_value = uh.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(uh));
+      }
+
+      SECTION("matrix sum")
+      {
+        DiscreteFunctionP0<Dimension, TinyMatrix<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3, 2 * x, 3 * x - 1};
+          });
+        const CellValue<const TinyMatrix<2>> cell_value = uh.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(uh));
+      }
+
+      SECTION("integrate scalar")
+      {
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<double> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            cell_value[cell_id] = cell_volume[cell_id] * positive_function[cell_id];
+          });
+
+        REQUIRE(integrate(positive_function) == Catch::Approx(sum(cell_value)));
+      }
+
+      SECTION("integrate vector")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<TinyVector<2>> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(),
+          PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_volume[cell_id] * uh[cell_id]; });
+
+        REQUIRE(integrate(uh)[0] == Catch::Approx(sum(cell_value)[0]));
+        REQUIRE(integrate(uh)[1] == Catch::Approx(sum(cell_value)[1]));
+      }
+
+      SECTION("integrate matrix")
+      {
+        DiscreteFunctionP0<Dimension, TinyMatrix<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3, 2 * x, 1 - x};
+          });
+
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<TinyMatrix<2>> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(),
+          PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_volume[cell_id] * uh[cell_id]; });
+
+        REQUIRE(integrate(uh)(0, 0) == Catch::Approx(sum(cell_value)(0, 0)));
+        REQUIRE(integrate(uh)(0, 1) == Catch::Approx(sum(cell_value)(0, 1)));
+        REQUIRE(integrate(uh)(1, 0) == Catch::Approx(sum(cell_value)(1, 0)));
+        REQUIRE(integrate(uh)(1, 1) == Catch::Approx(sum(cell_value)(1, 1)));
+      }
+    }
+
+    SECTION("2D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();
+
+      constexpr size_t Dimension = 2;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      DiscreteFunctionP0<Dimension, double> positive_function{mesh};
+
+      parallel_for(
+        mesh->numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { positive_function[cell_id] = 1 + std::abs(xj[cell_id][0]); });
+
+      const double min_value = min(positive_function);
+      SECTION("min")
+      {
+        double local_min = std::numeric_limits<double>::max();
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          local_min = std::min(local_min, positive_function[cell_id]);
+        }
+        REQUIRE(min_value == parallel::allReduceMin(local_min));
+      }
+
+      const double max_value = max(positive_function);
+      SECTION("max")
+      {
+        double local_max = -std::numeric_limits<double>::max();
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          local_max = std::max(local_max, positive_function[cell_id]);
+        }
+        REQUIRE(max_value == parallel::allReduceMax(local_max));
+      }
+
+      REQUIRE(min_value < max_value);
+
+      DiscreteFunctionP0 unsigned_function = positive_function - 0.5 * (min_value + max_value);
+
+      SECTION("sqrt")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, sqrt);
+      }
+
+      SECTION("abs")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, abs);
+      }
+
+      SECTION("cos")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, cos);
+      }
+
+      SECTION("sin")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, sin);
+      }
+
+      SECTION("tan")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, tan);
+      }
+
+      DiscreteFunctionP0<Dimension, double> unit_function{mesh};
+
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          unit_function[cell_id] = (2 * (positive_function[cell_id] - min_value) / (max_value - min_value) - 1) * 0.95;
+        });
+
+      SECTION("acos")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, acos);
+      }
+
+      SECTION("asin")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, asin);
+      }
+
+      SECTION("atan")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, atan);
+      }
+
+      SECTION("cosh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, cosh);
+      }
+
+      SECTION("sinh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, sinh);
+      }
+
+      SECTION("tanh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, tanh);
+      }
+
+      SECTION("acosh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, acosh);
+      }
+
+      SECTION("asinh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, asinh);
+      }
+
+      SECTION("atanh")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, atanh);
+      }
+
+      SECTION("exp")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, exp);
+      }
+
+      SECTION("log")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, log);
+      }
+
+      SECTION("max(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(cos(positive_function), sin(positive_function), max);
+      }
+
+      SECTION("max(0.2,vh)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.2, sin(positive_function), max);
+      }
+
+      SECTION("max(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(cos(positive_function), 0.2, max);
+      }
+
+      SECTION("atan2(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(positive_function, 2 + positive_function, atan2);
+      }
+
+      SECTION("atan2(0.5,uh)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.5, 2 + positive_function, atan2);
+      }
+
+      SECTION("atan2(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(2 + cos(positive_function), 0.2, atan2);
+      }
+
+      SECTION("pow(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(positive_function, 0.5 * positive_function, pow);
+      }
+
+      SECTION("pow(uh,0.5)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.5, positive_function, pow);
+      }
+
+      SECTION("pow(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(positive_function, 1.3, pow);
+      }
+
+      SECTION("min(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(sin(positive_function), cos(positive_function), min);
+      }
+
+      SECTION("min(uh,0.5)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.5, cos(positive_function), min);
+      }
+
+      SECTION("min(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(sin(positive_function), 0.5, min);
+      }
+
+      SECTION("max(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(sin(positive_function), cos(positive_function), max);
+      }
+
+      SECTION("min(uh,0.5)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.1, cos(positive_function), max);
+      }
+
+      SECTION("min(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(sin(positive_function), 0.1, max);
+      }
+
+      SECTION("dot(uh,hv)")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        DiscreteFunctionP0<Dimension, TinyVector<2>> vh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            vh[cell_id]    = TinyVector<2>{2.3 * x, 1 - x};
+          });
+
+        CHECK_STD_BINARY_MATH_FUNCTION(uh, vh, dot);
+      }
+
+      SECTION("dot(uh,v)")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        const TinyVector<2> v{1, 2};
+
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(uh, v, dot);
+      }
+
+      SECTION("dot(u,hv)")
+      {
+        const TinyVector<2> u{3, -2};
+
+        DiscreteFunctionP0<Dimension, TinyVector<2>> vh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            vh[cell_id]    = TinyVector<2>{2.3 * x, 1 - x};
+          });
+
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, dot);
+      }
+
+      SECTION("scalar sum")
+      {
+        const CellValue<const double> cell_value = positive_function.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(positive_function));
+      }
+
+      SECTION("vector sum")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+        const CellValue<const TinyVector<2>> cell_value = uh.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(uh));
+      }
+
+      SECTION("matrix sum")
+      {
+        DiscreteFunctionP0<Dimension, TinyMatrix<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3, 2 * x, 3 * x - 1};
+          });
+        const CellValue<const TinyMatrix<2>> cell_value = uh.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(uh));
+      }
+
+      SECTION("integrate scalar")
+      {
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<double> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            cell_value[cell_id] = cell_volume[cell_id] * positive_function[cell_id];
+          });
+
+        REQUIRE(integrate(positive_function) == Catch::Approx(sum(cell_value)));
+      }
+
+      SECTION("integrate vector")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<TinyVector<2>> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(),
+          PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_volume[cell_id] * uh[cell_id]; });
+
+        REQUIRE(integrate(uh)[0] == Catch::Approx(sum(cell_value)[0]));
+        REQUIRE(integrate(uh)[1] == Catch::Approx(sum(cell_value)[1]));
+      }
+
+      SECTION("integrate matrix")
+      {
+        DiscreteFunctionP0<Dimension, TinyMatrix<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3, 2 * x, 1 - x};
+          });
+
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<TinyMatrix<2>> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(),
+          PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_volume[cell_id] * uh[cell_id]; });
+
+        REQUIRE(integrate(uh)(0, 0) == Catch::Approx(sum(cell_value)(0, 0)));
+        REQUIRE(integrate(uh)(0, 1) == Catch::Approx(sum(cell_value)(0, 1)));
+        REQUIRE(integrate(uh)(1, 0) == Catch::Approx(sum(cell_value)(1, 0)));
+        REQUIRE(integrate(uh)(1, 1) == Catch::Approx(sum(cell_value)(1, 1)));
+      }
+    }
+
+    SECTION("3D")
+    {
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();
+
+      constexpr size_t Dimension = 3;
+
+      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      DiscreteFunctionP0<Dimension, double> positive_function{mesh};
+
+      parallel_for(
+        mesh->numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { positive_function[cell_id] = 1 + std::abs(xj[cell_id][0]); });
+
+      const double min_value = min(positive_function);
+      SECTION("min")
+      {
+        double local_min = std::numeric_limits<double>::max();
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          local_min = std::min(local_min, positive_function[cell_id]);
+        }
+        REQUIRE(min_value == parallel::allReduceMin(local_min));
+      }
+
+      const double max_value = max(positive_function);
+      SECTION("max")
+      {
+        double local_max = -std::numeric_limits<double>::max();
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          local_max = std::max(local_max, positive_function[cell_id]);
+        }
+        REQUIRE(max_value == parallel::allReduceMax(local_max));
+      }
+
+      REQUIRE(min_value < max_value);
+
+      DiscreteFunctionP0 unsigned_function = positive_function - 0.5 * (min_value + max_value);
+
+      SECTION("sqrt")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, sqrt);
+      }
+
+      SECTION("abs")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, abs);
+      }
+
+      SECTION("cos")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, cos);
+      }
+
+      SECTION("sin")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, sin);
+      }
+
+      SECTION("tan")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, tan);
+      }
+
+      DiscreteFunctionP0<Dimension, double> unit_function{mesh};
+
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          unit_function[cell_id] = (2 * (positive_function[cell_id] - min_value) / (max_value - min_value) - 1) * 0.95;
+        });
+
+      SECTION("acos")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, acos);
+      }
+
+      SECTION("asin")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, asin);
+      }
+
+      SECTION("atan")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, atan);
+      }
+
+      SECTION("cosh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, cosh);
+      }
+
+      SECTION("sinh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, sinh);
+      }
+
+      SECTION("tanh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, tanh);
+      }
+
+      SECTION("acosh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, acosh);
+      }
+
+      SECTION("asinh")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, asinh);
+      }
+
+      SECTION("atanh")
+      {
+        CHECK_STD_MATH_FUNCTION(unit_function, atanh);
+      }
+
+      SECTION("exp")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, exp);
+      }
+
+      SECTION("log")
+      {
+        CHECK_STD_MATH_FUNCTION(positive_function, log);
+      }
+
+      SECTION("max(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(cos(positive_function), sin(positive_function), max);
+      }
+
+      SECTION("max(0.2,vh)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.2, sin(positive_function), max);
+      }
+
+      SECTION("max(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(cos(positive_function), 0.2, max);
+      }
+
+      SECTION("atan2(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(positive_function, 2 + positive_function, atan2);
+      }
+
+      SECTION("atan2(0.5,uh)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.5, 2 + positive_function, atan2);
+      }
+
+      SECTION("atan2(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(2 + cos(positive_function), 0.2, atan2);
+      }
+
+      SECTION("pow(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(positive_function, 0.5 * positive_function, pow);
+      }
+
+      SECTION("pow(uh,0.5)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.5, positive_function, pow);
+      }
+
+      SECTION("pow(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(positive_function, 1.3, pow);
+      }
+
+      SECTION("min(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(sin(positive_function), cos(positive_function), min);
+      }
+
+      SECTION("min(uh,0.5)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.5, cos(positive_function), min);
+      }
+
+      SECTION("min(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(sin(positive_function), 0.5, min);
+      }
+
+      SECTION("max(uh,hv)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION(sin(positive_function), cos(positive_function), max);
+      }
+
+      SECTION("min(uh,0.5)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(0.1, cos(positive_function), max);
+      }
+
+      SECTION("min(uh,0.2)")
+      {
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(sin(positive_function), 0.1, max);
+      }
+
+      SECTION("dot(uh,hv)")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        DiscreteFunctionP0<Dimension, TinyVector<2>> vh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            vh[cell_id]    = TinyVector<2>{2.3 * x, 1 - x};
+          });
+
+        CHECK_STD_BINARY_MATH_FUNCTION(uh, vh, dot);
+      }
+
+      SECTION("dot(uh,v)")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        const TinyVector<2> v{1, 2};
+
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(uh, v, dot);
+      }
+
+      SECTION("dot(u,hv)")
+      {
+        const TinyVector<2> u{3, -2};
+
+        DiscreteFunctionP0<Dimension, TinyVector<2>> vh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            vh[cell_id]    = TinyVector<2>{2.3 * x, 1 - x};
+          });
+
+        CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, dot);
+      }
+
+      SECTION("scalar sum")
+      {
+        const CellValue<const double> cell_value = positive_function.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(positive_function));
+      }
+
+      SECTION("vector sum")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+        const CellValue<const TinyVector<2>> cell_value = uh.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(uh));
+      }
+
+      SECTION("matrix sum")
+      {
+        DiscreteFunctionP0<Dimension, TinyMatrix<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3, 2 * x, 3 * x - 1};
+          });
+        const CellValue<const TinyMatrix<2>> cell_value = uh.cellValues();
+
+        REQUIRE(sum(cell_value) == sum(uh));
+      }
+
+      SECTION("integrate scalar")
+      {
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<double> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            cell_value[cell_id] = cell_volume[cell_id] * positive_function[cell_id];
+          });
+
+        REQUIRE(integrate(positive_function) == Catch::Approx(sum(cell_value)));
+      }
+
+      SECTION("integrate vector")
+      {
+        DiscreteFunctionP0<Dimension, TinyVector<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyVector<2>{x + 1, 2 * x - 3};
+          });
+
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<TinyVector<2>> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(),
+          PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_volume[cell_id] * uh[cell_id]; });
+
+        REQUIRE(integrate(uh)[0] == Catch::Approx(sum(cell_value)[0]));
+        REQUIRE(integrate(uh)[1] == Catch::Approx(sum(cell_value)[1]));
+      }
+
+      SECTION("integrate matrix")
+      {
+        DiscreteFunctionP0<Dimension, TinyMatrix<2>> uh{mesh};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            const double x = xj[cell_id][0];
+            uh[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3, 2 * x, 1 - x};
+          });
+
+        const CellValue<const double> cell_volume = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+        CellValue<TinyMatrix<2>> cell_value{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(),
+          PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_volume[cell_id] * uh[cell_id]; });
+
+        REQUIRE(integrate(uh)(0, 0) == Catch::Approx(sum(cell_value)(0, 0)));
+        REQUIRE(integrate(uh)(0, 1) == Catch::Approx(sum(cell_value)(0, 1)));
+        REQUIRE(integrate(uh)(1, 0) == Catch::Approx(sum(cell_value)(1, 0)));
+        REQUIRE(integrate(uh)(1, 1) == Catch::Approx(sum(cell_value)(1, 1)));
+      }
     }
   }