diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index 68b22b552f932a1fb48baeb96d6742a74a4717c2..da49a03833ff1f0a5e9aa19c345855da12b1ce85 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -147,4 +147,54 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
                                                                     std::shared_ptr<const IDiscreteFunction>)>>(
                            [](std::shared_ptr<const IDiscreteFunction> a, std::shared_ptr<const IDiscreteFunction> b)
                              -> std::shared_ptr<const IDiscreteFunction> { return pow(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("dot",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>,
+                                                                    std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("dot",
+                         std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                           const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, const TinyVector<1>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a, const TinyVector<1> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("dot",
+                         std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                           const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, const TinyVector<2>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a, const TinyVector<2> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("dot",
+                         std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                           const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, const TinyVector<3>&)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a, const TinyVector<3>& b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("dot",
+                         std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                           const IDiscreteFunction>(const TinyVector<1>, std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](const TinyVector<1> a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("dot",
+                         std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                           const IDiscreteFunction>(const TinyVector<2>, std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](const TinyVector<2> a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("dot",
+                         std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                           const IDiscreteFunction>(const TinyVector<3>&, std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](const TinyVector<3>& a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); }));
 }
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index 2939e011af950fc414cde730e2c1206bf8f4ae01..e785033cbc27a2a52081dc1fd20efd8f4ad2d097 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -27,7 +27,7 @@
     }                                                                                                                 \
     }                                                                                                                 \
   } else {                                                                                                            \
-    throw UnexpectedError("invalid operand type " + operand_type_name(ARG));                                          \
+    throw NormalError("invalid operand type " + operand_type_name(ARG));                                              \
   }
 
 std::shared_ptr<const IDiscreteFunction>
@@ -82,7 +82,7 @@ std::shared_ptr<const IDiscreteFunction>
 atan2(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
 {
   if ((f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) and
-      (f->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
+      (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
     std::shared_ptr mesh = getCommonMesh({f, g});
 
     if (mesh.use_count() == 0) {
@@ -111,7 +111,7 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<c
     }
   } else {
     std::stringstream os;
-    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g);
+    os << "incompatible operands type " << operand_type_name(f) << " and " << operand_type_name(g);
     throw NormalError(os.str());
   }
 }
@@ -138,7 +138,9 @@ atan2(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
     }
     }
   } else {
-    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+    std::stringstream os;
+    os << "incompatible operands type " << operand_type_name(a) << " and " << operand_type_name(f);
+    throw NormalError(os.str());
   }
 }
 
@@ -164,7 +166,9 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
     }
     }
   } else {
-    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+    std::stringstream os;
+    os << "incompatible operands type " << operand_type_name(f) << " and " << operand_type_name(a);
+    throw NormalError(os.str());
   }
 }
 
@@ -220,7 +224,7 @@ std::shared_ptr<const IDiscreteFunction>
 pow(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
 {
   if ((f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) and
-      (f->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
+      (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
     std::shared_ptr mesh = getCommonMesh({f, g});
 
     if (mesh.use_count() == 0) {
@@ -276,7 +280,9 @@ pow(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
     }
     }
   } else {
-    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f);
+    throw NormalError(os.str());
   }
 }
 
@@ -302,6 +308,174 @@ pow(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
     }
     }
   } else {
-    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a);
+    throw NormalError(os.str());
+  }
+}
+
+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()));
+
+  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())) {
+    std::shared_ptr mesh = getCommonMesh({f, g});
+
+    if (mesh.use_count() == 0) {
+      throw NormalError("operands are defined on different meshes");
+    }
+
+    switch (mesh->dimension()) {
+    case 1: {
+      return dot<1>(f, g);
+    }
+    case 2: {
+      return dot<2>(f, g);
+    }
+    case 3: {
+      return dot<3>(f, g);
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g);
+    throw NormalError(os.str());
+  }
+}
+
+template <size_t Dimension, size_t VectorDimension>
+std::shared_ptr<const IDiscreteFunction>
+dot(const std::shared_ptr<const IDiscreteFunction>& f, const TinyVector<VectorDimension>& a)
+{
+  Assert((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+         (f->dataType().dimension() == a.dimension()));
+
+  using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+  using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>>;
+
+  return std::make_shared<const DiscreteFunctionResultType>(dot(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+}
+
+template <size_t VectorDimension>
+std::shared_ptr<const IDiscreteFunction>
+dot(const std::shared_ptr<const IDiscreteFunction>& f, const TinyVector<VectorDimension>& a)
+{
+  if ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+      (f->dataType().dimension() == a.dimension())) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      return dot<1, VectorDimension>(f, a);
+    }
+    case 2: {
+      return dot<2, VectorDimension>(f, a);
+    }
+    case 3: {
+      return dot<3, VectorDimension>(f, a);
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a);
+    throw NormalError(os.str());
   }
 }
+
+template <size_t Dimension, size_t VectorDimension>
+std::shared_ptr<const IDiscreteFunction>
+dot(const TinyVector<VectorDimension>& a, const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  Assert((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+         (f->dataType().dimension() == a.dimension()));
+
+  using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+  using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, TinyVector<VectorDimension>>;
+
+  return std::make_shared<const DiscreteFunctionResultType>(dot(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+}
+
+template <size_t VectorDimension>
+std::shared_ptr<const IDiscreteFunction>
+dot(const TinyVector<VectorDimension>& a, const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+      (f->dataType().dimension() == a.dimension())) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      return dot<1, VectorDimension>(a, f);
+    }
+    case 2: {
+      return dot<2, VectorDimension>(a, f);
+    }
+    case 3: {
+      return dot<3, VectorDimension>(a, f);
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    std::stringstream os;
+    os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f);
+    throw NormalError(os.str());
+  }
+}
+
+template std::shared_ptr<const IDiscreteFunction> dot(const std::shared_ptr<const IDiscreteFunction>&,
+                                                      const TinyVector<1>&);
+
+template std::shared_ptr<const IDiscreteFunction> dot(const std::shared_ptr<const IDiscreteFunction>&,
+                                                      const TinyVector<2>&);
+
+template std::shared_ptr<const IDiscreteFunction> dot(const std::shared_ptr<const IDiscreteFunction>&,
+                                                      const TinyVector<3>&);
+
+template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<1>&,
+                                                      const std::shared_ptr<const IDiscreteFunction>&);
+
+template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<2>&,
+                                                      const std::shared_ptr<const IDiscreteFunction>&);
+
+template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<3>&,
+                                                      const std::shared_ptr<const IDiscreteFunction>&);
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
index 0030f46759dfe683a09228131140e98ec205fc8b..51ece8b824ebbdd8004b06d61544cd506f77fce3 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -54,4 +54,15 @@ std::shared_ptr<const IDiscreteFunction> pow(const double, const std::shared_ptr
 
 std::shared_ptr<const IDiscreteFunction> pow(const std::shared_ptr<const IDiscreteFunction>&, const double);
 
+std::shared_ptr<const IDiscreteFunction> dot(const std::shared_ptr<const IDiscreteFunction>&,
+                                             const std::shared_ptr<const IDiscreteFunction>&);
+
+template <size_t VectorDimension>
+std::shared_ptr<const IDiscreteFunction> dot(const std::shared_ptr<const IDiscreteFunction>&,
+                                             const TinyVector<VectorDimension>&);
+
+template <size_t VectorDimension>
+std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&,
+                                             const std::shared_ptr<const IDiscreteFunction>&);
+
 #endif   // EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 0b77aaea481451c6662dd450756affc02ce0226d..1d9493a8bdc5032d81fb5e41d6d0ba5ee4b17603 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -523,6 +523,42 @@ class DiscreteFunctionP0 : public IDiscreteFunction
     return result;
   }
 
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  dot(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
+  {
+    static_assert(is_tiny_vector_v<DataType>);
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(f[cell_id], g[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  dot(const DiscreteFunctionP0& f, const DataType& a)
+  {
+    static_assert(is_tiny_vector_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(f[cell_id], a); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  dot(const DataType& a, const DiscreteFunctionP0& f)
+  {
+    static_assert(is_tiny_vector_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(a, f[cell_id]); });
+
+    return result;
+  }
+
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()} {}
 
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value)