diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index fa746269722229bff7fbd15f47edf926f14784cc..f49c4c033251bf46d7cd91e7c524af4d146db9d3 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -812,7 +812,7 @@ integral_of(const std::shared_ptr<const IDiscreteFunction>& f)
     default: {
       throw UnexpectedError("invalid mesh dimension");
     }
-      // LCOV_EXCL_START
+      // LCOV_EXCL_STOP
     }
   } else {
     throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
diff --git a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
index 7cf6866e96a6a9141be5771bd93768a1a6d1d415..d41be6566f5bcbc1205c5561abf474f02780b6af 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -5,6 +5,7 @@
 
 #include <language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -269,6 +270,44 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
       return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
     }();
 
+    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_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("sqrt Vh -> Vh")
     {
       CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, sqrt);
@@ -464,10 +503,30 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, dot);
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, dot);
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, dot);
+
+      {
+        auto p_UV = dot(p_Vector3_u, p_Vector3_v);
+        REQUIRE(p_UV.use_count() == 1);
+
+        auto UV        = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_UV);
+        auto direct_UV = dot(*p_Vector3_u, *p_Vector3_v);
+
+        bool is_same = true;
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          if (UV[cell_id] != direct_UV[cell_id]) {
+            is_same = false;
+            break;
+          }
+        }
+
+        REQUIRE(is_same);
+      }
+
       REQUIRE_THROWS_WITH(dot(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
     }
 
     SECTION("dot Vh*Rd -> Vh")
@@ -708,6 +767,44 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
       return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
     }();
 
+    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_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("sqrt Vh -> Vh")
     {
       CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, sqrt);
@@ -903,10 +1000,30 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, dot);
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, dot);
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, dot);
+
+      {
+        auto p_UV = dot(p_Vector3_u, p_Vector3_v);
+        REQUIRE(p_UV.use_count() == 1);
+
+        auto UV        = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_UV);
+        auto direct_UV = dot(*p_Vector3_u, *p_Vector3_v);
+
+        bool is_same = true;
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          if (UV[cell_id] != direct_UV[cell_id]) {
+            is_same = false;
+            break;
+          }
+        }
+
+        REQUIRE(is_same);
+      }
+
       REQUIRE_THROWS_WITH(dot(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
     }
 
     SECTION("dot Vh*Rd -> Vh")
@@ -1147,6 +1264,44 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
       return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
     }();
 
+    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_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("sqrt Vh -> Vh")
     {
       CHECK_EMBEDDED_VH_TO_VH_REAL_FUNCTION_EVALUATION(p_positive_u, sqrt);
@@ -1342,10 +1497,30 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, dot);
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, dot);
       CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, dot);
+
+      {
+        auto p_UV = dot(p_Vector3_u, p_Vector3_v);
+        REQUIRE(p_UV.use_count() == 1);
+
+        auto UV        = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_UV);
+        auto direct_UV = dot(*p_Vector3_u, *p_Vector3_v);
+
+        bool is_same = true;
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          if (UV[cell_id] != direct_UV[cell_id]) {
+            is_same = false;
+            break;
+          }
+        }
+
+        REQUIRE(is_same);
+      }
+
       REQUIRE_THROWS_WITH(dot(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes");
       REQUIRE_THROWS_WITH(dot(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+      REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
     }
 
     SECTION("dot Vh*Rd -> Vh")