diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index a688acf1372f93b10192fec99abe0bd15de52683..d4d10c6c7c531add6cd9f32c2ad376fcd97100a6 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -3,6 +3,7 @@
 #include <language/utils/EmbeddedIDiscreteFunctionUtils.hpp>
 #include <mesh/IMesh.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
@@ -318,43 +319,63 @@ template <size_t Dimension>
 std::shared_ptr<const IDiscreteFunction>
 dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
 {
-  Assert((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
-         (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
-         (f->dataType().dimension() == g->dataType().dimension()));
+  Assert(((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+          (g->descriptor().type() == DiscreteFunctionType::P0Vector)) or
+         ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+          (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
+          (f->dataType().dimension() == g->dataType().dimension())));
 
-  using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+  if ((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+      (g->descriptor().type() == DiscreteFunctionType::P0Vector)) {
+    using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+    using DiscreteFunctionType       = DiscreteFunctionP0Vector<Dimension, double>;
 
-  switch (f->dataType().dimension()) {
-  case 1: {
-    using Rd                   = TinyVector<1>;
-    using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
-    return std::make_shared<const DiscreteFunctionResultType>(
-      dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
-  }
-  case 2: {
-    using Rd                   = TinyVector<2>;
-    using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
-    return std::make_shared<const DiscreteFunctionResultType>(
-      dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
-  }
-  case 3: {
-    using Rd                   = TinyVector<3>;
-    using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
-    return std::make_shared<const DiscreteFunctionResultType>(
-      dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
-  }
-  default: {
-    throw UnexpectedError("invalid data dimension " + operand_type_name(f));
-  }
+    const DiscreteFunctionType& f_vector = dynamic_cast<const DiscreteFunctionType&>(*f);
+    const DiscreteFunctionType& g_vector = dynamic_cast<const DiscreteFunctionType&>(*g);
+
+    if (f_vector.size() != g_vector.size()) {
+      throw NormalError("operands have different dimension");
+    } else {
+      return std::make_shared<const DiscreteFunctionResultType>(dot(f_vector, g_vector));
+    }
+
+  } else {
+    using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+
+    switch (f->dataType().dimension()) {
+    case 1: {
+      using Rd                   = TinyVector<1>;
+      using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+      return std::make_shared<const DiscreteFunctionResultType>(
+        dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 2: {
+      using Rd                   = TinyVector<2>;
+      using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+      return std::make_shared<const DiscreteFunctionResultType>(
+        dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 3: {
+      using Rd                   = TinyVector<3>;
+      using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+      return std::make_shared<const DiscreteFunctionResultType>(
+        dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    default: {
+      throw UnexpectedError("invalid data dimension " + operand_type_name(f));
+    }
+    }
   }
 }
 
 std::shared_ptr<const IDiscreteFunction>
 dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
 {
-  if ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
-      (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
-      (f->dataType().dimension() == g->dataType().dimension())) {
+  if (((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+       (g->descriptor().type() == DiscreteFunctionType::P0Vector)) or
+      ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+       (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
+       (f->dataType().dimension() == g->dataType().dimension()))) {
     std::shared_ptr mesh = getCommonMesh({f, g});
 
     if (mesh.use_count() == 0) {
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index 103718b125e123692241715e697ab89869d078eb..5f00d5964d3b0e060549a6337b62e40fe2b97f21 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -3,6 +3,7 @@
 
 #include <scheme/IDiscreteFunction.hpp>
 
+#include <algebra/Vector.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemArray.hpp>
 #include <mesh/Mesh.hpp>
@@ -194,6 +195,19 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
     return product;
   }
 
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  dot(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
+  {
+    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
+    Assert(f.size() == g.size());
+    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(),
+      PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(Vector{f[cell_id]}, Vector{g[cell_id]}); });
+
+    return result;
+  }
+
   DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh, size_t size)
     : m_mesh{mesh}, m_cell_arrays{mesh->connectivity(), size}
   {}