diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index d3fa06d8e6a45c0b6f7453ab091b1ce9057d0a20..52e69c29af17a994d1a51d8d87de48c02249645c 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -7,7 +7,7 @@
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
-#define DISCRETE_FUNCTION_CALL(FUNCTION, ARG)                                                                         \
+#define DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(FUNCTION, ARG)                                                           \
   if (ARG->dataType() == ASTNodeDataType::double_t and ARG->descriptor().type() == DiscreteFunctionType::P0) {        \
     switch (ARG->mesh()->dimension()) {                                                                               \
     case 1: {                                                                                                         \
@@ -27,55 +27,55 @@
     }                                                                                                                 \
     }                                                                                                                 \
   } else {                                                                                                            \
-    throw NormalError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(ARG));             \
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(ARG));                                       \
   }
 
 std::shared_ptr<const IDiscreteFunction>
 sqrt(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(sqrt, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(sqrt, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 abs(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(abs, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(abs, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 sin(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(sin, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(sin, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 cos(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(cos, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(cos, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 tan(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(tan, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(tan, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 asin(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(asin, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(asin, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 acos(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(acos, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(acos, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 atan(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(atan, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(atan, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
@@ -105,15 +105,14 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<c
       return std::make_shared<const DiscreteFunctionType>(
         atan2(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -134,15 +133,14 @@ atan2(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(atan2(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
 }
 
@@ -163,64 +161,63 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(atan2(dynamic_cast<const DiscreteFunctionType&>(*f), a));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
 }
 
 std::shared_ptr<const IDiscreteFunction>
 sinh(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(sinh, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(sinh, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 cosh(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(cosh, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(cosh, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 tanh(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(tanh, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(tanh, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 asinh(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(asinh, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(asinh, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 acosh(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(acosh, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(acosh, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 atanh(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(atanh, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(atanh, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 exp(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(exp, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(exp, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
 log(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  DISCRETE_FUNCTION_CALL(log, f);
+  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(log, f);
 }
 
 std::shared_ptr<const IDiscreteFunction>
@@ -250,15 +247,14 @@ pow(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<con
       return std::make_shared<const DiscreteFunctionType>(
         pow(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -279,15 +275,14 @@ pow(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(pow(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
 }
 
@@ -308,15 +303,14 @@ pow(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(pow(dynamic_cast<const DiscreteFunctionType&>(*f), a));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
 }
 
@@ -349,9 +343,11 @@ dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<con
     return std::make_shared<const DiscreteFunctionResultType>(
       dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
   }
+    // LCOV_EXCL_START
   default: {
-    throw UnexpectedError("invalid data dimension " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -377,15 +373,14 @@ dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<con
     case 3: {
       return dot<3>(f, g);
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -418,15 +413,14 @@ dot(const std::shared_ptr<const IDiscreteFunction>& f, const TinyVector<VectorDi
     case 3: {
       return dot<3, VectorDimension>(f, a);
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
 }
 
@@ -459,15 +453,14 @@ dot(const TinyVector<VectorDimension>& a, const std::shared_ptr<const IDiscreteF
     case 3: {
       return dot<3, VectorDimension>(a, f);
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
 }
 
@@ -506,12 +499,14 @@ min(const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return min(dynamic_cast<const DiscreteFunctionType&>(*f));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    throw NormalError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
   }
 }
 
@@ -542,15 +537,14 @@ min(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<con
       return std::make_shared<const DiscreteFunctionType>(
         min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -571,15 +565,14 @@ min(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
 }
 
@@ -600,15 +593,14 @@ min(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
 }
 
@@ -629,12 +621,14 @@ max(const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return max(dynamic_cast<const DiscreteFunctionType&>(*f));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    throw NormalError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
   }
 }
 
@@ -665,15 +659,14 @@ max(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<con
       return std::make_shared<const DiscreteFunctionType>(
         max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -694,15 +687,14 @@ max(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
 }
 
@@ -723,15 +715,14 @@ max(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
       using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
       return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    std::stringstream os;
-    os << "incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-       << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(a);
-    throw NormalError(os.str());
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
 }
 
@@ -753,12 +744,14 @@ sum_of(const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, ValueT>;
       return sum(dynamic_cast<const DiscreteFunctionType&>(*f));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_STOP
     }
   } else {
-    throw NormalError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
   }
 }
 
@@ -794,12 +787,14 @@ integral_of(const std::shared_ptr<const IDiscreteFunction>& f)
       using DiscreteFunctionType = DiscreteFunctionP0<3, ValueT>;
       return integrate(dynamic_cast<const DiscreteFunctionType&>(*f));
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
+      // LCOV_EXCL_START
     }
   } else {
-    throw NormalError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
   }
 }
 
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp b/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
index 42150443939e990e93c17996a62e2b7630727f55..1e61f917dfa448add18b562dd403a70101dbec3b 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
@@ -5,6 +5,8 @@
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
+#include <string>
+
 struct EmbeddedIDiscreteFunctionUtils
 {
   template <typename T>
@@ -29,6 +31,20 @@ struct EmbeddedIDiscreteFunctionUtils
   {
     return isSameDiscretization(*f, *g);
   }
+
+  template <typename T1, typename T2>
+  PUGS_INLINE static std::string
+  incompatibleOperandTypes(const T1& t1, const T2& t2)
+  {
+    return "incompatible operand types " + getOperandTypeName(t1) + " and " + getOperandTypeName(t2);
+  }
+
+  template <typename T>
+  PUGS_INLINE static std::string
+  invalidOperandType(const T& t)
+  {
+    return "invalid operand type " + getOperandTypeName(t);
+  }
 };
 
 #endif   // EMBEDDED_I_DISCRETE_FUNCTION_UTILS_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 379526adfb7abe169b24408c16c24c27c4c4f2f3..63a95ae67d58d3028c9af1996955588f96992346 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -112,6 +112,7 @@ add_executable (mpi_unit_tests
   test_DiscreteFunctionP0.cpp
   test_DiscreteFunctionP0Vector.cpp
   test_DiscreteFunctionVectorInterpoler.cpp
+  test_EmbeddedIDiscreteFunctionMathFunctions.cpp
   test_InterpolateItemArray.cpp
   test_InterpolateItemValue.cpp
   test_ItemArray.cpp
diff --git a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7cf6866e96a6a9141be5771bd93768a1a6d1d415
--- /dev/null
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -0,0 +1,1421 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+#include <language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+#define CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(P_U, FCT)          \
+  {                                                                         \
+    using DiscreteFunctionType = const std::decay_t<decltype(*P_U)>;        \
+    std::shared_ptr p_fu       = ::FCT(P_U);                                \
+                                                                            \
+    REQUIRE(p_fu.use_count() > 0);                                          \
+    REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionType&>(*p_fu));      \
+                                                                            \
+    const auto& fu = dynamic_cast<const DiscreteFunctionType&>(*p_fu);      \
+                                                                            \
+    auto values  = P_U->cellValues();                                       \
+    bool is_same = true;                                                    \
+    for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) { \
+      if (fu[cell_id] != std::FCT(values[cell_id])) {                       \
+        is_same = false;                                                    \
+        break;                                                              \
+      }                                                                     \
+    }                                                                       \
+                                                                            \
+    REQUIRE(is_same);                                                       \
+  }
+
+#define CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(P_LHS, P_RHS, FCT)             \
+  {                                                                                 \
+    using DiscreteFunctionType = const std::decay_t<decltype(FCT(*P_LHS, *P_RHS))>; \
+    std::shared_ptr p_fuv      = ::FCT(P_LHS, P_RHS);                               \
+                                                                                    \
+    REQUIRE(p_fuv.use_count() > 0);                                                 \
+    REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionType&>(*p_fuv));             \
+                                                                                    \
+    const auto& fuv = dynamic_cast<const DiscreteFunctionType&>(*p_fuv);            \
+                                                                                    \
+    auto lhs_values = P_LHS->cellValues();                                          \
+    auto rhs_values = P_RHS->cellValues();                                          \
+    bool is_same    = true;                                                         \
+    for (CellId cell_id = 0; cell_id < lhs_values.numberOfItems(); ++cell_id) {     \
+      using namespace std;                                                          \
+      if (fuv[cell_id] != FCT(lhs_values[cell_id], rhs_values[cell_id])) {          \
+        is_same = false;                                                            \
+        break;                                                                      \
+      }                                                                             \
+    }                                                                               \
+                                                                                    \
+    REQUIRE(is_same);                                                               \
+  }
+
+#define CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(P_LHS, RHS, FCT)           \
+  {                                                                              \
+    using DiscreteFunctionType = const std::decay_t<decltype(FCT(*P_LHS, RHS))>; \
+    std::shared_ptr p_fuv      = ::FCT(P_LHS, RHS);                              \
+                                                                                 \
+    REQUIRE(p_fuv.use_count() > 0);                                              \
+    REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionType&>(*p_fuv));          \
+                                                                                 \
+    const auto& fuv = dynamic_cast<const DiscreteFunctionType&>(*p_fuv);         \
+                                                                                 \
+    auto lhs_values = P_LHS->cellValues();                                       \
+    bool is_same    = true;                                                      \
+    for (CellId cell_id = 0; cell_id < lhs_values.numberOfItems(); ++cell_id) {  \
+      using namespace std;                                                       \
+      if (fuv[cell_id] != FCT(lhs_values[cell_id], RHS)) {                       \
+        is_same = false;                                                         \
+        break;                                                                   \
+      }                                                                          \
+    }                                                                            \
+                                                                                 \
+    REQUIRE(is_same);                                                            \
+  }
+
+#define CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(LHS, P_RHS, FCT)           \
+  {                                                                              \
+    using DiscreteFunctionType = const std::decay_t<decltype(FCT(LHS, *P_RHS))>; \
+    std::shared_ptr p_fuv      = ::FCT(LHS, P_RHS);                              \
+                                                                                 \
+    REQUIRE(p_fuv.use_count() > 0);                                              \
+    REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionType&>(*p_fuv));          \
+                                                                                 \
+    const auto& fuv = dynamic_cast<const DiscreteFunctionType&>(*p_fuv);         \
+                                                                                 \
+    auto rhs_values = P_RHS->cellValues();                                       \
+    bool is_same    = true;                                                      \
+    for (CellId cell_id = 0; cell_id < rhs_values.numberOfItems(); ++cell_id) {  \
+      using namespace std;                                                       \
+      if (fuv[cell_id] != FCT(LHS, rhs_values[cell_id])) {                       \
+        is_same = false;                                                         \
+        break;                                                                   \
+      }                                                                          \
+    }                                                                            \
+                                                                                 \
+    REQUIRE(is_same);                                                            \
+  }
+
+TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
+{
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    using Rd = TinyVector<Dimension>;
+
+    std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();
+
+    std::shared_ptr other_mesh =
+      std::make_shared<Mesh<Connectivity<Dimension>>>(mesh->shared_connectivity(), mesh->xr());
+
+    CellValue<const Rd> xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+    CellValue<double> values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.2 + std::cos(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    CellValue<double> positive_values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 2 + std::sin(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    CellValue<double> bounded_values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.9 * std::sin(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    std::shared_ptr p_u            = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, values);
+    std::shared_ptr p_other_mesh_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(other_mesh, values);
+    std::shared_ptr p_positive_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, positive_values);
+    std::shared_ptr p_bounded_u  = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, bounded_values);
+
+    std::shared_ptr p_R1_u = [=] {
+      CellValue<TinyVector<1>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id][0] = 2 * xj[cell_id][0] + 1; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R1_v = [=] {
+      CellValue<TinyVector<1>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { vj[cell_id][0] = xj[cell_id][0] * xj[cell_id][0] + 1; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R1_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(other_mesh, p_R1_u->cellValues());
+
+    constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
+      if constexpr (Dimension == 1) {
+        return {x[0], 1 + x[0] * x[0]};
+      } else if constexpr (Dimension == 2) {
+        return {x[0], x[1]};
+      } else if constexpr (Dimension == 3) {
+        return {x[0], x[1] + x[2]};
+      }
+    };
+
+    std::shared_ptr p_R2_u = [=] {
+      CellValue<TinyVector<2>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+          uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R2_v = [=] {
+      CellValue<TinyVector<2>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+          vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R2_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(other_mesh, p_R2_u->cellValues());
+
+    constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
+      if constexpr (Dimension == 1) {
+        return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+      } else if constexpr (Dimension == 2) {
+        return {x[0], x[1], x[0] + x[1]};
+      } else if constexpr (Dimension == 3) {
+        return {x[0], x[1], x[2]};
+      }
+    };
+
+    std::shared_ptr p_R3_u = [=] {
+      CellValue<TinyVector<3>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+          uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R3_v = [=] {
+      CellValue<TinyVector<3>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+          vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R3_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(other_mesh, p_R3_u->cellValues());
+
+    std::shared_ptr p_R1x1_u = [=] {
+      CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R2x2_u = [=] {
+      CellValue<TinyMatrix<2>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+
+          uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
+                         2 * x[1], -x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R3x3_u = [=] {
+      CellValue<TinyMatrix<3>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+
+          uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
+                         2 * x[1],        -x[0],           x[0] - x[1],   //
+                         3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
+    }();
+
+    SECTION("sqrt Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, sqrt);
+      REQUIRE_THROWS_WITH(sqrt(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("abs Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, abs);
+      REQUIRE_THROWS_WITH(abs(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("sin Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, sin);
+      REQUIRE_THROWS_WITH(sin(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("cos Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, cos);
+      REQUIRE_THROWS_WITH(cos(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("tan Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, tan);
+      REQUIRE_THROWS_WITH(tan(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("asin Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, asin);
+      REQUIRE_THROWS_WITH(asin(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("acos Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, acos);
+      REQUIRE_THROWS_WITH(acos(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atan Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, atan);
+      REQUIRE_THROWS_WITH(atan(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("sinh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, sinh);
+      REQUIRE_THROWS_WITH(sinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("cosh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, cosh);
+      REQUIRE_THROWS_WITH(cosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("tanh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, tanh);
+      REQUIRE_THROWS_WITH(tanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("asinh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, asinh);
+      REQUIRE_THROWS_WITH(asinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("acosh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, acosh);
+      REQUIRE_THROWS_WITH(acosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atanh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, atanh);
+      REQUIRE_THROWS_WITH(atanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("exp Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, exp);
+      REQUIRE_THROWS_WITH(exp(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("log Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, log);
+      REQUIRE_THROWS_WITH(log(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atan2 Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, atan2);
+      REQUIRE_THROWS_WITH(atan2(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(atan2(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+      REQUIRE_THROWS_WITH(atan2(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+    }
+
+    SECTION("atan2 Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 3.6, atan2);
+      REQUIRE_THROWS_WITH(atan2(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("atan2 R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.4, p_u, atan2);
+      REQUIRE_THROWS_WITH(atan2(2.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("min Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, min);
+      REQUIRE_THROWS_WITH(::min(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(::min(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(::min(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("min Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, min);
+      REQUIRE_THROWS_WITH(min(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("min R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, min);
+      REQUIRE_THROWS_WITH(min(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("min Vh -> R")
+    {
+      REQUIRE(min(std::shared_ptr<const IDiscreteFunction>{p_u}) == min(p_u->cellValues()));
+      REQUIRE_THROWS_WITH(min(std::shared_ptr<const IDiscreteFunction>{p_R1_u}),
+                          "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("max Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, max);
+      REQUIRE_THROWS_WITH(::max(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(::max(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(::max(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("max Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, max);
+      REQUIRE_THROWS_WITH(max(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("max Vh -> R")
+    {
+      REQUIRE(max(std::shared_ptr<const IDiscreteFunction>{p_u}) == max(p_u->cellValues()));
+      REQUIRE_THROWS_WITH(max(std::shared_ptr<const IDiscreteFunction>{p_R1_u}),
+                          "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("max R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, max);
+      REQUIRE_THROWS_WITH(max(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("pow Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, pow);
+      REQUIRE_THROWS_WITH(pow(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(pow(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(pow(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("pow Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_positive_u, 3.3, pow);
+      REQUIRE_THROWS_WITH(pow(p_R1_u, 3.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("pow R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.1, p_u, pow);
+      REQUIRE_THROWS_WITH(pow(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("dot Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, dot);
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, dot);
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, dot);
+      REQUIRE_THROWS_WITH(dot(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+    }
+
+    SECTION("dot Vh*Rd -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), dot);
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), dot);
+      REQUIRE_THROWS_WITH(dot(p_R1_u, (TinyVector<2>{-6, 2})), "error: incompatible operand types Vh(P0:R^1) and R^2");
+      REQUIRE_THROWS_WITH(dot(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+                          "error: incompatible operand types Vh(P0:R^2) and R^3");
+      REQUIRE_THROWS_WITH(dot(p_R3_u, (TinyVector<1>{-1})), "error: incompatible operand types Vh(P0:R^3) and R^1");
+    }
+
+    SECTION("dot Rd*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, dot);
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, dot);
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, dot);
+      REQUIRE_THROWS_WITH(dot((TinyVector<2>{-6, 2}), p_R1_u), "error: incompatible operand types R^2 and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(dot((TinyVector<3>{-1, 5, 2}), p_R2_u),
+                          "error: incompatible operand types R^3 and Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
+    }
+
+    SECTION("sum_of_R* Vh -> R*")
+    {
+      REQUIRE(sum_of<double>(p_u) == sum(p_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<1>>(p_R1_u) == sum(p_R1_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<2>>(p_R2_u) == sum(p_R2_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<3>>(p_R3_u) == sum(p_R3_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<1>>(p_R1x1_u) == sum(p_R1x1_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<2>>(p_R2x2_u) == sum(p_R2x2_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<3>>(p_R3x3_u) == sum(p_R3x3_u->cellValues()));
+
+      REQUIRE_THROWS_WITH(sum_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+    }
+
+    SECTION("integral_of_R* Vh -> R*")
+    {
+      auto integrate_locally = [&](const auto& cell_values) {
+        const auto& Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+        using DataType = decltype(double{} * cell_values[CellId{0}]);
+        CellValue<DataType> local_integral{mesh->connectivity()};
+        parallel_for(
+          local_integral.numberOfItems(),
+          PUGS_LAMBDA(const CellId cell_id) { local_integral[cell_id] = Vj[cell_id] * cell_values[cell_id]; });
+        return local_integral;
+      };
+
+      REQUIRE(integral_of<double>(p_u) == sum(integrate_locally(p_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<1>>(p_R1_u) == sum(integrate_locally(p_R1_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<2>>(p_R2_u) == sum(integrate_locally(p_R2_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<3>>(p_R3_u) == sum(integrate_locally(p_R3_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<1>>(p_R1x1_u) == sum(integrate_locally(p_R1x1_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<2>>(p_R2x2_u) == sum(integrate_locally(p_R2x2_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<3>>(p_R3x3_u) == sum(integrate_locally(p_R3x3_u->cellValues())));
+
+      REQUIRE_THROWS_WITH(integral_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    using Rd = TinyVector<Dimension>;
+
+    std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();
+
+    std::shared_ptr other_mesh =
+      std::make_shared<Mesh<Connectivity<Dimension>>>(mesh->shared_connectivity(), mesh->xr());
+
+    CellValue<const Rd> xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+    CellValue<double> values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.2 + std::cos(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    CellValue<double> positive_values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 2 + std::sin(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    CellValue<double> bounded_values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.9 * std::sin(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    std::shared_ptr p_u            = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, values);
+    std::shared_ptr p_other_mesh_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(other_mesh, values);
+    std::shared_ptr p_positive_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, positive_values);
+    std::shared_ptr p_bounded_u  = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, bounded_values);
+
+    std::shared_ptr p_R1_u = [=] {
+      CellValue<TinyVector<1>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id][0] = 2 * xj[cell_id][0] + 1; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R1_v = [=] {
+      CellValue<TinyVector<1>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { vj[cell_id][0] = xj[cell_id][0] * xj[cell_id][0] + 1; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R1_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(other_mesh, p_R1_u->cellValues());
+
+    constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
+      if constexpr (Dimension == 1) {
+        return {x[0], 1 + x[0] * x[0]};
+      } else if constexpr (Dimension == 2) {
+        return {x[0], x[1]};
+      } else if constexpr (Dimension == 3) {
+        return {x[0], x[1] + x[2]};
+      }
+    };
+
+    std::shared_ptr p_R2_u = [=] {
+      CellValue<TinyVector<2>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+          uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R2_v = [=] {
+      CellValue<TinyVector<2>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+          vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R2_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(other_mesh, p_R2_u->cellValues());
+
+    constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
+      if constexpr (Dimension == 1) {
+        return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+      } else if constexpr (Dimension == 2) {
+        return {x[0], x[1], x[0] + x[1]};
+      } else if constexpr (Dimension == 3) {
+        return {x[0], x[1], x[2]};
+      }
+    };
+
+    std::shared_ptr p_R3_u = [=] {
+      CellValue<TinyVector<3>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+          uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R3_v = [=] {
+      CellValue<TinyVector<3>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+          vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R3_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(other_mesh, p_R3_u->cellValues());
+
+    std::shared_ptr p_R1x1_u = [=] {
+      CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R2x2_u = [=] {
+      CellValue<TinyMatrix<2>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+
+          uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
+                         2 * x[1], -x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R3x3_u = [=] {
+      CellValue<TinyMatrix<3>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+
+          uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
+                         2 * x[1],        -x[0],           x[0] - x[1],   //
+                         3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
+    }();
+
+    SECTION("sqrt Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, sqrt);
+      REQUIRE_THROWS_WITH(sqrt(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("abs Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, abs);
+      REQUIRE_THROWS_WITH(abs(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("sin Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, sin);
+      REQUIRE_THROWS_WITH(sin(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("cos Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, cos);
+      REQUIRE_THROWS_WITH(cos(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("tan Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, tan);
+      REQUIRE_THROWS_WITH(tan(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("asin Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, asin);
+      REQUIRE_THROWS_WITH(asin(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("acos Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, acos);
+      REQUIRE_THROWS_WITH(acos(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atan Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, atan);
+      REQUIRE_THROWS_WITH(atan(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("sinh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, sinh);
+      REQUIRE_THROWS_WITH(sinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("cosh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, cosh);
+      REQUIRE_THROWS_WITH(cosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("tanh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, tanh);
+      REQUIRE_THROWS_WITH(tanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("asinh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, asinh);
+      REQUIRE_THROWS_WITH(asinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("acosh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, acosh);
+      REQUIRE_THROWS_WITH(acosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atanh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, atanh);
+      REQUIRE_THROWS_WITH(atanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("exp Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, exp);
+      REQUIRE_THROWS_WITH(exp(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("log Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, log);
+      REQUIRE_THROWS_WITH(log(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atan2 Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, atan2);
+      REQUIRE_THROWS_WITH(atan2(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(atan2(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+      REQUIRE_THROWS_WITH(atan2(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+    }
+
+    SECTION("atan2 Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 3.6, atan2);
+      REQUIRE_THROWS_WITH(atan2(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("atan2 R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.4, p_u, atan2);
+      REQUIRE_THROWS_WITH(atan2(2.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("min Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, min);
+      REQUIRE_THROWS_WITH(::min(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(::min(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(::min(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("min Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, min);
+      REQUIRE_THROWS_WITH(min(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("min R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, min);
+      REQUIRE_THROWS_WITH(min(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("min Vh -> R")
+    {
+      REQUIRE(min(std::shared_ptr<const IDiscreteFunction>{p_u}) == min(p_u->cellValues()));
+      REQUIRE_THROWS_WITH(min(std::shared_ptr<const IDiscreteFunction>{p_R1_u}),
+                          "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("max Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, max);
+      REQUIRE_THROWS_WITH(::max(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(::max(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(::max(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("max Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, max);
+      REQUIRE_THROWS_WITH(max(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("max Vh -> R")
+    {
+      REQUIRE(max(std::shared_ptr<const IDiscreteFunction>{p_u}) == max(p_u->cellValues()));
+      REQUIRE_THROWS_WITH(max(std::shared_ptr<const IDiscreteFunction>{p_R1_u}),
+                          "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("max R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, max);
+      REQUIRE_THROWS_WITH(max(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("pow Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, pow);
+      REQUIRE_THROWS_WITH(pow(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(pow(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(pow(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("pow Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_positive_u, 3.3, pow);
+      REQUIRE_THROWS_WITH(pow(p_R1_u, 3.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("pow R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.1, p_u, pow);
+      REQUIRE_THROWS_WITH(pow(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("dot Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, dot);
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, dot);
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, dot);
+      REQUIRE_THROWS_WITH(dot(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+    }
+
+    SECTION("dot Vh*Rd -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), dot);
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), dot);
+      REQUIRE_THROWS_WITH(dot(p_R1_u, (TinyVector<2>{-6, 2})), "error: incompatible operand types Vh(P0:R^1) and R^2");
+      REQUIRE_THROWS_WITH(dot(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+                          "error: incompatible operand types Vh(P0:R^2) and R^3");
+      REQUIRE_THROWS_WITH(dot(p_R3_u, (TinyVector<1>{-1})), "error: incompatible operand types Vh(P0:R^3) and R^1");
+    }
+
+    SECTION("dot Rd*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, dot);
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, dot);
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, dot);
+      REQUIRE_THROWS_WITH(dot((TinyVector<2>{-6, 2}), p_R1_u), "error: incompatible operand types R^2 and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(dot((TinyVector<3>{-1, 5, 2}), p_R2_u),
+                          "error: incompatible operand types R^3 and Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
+    }
+
+    SECTION("sum_of_R* Vh -> R*")
+    {
+      REQUIRE(sum_of<double>(p_u) == sum(p_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<1>>(p_R1_u) == sum(p_R1_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<2>>(p_R2_u) == sum(p_R2_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<3>>(p_R3_u) == sum(p_R3_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<1>>(p_R1x1_u) == sum(p_R1x1_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<2>>(p_R2x2_u) == sum(p_R2x2_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<3>>(p_R3x3_u) == sum(p_R3x3_u->cellValues()));
+
+      REQUIRE_THROWS_WITH(sum_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+    }
+
+    SECTION("integral_of_R* Vh -> R*")
+    {
+      auto integrate_locally = [&](const auto& cell_values) {
+        const auto& Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+        using DataType = decltype(double{} * cell_values[CellId{0}]);
+        CellValue<DataType> local_integral{mesh->connectivity()};
+        parallel_for(
+          local_integral.numberOfItems(),
+          PUGS_LAMBDA(const CellId cell_id) { local_integral[cell_id] = Vj[cell_id] * cell_values[cell_id]; });
+        return local_integral;
+      };
+
+      REQUIRE(integral_of<double>(p_u) == sum(integrate_locally(p_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<1>>(p_R1_u) == sum(integrate_locally(p_R1_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<2>>(p_R2_u) == sum(integrate_locally(p_R2_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<3>>(p_R3_u) == sum(integrate_locally(p_R3_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<1>>(p_R1x1_u) == sum(integrate_locally(p_R1x1_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<2>>(p_R2x2_u) == sum(integrate_locally(p_R2x2_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<3>>(p_R3x3_u) == sum(integrate_locally(p_R3x3_u->cellValues())));
+
+      REQUIRE_THROWS_WITH(integral_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    using Rd = TinyVector<Dimension>;
+
+    std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();
+
+    std::shared_ptr other_mesh =
+      std::make_shared<Mesh<Connectivity<Dimension>>>(mesh->shared_connectivity(), mesh->xr());
+
+    CellValue<const Rd> xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+    CellValue<double> values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.2 + std::cos(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    CellValue<double> positive_values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 2 + std::sin(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    CellValue<double> bounded_values = [=] {
+      CellValue<double> build_values{mesh->connectivity()};
+      parallel_for(
+        build_values.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.9 * std::sin(l2Norm(xj[cell_id])); });
+      return build_values;
+    }();
+
+    std::shared_ptr p_u            = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, values);
+    std::shared_ptr p_other_mesh_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(other_mesh, values);
+    std::shared_ptr p_positive_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, positive_values);
+    std::shared_ptr p_bounded_u  = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, bounded_values);
+
+    std::shared_ptr p_R1_u = [=] {
+      CellValue<TinyVector<1>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id][0] = 2 * xj[cell_id][0] + 1; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R1_v = [=] {
+      CellValue<TinyVector<1>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(),
+        PUGS_LAMBDA(const CellId cell_id) { vj[cell_id][0] = xj[cell_id][0] * xj[cell_id][0] + 1; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R1_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(other_mesh, p_R1_u->cellValues());
+
+    constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
+      if constexpr (Dimension == 1) {
+        return {x[0], 1 + x[0] * x[0]};
+      } else if constexpr (Dimension == 2) {
+        return {x[0], x[1]};
+      } else if constexpr (Dimension == 3) {
+        return {x[0], x[1] + x[2]};
+      }
+    };
+
+    std::shared_ptr p_R2_u = [=] {
+      CellValue<TinyVector<2>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+          uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R2_v = [=] {
+      CellValue<TinyVector<2>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+          vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R2_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(other_mesh, p_R2_u->cellValues());
+
+    constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
+      if constexpr (Dimension == 1) {
+        return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+      } else if constexpr (Dimension == 2) {
+        return {x[0], x[1], x[0] + x[1]};
+      } else if constexpr (Dimension == 3) {
+        return {x[0], x[1], x[2]};
+      }
+    };
+
+    std::shared_ptr p_R3_u = [=] {
+      CellValue<TinyVector<3>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+          uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R3_v = [=] {
+      CellValue<TinyVector<3>> vj{mesh->connectivity()};
+      parallel_for(
+        vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+          vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
+    }();
+
+    std::shared_ptr p_other_mesh_R3_u =
+      std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(other_mesh, p_R3_u->cellValues());
+
+    std::shared_ptr p_R1x1_u = [=] {
+      CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R2x2_u = [=] {
+      CellValue<TinyMatrix<2>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<2> x = to_2d(xj[cell_id]);
+
+          uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
+                         2 * x[1], -x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
+    }();
+
+    std::shared_ptr p_R3x3_u = [=] {
+      CellValue<TinyMatrix<3>> uj{mesh->connectivity()};
+      parallel_for(
+        uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          const TinyVector<3> x = to_3d(xj[cell_id]);
+
+          uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
+                         2 * x[1],        -x[0],           x[0] - x[1],   //
+                         3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+        });
+
+      return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
+    }();
+
+    SECTION("sqrt Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, sqrt);
+      REQUIRE_THROWS_WITH(sqrt(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("abs Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, abs);
+      REQUIRE_THROWS_WITH(abs(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("sin Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, sin);
+      REQUIRE_THROWS_WITH(sin(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("cos Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, cos);
+      REQUIRE_THROWS_WITH(cos(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("tan Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, tan);
+      REQUIRE_THROWS_WITH(tan(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("asin Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, asin);
+      REQUIRE_THROWS_WITH(asin(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("acos Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, acos);
+      REQUIRE_THROWS_WITH(acos(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atan Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, atan);
+      REQUIRE_THROWS_WITH(atan(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("sinh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, sinh);
+      REQUIRE_THROWS_WITH(sinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("cosh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, cosh);
+      REQUIRE_THROWS_WITH(cosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("tanh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, tanh);
+      REQUIRE_THROWS_WITH(tanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("asinh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, asinh);
+      REQUIRE_THROWS_WITH(asinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("acosh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, acosh);
+      REQUIRE_THROWS_WITH(acosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atanh Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_bounded_u, atanh);
+      REQUIRE_THROWS_WITH(atanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("exp Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_u, exp);
+      REQUIRE_THROWS_WITH(exp(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("log Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, log);
+      REQUIRE_THROWS_WITH(log(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("atan2 Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, atan2);
+      REQUIRE_THROWS_WITH(atan2(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(atan2(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+      REQUIRE_THROWS_WITH(atan2(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+    }
+
+    SECTION("atan2 Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 3.6, atan2);
+      REQUIRE_THROWS_WITH(atan2(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("atan2 R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.4, p_u, atan2);
+      REQUIRE_THROWS_WITH(atan2(2.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("min Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, min);
+      REQUIRE_THROWS_WITH(::min(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(::min(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(::min(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("min Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, min);
+      REQUIRE_THROWS_WITH(min(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("min R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, min);
+      REQUIRE_THROWS_WITH(min(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("min Vh -> R")
+    {
+      REQUIRE(min(std::shared_ptr<const IDiscreteFunction>{p_u}) == min(p_u->cellValues()));
+      REQUIRE_THROWS_WITH(min(std::shared_ptr<const IDiscreteFunction>{p_R1_u}),
+                          "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("max Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, max);
+      REQUIRE_THROWS_WITH(::max(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(::max(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(::max(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("max Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, max);
+      REQUIRE_THROWS_WITH(max(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("max Vh -> R")
+    {
+      REQUIRE(max(std::shared_ptr<const IDiscreteFunction>{p_u}) == max(p_u->cellValues()));
+      REQUIRE_THROWS_WITH(max(std::shared_ptr<const IDiscreteFunction>{p_R1_u}),
+                          "error: invalid operand type Vh(P0:R^1)");
+    }
+
+    SECTION("max R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, max);
+      REQUIRE_THROWS_WITH(max(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("pow Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, pow);
+      REQUIRE_THROWS_WITH(pow(p_u, p_other_mesh_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(pow(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(pow(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)");
+    }
+
+    SECTION("pow Vh*R -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_positive_u, 3.3, pow);
+      REQUIRE_THROWS_WITH(pow(p_R1_u, 3.1), "error: incompatible operand types Vh(P0:R^1) and R");
+    }
+
+    SECTION("pow R*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.1, p_u, pow);
+      REQUIRE_THROWS_WITH(pow(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)");
+    }
+
+    SECTION("dot Vh*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, dot);
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, dot);
+      CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, dot);
+      REQUIRE_THROWS_WITH(dot(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes");
+      REQUIRE_THROWS_WITH(dot(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+    }
+
+    SECTION("dot Vh*Rd -> Vh")
+    {
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), dot);
+      CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), dot);
+      REQUIRE_THROWS_WITH(dot(p_R1_u, (TinyVector<2>{-6, 2})), "error: incompatible operand types Vh(P0:R^1) and R^2");
+      REQUIRE_THROWS_WITH(dot(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+                          "error: incompatible operand types Vh(P0:R^2) and R^3");
+      REQUIRE_THROWS_WITH(dot(p_R3_u, (TinyVector<1>{-1})), "error: incompatible operand types Vh(P0:R^3) and R^1");
+    }
+
+    SECTION("dot Rd*Vh -> Vh")
+    {
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, dot);
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, dot);
+      CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, dot);
+      REQUIRE_THROWS_WITH(dot((TinyVector<2>{-6, 2}), p_R1_u), "error: incompatible operand types R^2 and Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(dot((TinyVector<3>{-1, 5, 2}), p_R2_u),
+                          "error: incompatible operand types R^3 and Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
+    }
+
+    SECTION("sum_of_R* Vh -> R*")
+    {
+      REQUIRE(sum_of<double>(p_u) == sum(p_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<1>>(p_R1_u) == sum(p_R1_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<2>>(p_R2_u) == sum(p_R2_u->cellValues()));
+      REQUIRE(sum_of<TinyVector<3>>(p_R3_u) == sum(p_R3_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<1>>(p_R1x1_u) == sum(p_R1x1_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<2>>(p_R2x2_u) == sum(p_R2x2_u->cellValues()));
+      REQUIRE(sum_of<TinyMatrix<3>>(p_R3x3_u) == sum(p_R3x3_u->cellValues()));
+
+      REQUIRE_THROWS_WITH(sum_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+      REQUIRE_THROWS_WITH(sum_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+    }
+
+    SECTION("integral_of_R* Vh -> R*")
+    {
+      auto integrate_locally = [&](const auto& cell_values) {
+        const auto& Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+        using DataType = decltype(double{} * cell_values[CellId{0}]);
+        CellValue<DataType> local_integral{mesh->connectivity()};
+        parallel_for(
+          local_integral.numberOfItems(),
+          PUGS_LAMBDA(const CellId cell_id) { local_integral[cell_id] = Vj[cell_id] * cell_values[cell_id]; });
+        return local_integral;
+      };
+
+      REQUIRE(integral_of<double>(p_u) == sum(integrate_locally(p_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<1>>(p_R1_u) == sum(integrate_locally(p_R1_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<2>>(p_R2_u) == sum(integrate_locally(p_R2_u->cellValues())));
+      REQUIRE(integral_of<TinyVector<3>>(p_R3_u) == sum(integrate_locally(p_R3_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<1>>(p_R1x1_u) == sum(integrate_locally(p_R1x1_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<2>>(p_R2x2_u) == sum(integrate_locally(p_R2x2_u->cellValues())));
+      REQUIRE(integral_of<TinyMatrix<3>>(p_R3x3_u) == sum(integrate_locally(p_R3x3_u->cellValues())));
+
+      REQUIRE_THROWS_WITH(integral_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+      REQUIRE_THROWS_WITH(integral_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+    }
+  }
+}