diff --git a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
index d17b53a0dcc29ad8a5cd6fc1f7774b9a9af088ba..656164d6fc012ef2469f846d3e54f32a9c1d4927 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
@@ -9,17 +9,6 @@
 #include <scheme/IDiscreteFunction.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename LHS_T, typename RHS_T>
-PUGS_INLINE std::string
-invalid_operands(const LHS_T& f, const RHS_T& g)
-{
-  std::ostringstream os;
-  os << "undefined binary operator\n";
-  os << "note: incompatible operand types " << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f) << " and "
-     << EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g);
-  return os.str();
-}
-
 // unary operators
 template <typename UnaryOperatorT, typename DiscreteFunctionT>
 std::shared_ptr<const IDiscreteFunction>
@@ -53,9 +42,11 @@ applyUnaryOperation(const std::shared_ptr<const IDiscreteFunction>& f)
         auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>&>(*f);
         return applyUnaryOperation<UnaryOperatorT>(fh);
       }
+        // LCOV_EXCL_START
       default: {
-        throw UnexpectedError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+        throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
       }
+        // LCOV_EXCL_STOP
       }
     }
     case ASTNodeDataType::matrix_t: {
@@ -73,14 +64,18 @@ applyUnaryOperation(const std::shared_ptr<const IDiscreteFunction>& f)
         auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
         return applyUnaryOperation<UnaryOperatorT>(fh);
       }
+        // LCOV_EXCL_START
       default: {
-        throw UnexpectedError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+        throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
       }
+        // LCOV_EXCL_STOP
       }
     }
+      // LCOV_EXCL_START
     default: {
-      throw UnexpectedError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+      throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
     }
+      // LCOV_EXCL_STOP
     }
     break;
   }
@@ -90,15 +85,19 @@ applyUnaryOperation(const std::shared_ptr<const IDiscreteFunction>& f)
       auto fh = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*f);
       return applyUnaryOperation<UnaryOperatorT>(fh);
     }
+      // LCOV_EXCL_START
     default: {
-      throw UnexpectedError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+      throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
     }
+      // LCOV_EXCL_STOP
     }
     break;
   }
+    // LCOV_EXCL_START
   default: {
-    throw UnexpectedError("invalid operand type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -116,9 +115,11 @@ applyUnaryOperation(const std::shared_ptr<const IDiscreteFunction>& f)
   case 3: {
     return applyUnaryOperation<UnaryOperatorT, 3>(f);
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid mesh dimension");
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -138,7 +139,7 @@ innerCompositionLaw(const DiscreteFunctionT& lhs, const DiscreteFunctionT& rhs)
   using data_type = typename DiscreteFunctionT::data_type;
   if constexpr ((std::is_same_v<language::multiply_op, BinOperatorT> and is_tiny_vector_v<data_type>) or
                 (std::is_same_v<language::divide_op, BinOperatorT> and not std::is_arithmetic_v<data_type>)) {
-    throw NormalError(invalid_operands(lhs, rhs));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(lhs, rhs));
   } else {
     return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(lhs, rhs))>(BinOp<BinOperatorT>{}.eval(lhs, rhs));
   }
@@ -172,10 +173,12 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
 
         return innerCompositionLaw<BinOperatorT>(fh, gh);
       } else {
-        throw NormalError(invalid_operands(f, g));
+        throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
       }
     } else {
-      throw UnexpectedError(invalid_operands(f, g));
+      // LCOV_EXCL_START
+      throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
+      // LCOV_EXCL_STOP
     }
   }
   case ASTNodeDataType::vector_t: {
@@ -199,9 +202,11 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
 
       return innerCompositionLaw<BinOperatorT>(fh, gh);
     }
+      // LCOV_EXCL_START
     default: {
-      throw NormalError(invalid_operands(f, g));
+      throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
     }
+      // LCOV_EXCL_STOP
     }
   }
   case ASTNodeDataType::matrix_t: {
@@ -226,14 +231,18 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
 
       return innerCompositionLaw<BinOperatorT>(fh, gh);
     }
+      // LCOV_EXCL_START
     default: {
-      throw UnexpectedError("invalid data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+      throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
     }
+      // LCOV_EXCL_STOP
     }
   }
+    // LCOV_EXCL_START
   default: {
-    throw UnexpectedError("invalid data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -247,9 +256,7 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
     throw NormalError("operands are defined on different meshes");
   }
 
-  if (not EmbeddedIDiscreteFunctionUtils::isSameDiscretization(f, g)) {
-    throw NormalError(invalid_operands(f, g));
-  }
+  Assert(EmbeddedIDiscreteFunctionUtils::isSameDiscretization(f, g));
 
   switch (mesh->dimension()) {
   case 1: {
@@ -261,9 +268,11 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
   case 3: {
     return innerCompositionLaw<BinOperatorT, 3>(f, g);
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid mesh dimension");
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -289,24 +298,18 @@ applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const ID
 
   switch (g->dataType()) {
   case ASTNodeDataType::double_t: {
-    if constexpr (not std::is_same_v<lhs_data_type, double>) {
-      if constexpr (not is_tiny_matrix_v<lhs_data_type>) {
-        auto gh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*g);
-
-        return applyBinaryOperation<BinOperatorT>(fh, gh);
-      } else {
-        throw NormalError(invalid_operands(fh, g));
-      }
-    } else if constexpr (std::is_same_v<BinOperatorT, language::multiply_op> and
-                         std::is_same_v<DiscreteFunctionT, DiscreteFunctionP0<Dimension, double>>) {
+    if constexpr (std::is_same_v<BinOperatorT, language::multiply_op> and
+                  std::is_same_v<DiscreteFunctionT, DiscreteFunctionP0<Dimension, double>>) {
       if (g->descriptor().type() == DiscreteFunctionType::P0Vector) {
         auto gh = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*g);
         return applyBinaryOperation<BinOperatorT>(fh, gh);
       } else {
-        throw NormalError(invalid_operands(fh, g));
+        // LCOV_EXCL_START
+        throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(fh, g));
+        // LCOV_EXCL_STOP
       }
     } else {
-      throw UnexpectedError("should have called innerCompositionLaw");
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(fh, g));
     }
   }
   case ASTNodeDataType::vector_t: {
@@ -319,7 +322,7 @@ applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const ID
 
           return applyBinaryOperation<BinOperatorT>(fh, gh);
         } else {
-          throw NormalError(invalid_operands(fh, g));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(fh, g));
         }
       }
       case 2: {
@@ -329,7 +332,7 @@ applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const ID
 
           return applyBinaryOperation<BinOperatorT>(fh, gh);
         } else {
-          throw NormalError(invalid_operands(fh, g));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(fh, g));
         }
       }
       case 3: {
@@ -339,15 +342,17 @@ applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const ID
 
           return applyBinaryOperation<BinOperatorT>(fh, gh);
         } else {
-          throw NormalError(invalid_operands(fh, g));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(fh, g));
         }
       }
+        // LCOV_EXCL_START
       default: {
         throw UnexpectedError("invalid rhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g));
       }
+        // LCOV_EXCL_STOP
       }
     } else {
-      throw NormalError(invalid_operands(fh, g));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(fh, g));
     }
   }
   case ASTNodeDataType::matrix_t: {
@@ -369,17 +374,21 @@ applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const ID
 
         return applyBinaryOperation<BinOperatorT>(fh, gh);
       }
+        // LCOV_EXCL_START
       default: {
         throw UnexpectedError("invalid rhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g));
       }
+        // LCOV_EXCL_STOP
       }
     } else {
-      throw NormalError(invalid_operands(fh, g));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(fh, g));
     }
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid rhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(g));
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -391,38 +400,43 @@ applyBinaryOperation(const std::shared_ptr<const IDiscreteFunction>& f,
   Assert(f->mesh() == g->mesh());
   Assert(not EmbeddedIDiscreteFunctionUtils::isSameDiscretization(f, g));
 
-  switch (f->dataType()) {
-  case ASTNodeDataType::double_t: {
-    auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*f);
-
-    return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
-  }
-  case ASTNodeDataType::matrix_t: {
-    Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
-    switch (f->dataType().numberOfRows()) {
-    case 1: {
-      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
-
+  if (f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->dataType()) {
+    case ASTNodeDataType::double_t: {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*f);
       return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
     }
-    case 2: {
-      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>&>(*f);
+    case ASTNodeDataType::matrix_t: {
+      Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
+      switch (f->dataType().numberOfRows()) {
+      case 1: {
+        auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
 
-      return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
-    }
-    case 3: {
-      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
+        return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
+      }
+      case 2: {
+        auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>&>(*f);
 
-      return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
+        return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
+      }
+      case 3: {
+        auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
+
+        return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
+      }
+        // LCOV_EXCL_START
+      default: {
+        throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+      }
+        // LCOV_EXCL_STOP
+      }
     }
     default: {
-      throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
     }
     }
-  }
-  default: {
-    throw NormalError(invalid_operands(f, g));
-  }
+  } else {
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -448,9 +462,11 @@ applyBinaryOperation(const std::shared_ptr<const IDiscreteFunction>& f,
   case 3: {
     return applyBinaryOperation<BinOperatorT, 3>(f, g);
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid mesh dimension");
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -460,7 +476,7 @@ operator+(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_p
   if (EmbeddedIDiscreteFunctionUtils::isSameDiscretization(f, g)) {
     return innerCompositionLaw<language::plus_op>(f, g);
   } else {
-    throw NormalError(invalid_operands(f, g));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -470,7 +486,7 @@ operator-(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_p
   if (EmbeddedIDiscreteFunctionUtils::isSameDiscretization(f, g)) {
     return innerCompositionLaw<language::minus_op>(f, g);
   } else {
-    throw NormalError(invalid_operands(f, g));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
   }
 }
 
@@ -508,7 +524,7 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const DiscreteFunctionT&
                          (is_tiny_matrix_v<rhs_data_type> or is_tiny_vector_v<rhs_data_type>)) {
       return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(a, f))>(BinOp<BinOperatorT>{}.eval(a, f));
     } else {
-      throw NormalError(invalid_operands(a, f));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
     }
   } else if constexpr (std::is_same_v<language::plus_op, BinOperatorT> or
                        std::is_same_v<language::minus_op, BinOperatorT>) {
@@ -517,16 +533,18 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const DiscreteFunctionT&
     } else if constexpr (std::is_same_v<lhs_data_type, rhs_data_type>) {
       return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(a, f))>(BinOp<BinOperatorT>{}.eval(a, f));
     } else {
-      throw NormalError(invalid_operands(a, f));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
     }
   } else if constexpr (std::is_same_v<language::divide_op, BinOperatorT>) {
     if constexpr (std::is_same_v<lhs_data_type, double> and std::is_arithmetic_v<rhs_data_type>) {
       return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(a, f))>(BinOp<BinOperatorT>{}.eval(a, f));
     } else {
-      throw NormalError(invalid_operands(a, f));
+      // LCOV_EXCL_START
+      throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
+      // LCOV_EXCL_STOP
     }
   } else {
-    throw NormalError(invalid_operands(a, f));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
 }
 
@@ -544,10 +562,10 @@ applyBinaryOperationToVectorWithLeftConstant(const DataType& a, const DiscreteFu
                          (is_tiny_matrix_v<rhs_data_type> or is_tiny_vector_v<rhs_data_type>)) {
       return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(a, f))>(BinOp<BinOperatorT>{}.eval(a, f));
     } else {
-      throw NormalError(invalid_operands(a, f));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
     }
   } else {
-    throw NormalError(invalid_operands(a, f));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
 }
 
@@ -567,7 +585,9 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
       auto fh = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*f);
       return applyBinaryOperationToVectorWithLeftConstant<BinOperatorT>(a, fh);
     } else {
-      throw NormalError(invalid_operands(a, f));
+      // LCOV_EXCL_START
+      throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
+      // LCOV_EXCL_STOP
     }
   }
   case ASTNodeDataType::vector_t: {
@@ -578,7 +598,7 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
           auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<1>>&>(*f);
           return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
         } else {
-          throw NormalError(invalid_operands(a, f));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
         }
       }
       case 2: {
@@ -586,7 +606,7 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
           auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<2>>&>(*f);
           return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
         } else {
-          throw NormalError(invalid_operands(a, f));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
         }
       }
       case 3: {
@@ -594,12 +614,14 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
           auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>&>(*f);
           return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
         } else {
-          throw NormalError(invalid_operands(a, f));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
         }
       }
+        // LCOV_EXCL_START
       default: {
         throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
       }
+        // LCOV_EXCL_STOP
       }
     } else {
       switch (f->dataType().dimension()) {
@@ -615,9 +637,11 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
         auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>&>(*f);
         return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
       }
+        // LCOV_EXCL_START
       default: {
         throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
       }
+        // LCOV_EXCL_STOP
       }
     }
   }
@@ -630,7 +654,7 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
           auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
           return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
         } else {
-          throw NormalError(invalid_operands(a, f));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
         }
       }
       case 2: {
@@ -638,7 +662,7 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
           auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>&>(*f);
           return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
         } else {
-          throw NormalError(invalid_operands(a, f));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
         }
       }
       case 3: {
@@ -646,12 +670,14 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
           auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
           return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
         } else {
-          throw NormalError(invalid_operands(a, f));
+          throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
         }
       }
+        // LCOV_EXCL_START
       default: {
         throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
       }
+        // LCOV_EXCL_STOP
       }
     } else {
       switch (f->dataType().numberOfRows()) {
@@ -667,15 +693,19 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
         auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
         return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
       }
+        // LCOV_EXCL_START
       default: {
         throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
       }
+        // LCOV_EXCL_STOP
       }
     }
   }
+    // LCOV_EXCL_START
   default: {
-    throw NormalError(invalid_operands(a, f));
+    throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -693,9 +723,11 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
   case 3: {
     return applyBinaryOperationWithLeftConstant<BinOperatorT, 3>(a, f);
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid mesh dimension");
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -710,13 +742,17 @@ applyBinaryOperationWithRightConstant(const DiscreteFunctionT& f, const DataType
 
   if constexpr (std::is_same_v<language::multiply_op, BinOperatorT>) {
     if constexpr (is_tiny_matrix_v<lhs_data_type> and is_tiny_matrix_v<rhs_data_type>) {
-      return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(f, a))>(BinOp<BinOperatorT>{}.eval(f, a));
+      if constexpr (lhs_data_type::NumberOfColumns == rhs_data_type::NumberOfRows) {
+        return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(f, a))>(BinOp<BinOperatorT>{}.eval(f, a));
+      } else {
+        throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
+      }
     } else if constexpr (std::is_same_v<lhs_data_type, double> and
                          (is_tiny_matrix_v<rhs_data_type> or is_tiny_vector_v<rhs_data_type> or
                           std::is_arithmetic_v<rhs_data_type>)) {
       return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(f, a))>(BinOp<BinOperatorT>{}.eval(f, a));
     } else {
-      throw NormalError(invalid_operands(f, a));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
     }
   } else if constexpr (std::is_same_v<language::plus_op, BinOperatorT> or
                        std::is_same_v<language::minus_op, BinOperatorT>) {
@@ -724,10 +760,10 @@ applyBinaryOperationWithRightConstant(const DiscreteFunctionT& f, const DataType
                   (std::is_arithmetic_v<lhs_data_type> and std::is_arithmetic_v<rhs_data_type>)) {
       return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(f, a))>(BinOp<BinOperatorT>{}.eval(f, a));
     } else {
-      throw NormalError(invalid_operands(f, a));
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
     }
   } else {
-    throw NormalError(invalid_operands(f, a));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
 }
 
@@ -736,7 +772,7 @@ std::shared_ptr<const IDiscreteFunction>
 applyBinaryOperationWithRightConstant(const std::shared_ptr<const IDiscreteFunction>& f, const DataType& a)
 {
   if (f->descriptor().type() != DiscreteFunctionType::P0) {
-    throw NormalError(invalid_operands(f, a));
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
 
   switch (f->dataType()) {
@@ -747,61 +783,54 @@ applyBinaryOperationWithRightConstant(const std::shared_ptr<const IDiscreteFunct
     auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*f);
     return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
   }
+  case ASTNodeDataType::vector_t: {
+    switch (f->dataType().dimension()) {
+    case 1: {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<1>>&>(*f);
+      return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
+    }
+    case 2: {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<2>>&>(*f);
+      return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
+    }
+    case 3: {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>&>(*f);
+      return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
   case ASTNodeDataType::matrix_t: {
     Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
-    if constexpr (is_tiny_matrix_v<DataType>) {
-      switch (f->dataType().numberOfRows()) {
-      case 1: {
-        if constexpr (std::is_same_v<DataType, TinyMatrix<1>>) {
-          auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
-          return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
-        } else {
-          throw NormalError(invalid_operands(f, a));
-        }
-      }
-      case 2: {
-        if constexpr (std::is_same_v<DataType, TinyMatrix<2>>) {
-          auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>&>(*f);
-          return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
-        } else {
-          throw NormalError(invalid_operands(f, a));
-        }
-      }
-      case 3: {
-        if constexpr (std::is_same_v<DataType, TinyMatrix<3>>) {
-          auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
-          return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
-        } else {
-          throw NormalError(invalid_operands(f, a));
-        }
-      }
-      default: {
-        throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
-      }
-      }
-    } else {
-      switch (f->dataType().numberOfRows()) {
-      case 1: {
-        auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
-        return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
-      }
-      case 2: {
-        auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>&>(*f);
-        return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
-      }
-      case 3: {
-        auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
-        return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
-      }
-      default: {
-        throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
-      }
-      }
+    switch (f->dataType().numberOfRows()) {
+    case 1: {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
+      return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
+    }
+    case 2: {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>&>(*f);
+      return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
+    }
+    case 3: {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>&>(*f);
+      return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("invalid lhs data type " + EmbeddedIDiscreteFunctionUtils::getOperandTypeName(f));
+    }
+      // LCOV_EXCL_STOP
     }
   }
+    // LCOV_EXCL_START
   default: {
-    throw NormalError(invalid_operands(f, a));
+    throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -819,9 +848,11 @@ applyBinaryOperationWithRightConstant(const std::shared_ptr<const IDiscreteFunct
   case 3: {
     return applyBinaryOperationWithRightConstant<BinOperatorT, 3>(f, a);
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid mesh dimension");
   }
+    // LCOV_EXCL_STOP
   }
 }
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 714479919d1c6eed6d8ef3511694f966143d2b54..25bcaa008af00dd141bfc4a0ba1fc65bb11a2003 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -118,6 +118,7 @@ add_executable (mpi_unit_tests
   test_DiscreteFunctionP0Vector.cpp
   test_DiscreteFunctionVectorInterpoler.cpp
   test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+  test_EmbeddedIDiscreteFunctionOperators.cpp
   test_InterpolateItemArray.cpp
   test_InterpolateItemValue.cpp
   test_ItemArray.cpp
diff --git a/tests/test_EmbeddedIDiscreteFunctionOperators.cpp b/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dc413efce603bf8512d2b1fcaaee35a7a3170fa0
--- /dev/null
+++ b/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
@@ -0,0 +1,2153 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+#include <language/utils/EmbeddedIDiscreteFunctionOperators.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+#define CHECK_SCALAR_VH2_TO_VH(P_LHS, OPERATOR, P_RHS)                                  \
+  {                                                                                     \
+    using DiscreteFunctionType = const std::decay_t<decltype(*P_LHS OPERATOR * P_RHS)>; \
+                                                                                        \
+    std::shared_ptr p_fuv = P_LHS OPERATOR 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) {         \
+      if (fuv[cell_id] != (lhs_values[cell_id] OPERATOR rhs_values[cell_id])) {         \
+        is_same = false;                                                                \
+        break;                                                                          \
+      }                                                                                 \
+    }                                                                                   \
+                                                                                        \
+    REQUIRE(is_same);                                                                   \
+  }
+
+#define CHECK_SCALAR_VHxX_TO_VH(P_LHS, OPERATOR, RHS)                               \
+  {                                                                                 \
+    using DiscreteFunctionType = const std::decay_t<decltype(*P_LHS OPERATOR RHS)>; \
+                                                                                    \
+    std::shared_ptr p_fuv = P_LHS OPERATOR 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) {     \
+      if (fuv[cell_id] != (lhs_values[cell_id] OPERATOR RHS)) {                     \
+        is_same = false;                                                            \
+        break;                                                                      \
+      }                                                                             \
+    }                                                                               \
+                                                                                    \
+    REQUIRE(is_same);                                                               \
+  }
+
+#define CHECK_SCALAR_XxVH_TO_VH(LHS, OPERATOR, P_RHS)                                \
+  {                                                                                  \
+    using DiscreteFunctionType = const std::decay_t<decltype(LHS OPERATOR * P_RHS)>; \
+                                                                                     \
+    std::shared_ptr p_fuv = LHS OPERATOR 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) {      \
+      if (fuv[cell_id] != (LHS OPERATOR rhs_values[cell_id])) {                      \
+        is_same = false;                                                             \
+        break;                                                                       \
+      }                                                                              \
+    }                                                                                \
+                                                                                     \
+    REQUIRE(is_same);                                                                \
+  }
+
+#define CHECK_VECTOR_VH2_TO_VH(P_LHS, OPERATOR, P_RHS)                                     \
+  {                                                                                        \
+    using DiscreteFunctionType = const std::decay_t<decltype(*P_LHS OPERATOR * P_RHS)>;    \
+                                                                                           \
+    std::shared_ptr p_fuv = P_LHS OPERATOR 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->cellArrays();                                                 \
+    auto rhs_values = P_RHS->cellArrays();                                                 \
+    bool is_same    = true;                                                                \
+    for (CellId cell_id = 0; cell_id < lhs_values.numberOfItems(); ++cell_id) {            \
+      for (size_t i = 0; i < fuv.size(); ++i) {                                            \
+        if (fuv[cell_id][i] != (lhs_values[cell_id][i] OPERATOR rhs_values[cell_id][i])) { \
+          is_same = false;                                                                 \
+          break;                                                                           \
+        }                                                                                  \
+      }                                                                                    \
+    }                                                                                      \
+                                                                                           \
+    REQUIRE(is_same);                                                                      \
+  }
+
+#define CHECK_VECTOR_XxVH_TO_VH(LHS, OPERATOR, P_RHS)                                \
+  {                                                                                  \
+    using DiscreteFunctionType = const std::decay_t<decltype(LHS OPERATOR * P_RHS)>; \
+                                                                                     \
+    std::shared_ptr p_fuv = LHS OPERATOR 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->cellArrays();                                           \
+    bool is_same    = true;                                                          \
+    for (CellId cell_id = 0; cell_id < rhs_values.numberOfItems(); ++cell_id) {      \
+      for (size_t i = 0; i < fuv.size(); ++i) {                                      \
+        if (fuv[cell_id][i] != (LHS OPERATOR rhs_values[cell_id][i])) {              \
+          is_same = false;                                                           \
+          break;                                                                     \
+        }                                                                            \
+      }                                                                              \
+    }                                                                                \
+                                                                                     \
+    REQUIRE(is_same);                                                                \
+  }
+
+TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
+{
+  SECTION("binary operators")
+  {
+    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> u_R_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> v_R_values = [=] {
+        CellValue<double> build_values{mesh->connectivity()};
+        parallel_for(
+          build_values.numberOfItems(),
+          PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.6 + std::sin(l2Norm(xj[cell_id])); });
+        return build_values;
+      }();
+
+      std::shared_ptr p_R_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, u_R_values);
+      std::shared_ptr p_other_mesh_R_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, double>>(other_mesh, u_R_values);
+      std::shared_ptr p_R_v = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, v_R_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_other_mesh_R1x1_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(other_mesh, p_R1x1_u->cellValues());
+
+      std::shared_ptr p_R1x1_v = [=] {
+        CellValue<TinyMatrix<1>> vj{mesh->connectivity()};
+        parallel_for(
+          vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = {0.3 - xj[cell_id][0]}; });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, vj);
+      }();
+
+      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_other_mesh_R2x2_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(other_mesh, p_R2x2_u->cellValues());
+
+      std::shared_ptr p_R2x2_v = [=] {
+        CellValue<TinyMatrix<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] + 0.3, 1 - x[1] - x[0],   //
+                           2 * x[1] + x[0], x[1] - x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, vj);
+      }();
+
+      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);
+      }();
+
+      std::shared_ptr p_other_mesh_R3x3_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(other_mesh, p_R3x3_u->cellValues());
+
+      std::shared_ptr p_R3x3_v = [=] {
+        CellValue<TinyMatrix<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] = {0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
+                           2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
+                           2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, vj);
+      }();
+
+      std::shared_ptr p_Vector3_u = [=] {
+        CellArray<double> uj_vector{mesh->connectivity(), 3};
+        parallel_for(
+          uj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            uj_vector[cell_id][0] = 2 * x[0] + 1;
+            uj_vector[cell_id][1] = 1 - x[1] * x[2];
+            uj_vector[cell_id][2] = x[0] + x[2];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, uj_vector);
+      }();
+
+      std::shared_ptr p_other_mesh_Vector3_u =
+        std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(other_mesh, p_Vector3_u->cellArrays());
+
+      std::shared_ptr p_Vector3_v = [=] {
+        CellArray<double> vj_vector{mesh->connectivity(), 3};
+        parallel_for(
+          vj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            vj_vector[cell_id][0] = x[0] * x[1] + 1;
+            vj_vector[cell_id][1] = 2 * x[1];
+            vj_vector[cell_id][2] = x[2] * x[0];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, vj_vector);
+      }();
+
+      std::shared_ptr p_Vector2_w = [=] {
+        CellArray<double> wj_vector{mesh->connectivity(), 2};
+        parallel_for(
+          wj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            wj_vector[cell_id][0] = x[0] + x[1] * 2;
+            wj_vector[cell_id][1] = x[0] * x[1];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, wj_vector);
+      }();
+
+      SECTION("sum")
+      {
+        SECTION("Vh + Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, +, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1_u, +, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2_u, +, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3_u, +, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, +, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, +, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, +, p_R3x3_v);
+
+          CHECK_VECTOR_VH2_TO_VH(p_Vector3_u, +, p_Vector3_v);
+
+          REQUIRE_THROWS_WITH(p_R_u + p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u + p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u + p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u + p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_R_u + p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_R_v, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_Vector2_w, "error: Vh(P0Vector:R) spaces have different sizes");
+
+          REQUIRE_THROWS_WITH(p_R_u + p_other_mesh_R_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1_u + p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2_u + p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3_u + p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u + p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u + p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u + p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh + X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1_u, +, (TinyVector<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2_u, +, (TinyVector<2>{1.2, 2.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3_u, +, (TinyVector<3>{3.2, 7.1, 5.2}));
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, +, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, +, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, +,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R) and R^1");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R) and R^2");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<3>{2, 3, 2}), "error: incompatible operand types Vh(P0:R) and R^3");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u + (double{1}), "error: incompatible operand types Vh(P0Vector:R) and R");
+          REQUIRE_THROWS_WITH(p_Vector3_u + (TinyVector<1>{1}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^1");
+          REQUIRE_THROWS_WITH(p_Vector3_u + (TinyVector<2>{1, 2}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^2");
+        }
+
+        SECTION("X + Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, +, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<1>{1.3}), +, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<2>{1.2, 2.3}), +, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<3>{3.2, 7.1, 5.2}), +, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), +, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), +, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}),
+                                  +, p_R3x3_u);
+
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) + p_R_u, "error: incompatible operand types R^1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) + p_R_u, "error: incompatible operand types R^2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<3>{2, 3, 2}) + p_R_u, "error: incompatible operand types R^3 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) + p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) + p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) + p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((double{1}) + p_Vector3_u, "error: incompatible operand types R and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) + p_Vector3_u,
+                              "error: incompatible operand types R^1 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) + p_Vector3_u,
+                              "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("difference")
+      {
+        SECTION("Vh - Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, -, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1_u, -, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2_u, -, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3_u, -, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, -, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, -, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, -, p_R3x3_v);
+
+          CHECK_VECTOR_VH2_TO_VH(p_Vector3_u, -, p_Vector3_v);
+
+          REQUIRE_THROWS_WITH(p_R_u - p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u - p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u - p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u - p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_R_v, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_Vector2_w, "error: Vh(P0Vector:R) spaces have different sizes");
+
+          REQUIRE_THROWS_WITH(p_R_u - p_other_mesh_R_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1_u - p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2_u - p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3_u - p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u - p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u - p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u - p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh - X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1_u, -, (TinyVector<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2_u, -, (TinyVector<2>{1.2, 2.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3_u, -, (TinyVector<3>{3.2, 7.1, 5.2}));
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, -, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, -, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, -,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R) and R^1");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R) and R^2");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<3>{2, 3, 2}), "error: incompatible operand types Vh(P0:R) and R^3");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u - (double{1}), "error: incompatible operand types Vh(P0Vector:R) and R");
+          REQUIRE_THROWS_WITH(p_Vector3_u - (TinyVector<1>{1}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^1");
+          REQUIRE_THROWS_WITH(p_Vector3_u - (TinyVector<2>{1, 2}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^2");
+        }
+
+        SECTION("X - Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, -, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<1>{1.3}), -, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<2>{1.2, 2.3}), -, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<3>{3.2, 7.1, 5.2}), -, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), -, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), -, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}),
+                                  -, p_R3x3_u);
+
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) - p_R_u, "error: incompatible operand types R^1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) - p_R_u, "error: incompatible operand types R^2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<3>{2, 3, 2}) - p_R_u, "error: incompatible operand types R^3 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) - p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) - p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) - p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((double{1}) - p_Vector3_u, "error: incompatible operand types R and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) - p_Vector3_u,
+                              "error: incompatible operand types R^1 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) - p_Vector3_u,
+                              "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("product")
+      {
+        SECTION("Vh * Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, *, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, *, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, *, p_R3x3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R3x3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, *, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, *, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, *, p_R3_v);
+
+          {
+            std::shared_ptr p_fuv = p_R_u * p_Vector3_v;
+
+            REQUIRE(p_fuv.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_fuv));
+
+            const auto& fuv = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_fuv);
+
+            auto lhs_values = p_R_u->cellValues();
+            auto rhs_arrays = p_Vector3_v->cellArrays();
+            bool is_same    = true;
+            for (CellId cell_id = 0; cell_id < lhs_values.numberOfItems(); ++cell_id) {
+              for (size_t i = 0; i < fuv.size(); ++i) {
+                if (fuv[cell_id][i] != (lhs_values[cell_id] * rhs_arrays[cell_id][i])) {
+                  is_same = false;
+                  break;
+                }
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(p_R1_u * p_R1_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u * p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u * p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R1_u * p_R2x2_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^2x2)");
+
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_R2x2_v, "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_R3x3_v, "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0:R^3x3)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_R1x1_v, "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0:R^1x1)");
+
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_R2_v, "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_R3_v, "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_R1_v, "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0:R^1)");
+
+          REQUIRE_THROWS_WITH(p_R1_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R2_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R3_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0Vector:R)");
+
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R1_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R2_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R3_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R1x1_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R2x2_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R3x3_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^3x3)");
+
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh * X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, *, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, *, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, *,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R1_u * (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R^1) and R^1");
+          REQUIRE_THROWS_WITH(p_R2_u * (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R^2) and R^2");
+          REQUIRE_THROWS_WITH(p_R3_u * (TinyVector<3>{2, 3, 2}),
+                              "error: incompatible operand types Vh(P0:R^3) and R^3");
+          REQUIRE_THROWS_WITH(p_R1_u * (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R^1) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R2_u * (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R^2) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R3_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R^3) and R^3x3");
+          REQUIRE_THROWS_WITH(p_R2x2_u * (TinyMatrix<1>{2}),
+                              "error: incompatible operand types Vh(P0:R^2x2) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R1x1_u * (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R^1x1) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R2x2_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R^2x2) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^3x3");
+          REQUIRE_THROWS_WITH(p_Vector3_u * (double{2}), "error: incompatible operand types Vh(P0Vector:R) and R");
+        }
+
+        SECTION("X * Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R1x1_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R2x2_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R3x3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), *, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), *, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                                     4.7, 2.3, 7.1,   //
+                                                                     9.7, 3.2, 6.8}),
+                                                      *, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                                     4.7, 2.3, 7.1,   //
+                                                                     9.7, 3.2, 6.8}),
+                                                      *, p_R3x3_u);
+
+          CHECK_VECTOR_XxVH_TO_VH(bool{true}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(uint64_t{1}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(int64_t{2}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(double{1.3}, *, p_Vector3_u);
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R2_u, "error: incompatible operand types R^1x1 and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R3_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R2_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R1_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^1)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R2x2_u,
+                              "error: incompatible operand types R^1x1 and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R3x3_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^3x3)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R2x2_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R1x1_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^1x1)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_Vector3_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_Vector3_u,
+                              "error: incompatible operand types R^1x1 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("ratio")
+      {
+        SECTION("Vh / Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, /, p_R_v);
+
+          REQUIRE_THROWS_WITH(p_R_u / p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u / p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u / p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u / p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+
+          REQUIRE_THROWS_WITH(p_R_u / p_other_mesh_R_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("X / Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, /, p_R_u);
+        }
+      }
+    }
+
+    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> u_R_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> v_R_values = [=] {
+        CellValue<double> build_values{mesh->connectivity()};
+        parallel_for(
+          build_values.numberOfItems(),
+          PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.6 + std::sin(l2Norm(xj[cell_id])); });
+        return build_values;
+      }();
+
+      std::shared_ptr p_R_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, u_R_values);
+      std::shared_ptr p_other_mesh_R_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, double>>(other_mesh, u_R_values);
+      std::shared_ptr p_R_v = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, v_R_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_other_mesh_R1x1_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(other_mesh, p_R1x1_u->cellValues());
+
+      std::shared_ptr p_R1x1_v = [=] {
+        CellValue<TinyMatrix<1>> vj{mesh->connectivity()};
+        parallel_for(
+          vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = {0.3 - xj[cell_id][0]}; });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, vj);
+      }();
+
+      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_other_mesh_R2x2_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(other_mesh, p_R2x2_u->cellValues());
+
+      std::shared_ptr p_R2x2_v = [=] {
+        CellValue<TinyMatrix<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] + 0.3, 1 - x[1] - x[0],   //
+                           2 * x[1] + x[0], x[1] - x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, vj);
+      }();
+
+      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);
+      }();
+
+      std::shared_ptr p_other_mesh_R3x3_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(other_mesh, p_R3x3_u->cellValues());
+
+      std::shared_ptr p_R3x3_v = [=] {
+        CellValue<TinyMatrix<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] = {0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
+                           2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
+                           2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, vj);
+      }();
+
+      std::shared_ptr p_Vector3_u = [=] {
+        CellArray<double> uj_vector{mesh->connectivity(), 3};
+        parallel_for(
+          uj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            uj_vector[cell_id][0] = 2 * x[0] + 1;
+            uj_vector[cell_id][1] = 1 - x[1] * x[2];
+            uj_vector[cell_id][2] = x[0] + x[2];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, uj_vector);
+      }();
+
+      std::shared_ptr p_other_mesh_Vector3_u =
+        std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(other_mesh, p_Vector3_u->cellArrays());
+
+      std::shared_ptr p_Vector3_v = [=] {
+        CellArray<double> vj_vector{mesh->connectivity(), 3};
+        parallel_for(
+          vj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            vj_vector[cell_id][0] = x[0] * x[1] + 1;
+            vj_vector[cell_id][1] = 2 * x[1];
+            vj_vector[cell_id][2] = x[2] * x[0];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, vj_vector);
+      }();
+
+      std::shared_ptr p_Vector2_w = [=] {
+        CellArray<double> wj_vector{mesh->connectivity(), 2};
+        parallel_for(
+          wj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            wj_vector[cell_id][0] = x[0] + x[1] * 2;
+            wj_vector[cell_id][1] = x[0] * x[1];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, wj_vector);
+      }();
+
+      SECTION("sum")
+      {
+        SECTION("Vh + Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, +, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1_u, +, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2_u, +, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3_u, +, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, +, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, +, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, +, p_R3x3_v);
+
+          CHECK_VECTOR_VH2_TO_VH(p_Vector3_u, +, p_Vector3_v);
+
+          REQUIRE_THROWS_WITH(p_R_u + p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u + p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u + p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u + p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_R_u + p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_R_v, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_Vector2_w, "error: Vh(P0Vector:R) spaces have different sizes");
+
+          REQUIRE_THROWS_WITH(p_R_u + p_other_mesh_R_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1_u + p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2_u + p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3_u + p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u + p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u + p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u + p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh + X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1_u, +, (TinyVector<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2_u, +, (TinyVector<2>{1.2, 2.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3_u, +, (TinyVector<3>{3.2, 7.1, 5.2}));
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, +, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, +, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, +,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R) and R^1");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R) and R^2");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<3>{2, 3, 2}), "error: incompatible operand types Vh(P0:R) and R^3");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u + (double{1}), "error: incompatible operand types Vh(P0Vector:R) and R");
+          REQUIRE_THROWS_WITH(p_Vector3_u + (TinyVector<1>{1}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^1");
+          REQUIRE_THROWS_WITH(p_Vector3_u + (TinyVector<2>{1, 2}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^2");
+        }
+
+        SECTION("X + Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, +, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<1>{1.3}), +, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<2>{1.2, 2.3}), +, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<3>{3.2, 7.1, 5.2}), +, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), +, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), +, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}),
+                                  +, p_R3x3_u);
+
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) + p_R_u, "error: incompatible operand types R^1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) + p_R_u, "error: incompatible operand types R^2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<3>{2, 3, 2}) + p_R_u, "error: incompatible operand types R^3 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) + p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) + p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) + p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((double{1}) + p_Vector3_u, "error: incompatible operand types R and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) + p_Vector3_u,
+                              "error: incompatible operand types R^1 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) + p_Vector3_u,
+                              "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("difference")
+      {
+        SECTION("Vh - Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, -, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1_u, -, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2_u, -, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3_u, -, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, -, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, -, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, -, p_R3x3_v);
+
+          CHECK_VECTOR_VH2_TO_VH(p_Vector3_u, -, p_Vector3_v);
+
+          REQUIRE_THROWS_WITH(p_R_u - p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u - p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u - p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u - p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_R_v, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_Vector2_w, "error: Vh(P0Vector:R) spaces have different sizes");
+
+          REQUIRE_THROWS_WITH(p_R_u - p_other_mesh_R_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1_u - p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2_u - p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3_u - p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u - p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u - p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u - p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh - X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1_u, -, (TinyVector<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2_u, -, (TinyVector<2>{1.2, 2.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3_u, -, (TinyVector<3>{3.2, 7.1, 5.2}));
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, -, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, -, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, -,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R) and R^1");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R) and R^2");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<3>{2, 3, 2}), "error: incompatible operand types Vh(P0:R) and R^3");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u - (double{1}), "error: incompatible operand types Vh(P0Vector:R) and R");
+          REQUIRE_THROWS_WITH(p_Vector3_u - (TinyVector<1>{1}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^1");
+          REQUIRE_THROWS_WITH(p_Vector3_u - (TinyVector<2>{1, 2}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^2");
+        }
+
+        SECTION("X - Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, -, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<1>{1.3}), -, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<2>{1.2, 2.3}), -, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<3>{3.2, 7.1, 5.2}), -, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), -, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), -, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}),
+                                  -, p_R3x3_u);
+
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) - p_R_u, "error: incompatible operand types R^1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) - p_R_u, "error: incompatible operand types R^2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<3>{2, 3, 2}) - p_R_u, "error: incompatible operand types R^3 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) - p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) - p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) - p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((double{1}) - p_Vector3_u, "error: incompatible operand types R and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) - p_Vector3_u,
+                              "error: incompatible operand types R^1 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) - p_Vector3_u,
+                              "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("product")
+      {
+        SECTION("Vh * Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, *, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, *, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, *, p_R3x3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R3x3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, *, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, *, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, *, p_R3_v);
+
+          {
+            std::shared_ptr p_fuv = p_R_u * p_Vector3_v;
+
+            REQUIRE(p_fuv.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_fuv));
+
+            const auto& fuv = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_fuv);
+
+            auto lhs_values = p_R_u->cellValues();
+            auto rhs_arrays = p_Vector3_v->cellArrays();
+            bool is_same    = true;
+            for (CellId cell_id = 0; cell_id < lhs_values.numberOfItems(); ++cell_id) {
+              for (size_t i = 0; i < fuv.size(); ++i) {
+                if (fuv[cell_id][i] != (lhs_values[cell_id] * rhs_arrays[cell_id][i])) {
+                  is_same = false;
+                  break;
+                }
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(p_R1_u * p_R1_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u * p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u * p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R1_u * p_R2x2_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^2x2)");
+
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_R2x2_v, "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_R3x3_v, "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0:R^3x3)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_R1x1_v, "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0:R^1x1)");
+
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_R2_v, "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_R3_v, "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_R1_v, "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0:R^1)");
+
+          REQUIRE_THROWS_WITH(p_R1_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R2_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R3_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0Vector:R)");
+
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R1_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R2_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R3_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R1x1_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R2x2_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R3x3_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^3x3)");
+
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh * X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, *, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, *, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, *,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R1_u * (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R^1) and R^1");
+          REQUIRE_THROWS_WITH(p_R2_u * (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R^2) and R^2");
+          REQUIRE_THROWS_WITH(p_R3_u * (TinyVector<3>{2, 3, 2}),
+                              "error: incompatible operand types Vh(P0:R^3) and R^3");
+          REQUIRE_THROWS_WITH(p_R1_u * (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R^1) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R2_u * (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R^2) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R3_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R^3) and R^3x3");
+          REQUIRE_THROWS_WITH(p_R2x2_u * (TinyMatrix<1>{2}),
+                              "error: incompatible operand types Vh(P0:R^2x2) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R1x1_u * (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R^1x1) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R2x2_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R^2x2) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^3x3");
+          REQUIRE_THROWS_WITH(p_Vector3_u * (double{2}), "error: incompatible operand types Vh(P0Vector:R) and R");
+        }
+
+        SECTION("X * Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R1x1_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R2x2_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R3x3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), *, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), *, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                                     4.7, 2.3, 7.1,   //
+                                                                     9.7, 3.2, 6.8}),
+                                                      *, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                                     4.7, 2.3, 7.1,   //
+                                                                     9.7, 3.2, 6.8}),
+                                                      *, p_R3x3_u);
+
+          CHECK_VECTOR_XxVH_TO_VH(bool{true}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(uint64_t{1}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(int64_t{2}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(double{1.3}, *, p_Vector3_u);
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R2_u, "error: incompatible operand types R^1x1 and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R3_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R2_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R1_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^1)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R2x2_u,
+                              "error: incompatible operand types R^1x1 and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R3x3_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^3x3)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R2x2_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R1x1_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^1x1)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_Vector3_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_Vector3_u,
+                              "error: incompatible operand types R^1x1 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("ratio")
+      {
+        SECTION("Vh / Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, /, p_R_v);
+
+          REQUIRE_THROWS_WITH(p_R_u / p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u / p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u / p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u / p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+
+          REQUIRE_THROWS_WITH(p_R_u / p_other_mesh_R_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("X / Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, /, p_R_u);
+        }
+      }
+    }
+
+    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> u_R_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> v_R_values = [=] {
+        CellValue<double> build_values{mesh->connectivity()};
+        parallel_for(
+          build_values.numberOfItems(),
+          PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.6 + std::sin(l2Norm(xj[cell_id])); });
+        return build_values;
+      }();
+
+      std::shared_ptr p_R_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, u_R_values);
+      std::shared_ptr p_other_mesh_R_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, double>>(other_mesh, u_R_values);
+      std::shared_ptr p_R_v = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, v_R_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_other_mesh_R1x1_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(other_mesh, p_R1x1_u->cellValues());
+
+      std::shared_ptr p_R1x1_v = [=] {
+        CellValue<TinyMatrix<1>> vj{mesh->connectivity()};
+        parallel_for(
+          vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = {0.3 - xj[cell_id][0]}; });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, vj);
+      }();
+
+      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_other_mesh_R2x2_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(other_mesh, p_R2x2_u->cellValues());
+
+      std::shared_ptr p_R2x2_v = [=] {
+        CellValue<TinyMatrix<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] + 0.3, 1 - x[1] - x[0],   //
+                           2 * x[1] + x[0], x[1] - x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, vj);
+      }();
+
+      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);
+      }();
+
+      std::shared_ptr p_other_mesh_R3x3_u =
+        std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(other_mesh, p_R3x3_u->cellValues());
+
+      std::shared_ptr p_R3x3_v = [=] {
+        CellValue<TinyMatrix<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] = {0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
+                           2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
+                           2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, vj);
+      }();
+
+      std::shared_ptr p_Vector3_u = [=] {
+        CellArray<double> uj_vector{mesh->connectivity(), 3};
+        parallel_for(
+          uj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            uj_vector[cell_id][0] = 2 * x[0] + 1;
+            uj_vector[cell_id][1] = 1 - x[1] * x[2];
+            uj_vector[cell_id][2] = x[0] + x[2];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, uj_vector);
+      }();
+
+      std::shared_ptr p_other_mesh_Vector3_u =
+        std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(other_mesh, p_Vector3_u->cellArrays());
+
+      std::shared_ptr p_Vector3_v = [=] {
+        CellArray<double> vj_vector{mesh->connectivity(), 3};
+        parallel_for(
+          vj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            vj_vector[cell_id][0] = x[0] * x[1] + 1;
+            vj_vector[cell_id][1] = 2 * x[1];
+            vj_vector[cell_id][2] = x[2] * x[0];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, vj_vector);
+      }();
+
+      std::shared_ptr p_Vector2_w = [=] {
+        CellArray<double> wj_vector{mesh->connectivity(), 2};
+        parallel_for(
+          wj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<3> x = to_3d(xj[cell_id]);
+            wj_vector[cell_id][0] = x[0] + x[1] * 2;
+            wj_vector[cell_id][1] = x[0] * x[1];
+          });
+
+        return std::make_shared<const DiscreteFunctionP0Vector<Dimension, double>>(mesh, wj_vector);
+      }();
+
+      SECTION("sum")
+      {
+        SECTION("Vh + Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, +, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1_u, +, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2_u, +, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3_u, +, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, +, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, +, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, +, p_R3x3_v);
+
+          CHECK_VECTOR_VH2_TO_VH(p_Vector3_u, +, p_Vector3_v);
+
+          REQUIRE_THROWS_WITH(p_R_u + p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u + p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u + p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u + p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_R_u + p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_R_v, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_Vector2_w, "error: Vh(P0Vector:R) spaces have different sizes");
+
+          REQUIRE_THROWS_WITH(p_R_u + p_other_mesh_R_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1_u + p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2_u + p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3_u + p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u + p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u + p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u + p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_Vector3_u + p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh + X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, +, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1_u, +, (TinyVector<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2_u, +, (TinyVector<2>{1.2, 2.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3_u, +, (TinyVector<3>{3.2, 7.1, 5.2}));
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, +, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, +, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, +,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R) and R^1");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R) and R^2");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyVector<3>{2, 3, 2}), "error: incompatible operand types Vh(P0:R) and R^3");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R_u + (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u + (double{1}), "error: incompatible operand types Vh(P0Vector:R) and R");
+          REQUIRE_THROWS_WITH(p_Vector3_u + (TinyVector<1>{1}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^1");
+          REQUIRE_THROWS_WITH(p_Vector3_u + (TinyVector<2>{1, 2}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^2");
+        }
+
+        SECTION("X + Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, +, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, +, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<1>{1.3}), +, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<2>{1.2, 2.3}), +, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<3>{3.2, 7.1, 5.2}), +, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), +, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), +, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}),
+                                  +, p_R3x3_u);
+
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) + p_R_u, "error: incompatible operand types R^1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) + p_R_u, "error: incompatible operand types R^2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<3>{2, 3, 2}) + p_R_u, "error: incompatible operand types R^3 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) + p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) + p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) + p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((double{1}) + p_Vector3_u, "error: incompatible operand types R and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) + p_Vector3_u,
+                              "error: incompatible operand types R^1 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) + p_Vector3_u,
+                              "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("difference")
+      {
+        SECTION("Vh - Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, -, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1_u, -, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2_u, -, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3_u, -, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, -, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, -, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, -, p_R3x3_v);
+
+          CHECK_VECTOR_VH2_TO_VH(p_Vector3_u, -, p_Vector3_v);
+
+          REQUIRE_THROWS_WITH(p_R_u - p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u - p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u - p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u - p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_R_v, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_Vector2_w, "error: Vh(P0Vector:R) spaces have different sizes");
+
+          REQUIRE_THROWS_WITH(p_R_u - p_other_mesh_R_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1_u - p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2_u - p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3_u - p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u - p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u - p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u - p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_Vector3_u - p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh - X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, -, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1_u, -, (TinyVector<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2_u, -, (TinyVector<2>{1.2, 2.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3_u, -, (TinyVector<3>{3.2, 7.1, 5.2}));
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, -, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, -, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, -,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R) and R^1");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R) and R^2");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyVector<3>{2, 3, 2}), "error: incompatible operand types Vh(P0:R) and R^3");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R_u - (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u - (double{1}), "error: incompatible operand types Vh(P0Vector:R) and R");
+          REQUIRE_THROWS_WITH(p_Vector3_u - (TinyVector<1>{1}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^1");
+          REQUIRE_THROWS_WITH(p_Vector3_u - (TinyVector<2>{1, 2}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^2");
+        }
+
+        SECTION("X - Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, -, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, -, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<1>{1.3}), -, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<2>{1.2, 2.3}), -, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyVector<3>{3.2, 7.1, 5.2}), -, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), -, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), -, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}),
+                                  -, p_R3x3_u);
+
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) - p_R_u, "error: incompatible operand types R^1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) - p_R_u, "error: incompatible operand types R^2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyVector<3>{2, 3, 2}) - p_R_u, "error: incompatible operand types R^3 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) - p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) - p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) - p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((double{1}) - p_Vector3_u, "error: incompatible operand types R and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<1>{1}) - p_Vector3_u,
+                              "error: incompatible operand types R^1 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyVector<2>{1, 2}) - p_Vector3_u,
+                              "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("product")
+      {
+        SECTION("Vh * Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, *, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, *, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, *, p_R3x3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R1x1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R2x2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, *, p_R3x3_v);
+
+          CHECK_SCALAR_VH2_TO_VH(p_R1x1_u, *, p_R1_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R2x2_u, *, p_R2_v);
+          CHECK_SCALAR_VH2_TO_VH(p_R3x3_u, *, p_R3_v);
+
+          {
+            std::shared_ptr p_fuv = p_R_u * p_Vector3_v;
+
+            REQUIRE(p_fuv.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_fuv));
+
+            const auto& fuv = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_fuv);
+
+            auto lhs_values = p_R_u->cellValues();
+            auto rhs_arrays = p_Vector3_v->cellArrays();
+            bool is_same    = true;
+            for (CellId cell_id = 0; cell_id < lhs_values.numberOfItems(); ++cell_id) {
+              for (size_t i = 0; i < fuv.size(); ++i) {
+                if (fuv[cell_id][i] != (lhs_values[cell_id] * rhs_arrays[cell_id][i])) {
+                  is_same = false;
+                  break;
+                }
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(p_R1_u * p_R1_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u * p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u * p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R1_u * p_R2x2_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^2x2)");
+
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_R2x2_v, "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_R3x3_v, "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0:R^3x3)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_R1x1_v, "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0:R^1x1)");
+
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_R2_v, "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_R3_v, "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_R1_v, "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0:R^1)");
+
+          REQUIRE_THROWS_WITH(p_R1_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^1) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R2_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R3_u * p_Vector3_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^2x2) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0:R^3x3) and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_u * p_Vector3_v,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0Vector:R)");
+
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R1_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R2_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R3_u, "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R1x1_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R2x2_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH(p_Vector3_v * p_R3x3_u,
+                              "error: incompatible operand types Vh(P0Vector:R) and Vh(P0:R^3x3)");
+
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R1x1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R2x2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_R3x3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R1x1_u * p_other_mesh_R1_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R2x2_u * p_other_mesh_R2_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R3x3_u * p_other_mesh_R3_u, "error: operands are defined on different meshes");
+          REQUIRE_THROWS_WITH(p_R_u * p_other_mesh_Vector3_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("Vh * X -> Vh")
+        {
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, bool{true});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, uint64_t{1});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, int64_t{2});
+          CHECK_SCALAR_VHxX_TO_VH(p_R_u, *, double{1.3});
+
+          CHECK_SCALAR_VHxX_TO_VH(p_R1x1_u, *, (TinyMatrix<1>{1.3}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R2x2_u, *, (TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}));
+          CHECK_SCALAR_VHxX_TO_VH(p_R3x3_u, *,
+                                  (TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                 4.7, 2.3, 7.1,   //
+                                                 9.7, 3.2, 6.8}));
+
+          REQUIRE_THROWS_WITH(p_R1_u * (TinyVector<1>{1}), "error: incompatible operand types Vh(P0:R^1) and R^1");
+          REQUIRE_THROWS_WITH(p_R2_u * (TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0:R^2) and R^2");
+          REQUIRE_THROWS_WITH(p_R3_u * (TinyVector<3>{2, 3, 2}),
+                              "error: incompatible operand types Vh(P0:R^3) and R^3");
+          REQUIRE_THROWS_WITH(p_R1_u * (TinyMatrix<1>{2}), "error: incompatible operand types Vh(P0:R^1) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R2_u * (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R^2) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R3_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R^3) and R^3x3");
+          REQUIRE_THROWS_WITH(p_R2x2_u * (TinyMatrix<1>{2}),
+                              "error: incompatible operand types Vh(P0:R^2x2) and R^1x1");
+          REQUIRE_THROWS_WITH(p_R1x1_u * (TinyMatrix<2>{2, 3, 1, 4}),
+                              "error: incompatible operand types Vh(P0:R^1x1) and R^2x2");
+          REQUIRE_THROWS_WITH(p_R2x2_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0:R^2x2) and R^3x3");
+
+          REQUIRE_THROWS_WITH(p_Vector3_u * (TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}),
+                              "error: incompatible operand types Vh(P0Vector:R) and R^3x3");
+          REQUIRE_THROWS_WITH(p_Vector3_u * (double{2}), "error: incompatible operand types Vh(P0Vector:R) and R");
+        }
+
+        SECTION("X * Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R1x1_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R2x2_u);
+
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, *, p_R3x3_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, *, p_R3x3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), *, p_R1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), *, p_R2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                                     4.7, 2.3, 7.1,   //
+                                                                     9.7, 3.2, 6.8}),
+                                                      *, p_R3_u);
+
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<1>{1.3}), *, p_R1x1_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<2>{1.2, 2.3, 4.2, 5.1}), *, p_R2x2_u);
+          CHECK_SCALAR_XxVH_TO_VH((TinyMatrix<3>{3.2, 7.1, 5.2,   //
+                                                                     4.7, 2.3, 7.1,   //
+                                                                     9.7, 3.2, 6.8}),
+                                                      *, p_R3x3_u);
+
+          CHECK_VECTOR_XxVH_TO_VH(bool{true}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(uint64_t{1}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(int64_t{2}, *, p_Vector3_u);
+          CHECK_VECTOR_XxVH_TO_VH(double{1.3}, *, p_Vector3_u);
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R_u, "error: incompatible operand types R^1x1 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R2_u, "error: incompatible operand types R^1x1 and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R3_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^3)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R2_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R1_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^1)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_R2x2_u,
+                              "error: incompatible operand types R^1x1 and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R3x3_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^3x3)");
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_R2x2_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0:R^2x2)");
+          REQUIRE_THROWS_WITH((TinyMatrix<2>{2, 3, 1, 4}) * p_R1x1_u,
+                              "error: incompatible operand types R^2x2 and Vh(P0:R^1x1)");
+
+          REQUIRE_THROWS_WITH((TinyMatrix<3>{1, 3, 6, 4, 7, 2, 5, 9, 8}) * p_Vector3_u,
+                              "error: incompatible operand types R^3x3 and Vh(P0Vector:R)");
+          REQUIRE_THROWS_WITH((TinyMatrix<1>{2}) * p_Vector3_u,
+                              "error: incompatible operand types R^1x1 and Vh(P0Vector:R)");
+        }
+      }
+
+      SECTION("ratio")
+      {
+        SECTION("Vh / Vh -> Vh")
+        {
+          CHECK_SCALAR_VH2_TO_VH(p_R_u, /, p_R_v);
+
+          REQUIRE_THROWS_WITH(p_R_u / p_R1_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R2_u / p_R1_v, "error: incompatible operand types Vh(P0:R^2) and Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(p_R3_u / p_R1x1_v, "error: incompatible operand types Vh(P0:R^3) and Vh(P0:R^1x1)");
+          REQUIRE_THROWS_WITH(p_R_u / p_R2x2_v, "error: incompatible operand types Vh(P0:R) and Vh(P0:R^2x2)");
+
+          REQUIRE_THROWS_WITH(p_R_u / p_other_mesh_R_u, "error: operands are defined on different meshes");
+        }
+
+        SECTION("X / Vh -> Vh")
+        {
+          CHECK_SCALAR_XxVH_TO_VH(bool{true}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(uint64_t{1}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(int64_t{2}, /, p_R_u);
+          CHECK_SCALAR_XxVH_TO_VH(double{1.3}, /, p_R_u);
+        }
+      }
+    }
+  }
+}