diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index ff34b89e148d844215cf7fffc21c2ad78bca5ec9..dc62e370e6ec304cb2e4fd44d220218d9bd94283 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -27,7 +27,7 @@ struct BiCGStab
 
     r_k_1 = b - A * x;
 
-    double residu = std::sqrt((r_k_1, r_k_1));   // Norm(r_k_1);
+    double residu = std::sqrt(dot(r_k_1, r_k_1));   // Norm(r_k_1);
 
     if (residu != 0) {
       double resid0 = residu;
@@ -53,23 +53,23 @@ struct BiCGStab
 
         Ap_k = A * p_k;
 
-        const double alpha_k = (r_k_1, rTilda_0) / (Ap_k, rTilda_0);
+        const double alpha_k = dot(r_k_1, rTilda_0) / dot(Ap_k, rTilda_0);
 
         s_k  = r_k_1 - alpha_k * Ap_k;
         As_k = A * s_k;
 
-        const double w_k = (As_k, s_k) / (As_k, As_k);
+        const double w_k = dot(As_k, s_k) / dot(As_k, As_k);
 
         x += alpha_k * p_k + w_k * s_k;
         r_k = s_k - w_k * As_k;
 
-        const double beta_k = (r_k, rTilda_0) / (r_k_1, rTilda_0) * (alpha_k / w_k);
+        const double beta_k = dot(r_k, rTilda_0) / dot(r_k_1, rTilda_0) * (alpha_k / w_k);
 
         p_k -= w_k * Ap_k;
         p_k *= beta_k;
         p_k += r_k;
 
-        if ((residu = std::sqrt((r_k, r_k))) / resid0 < epsilon) {
+        if ((residu = std::sqrt(dot(r_k, r_k))) / resid0 < epsilon) {
           break;
         }
 
diff --git a/src/algebra/CG.hpp b/src/algebra/CG.hpp
index 538c77571dac4d551d7476ea79bf855e8026fb53..7e653d8f95cddab4d5d1b7eb36481fd7690ef170 100644
--- a/src/algebra/CG.hpp
+++ b/src/algebra/CG.hpp
@@ -28,7 +28,7 @@ struct CG
     if (verbose) {
       h = A * x;
       h -= f;
-      std::cout << "- initial " << rang::style::bold << "real" << rang::style::reset << " residu :   " << (h, h)
+      std::cout << "- initial " << rang::style::bold << "real" << rang::style::reset << " residu :   " << dot(h, h)
                 << '\n';
     }
 
@@ -48,14 +48,14 @@ struct CG
 
         g = copy(cg);   // TODO: precond: g = g/C
 
-        gcg = (g, cg);
+        gcg = dot(g, cg);
 
         h = copy(g);
       }
 
       b = A * h;
 
-      double hAh = (h, b);
+      double hAh = dot(h, b);
 
       if (hAh == 0) {
         hAh = 1.;
@@ -69,7 +69,7 @@ struct CG
       g -= ro * b;
 
       double gamma = gcg;
-      gcg          = (g, cg);
+      gcg          = dot(g, cg);
 
       if ((i == 1) && (gcg != 0)) {
         relativeEpsilon = epsilon * gcg;
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index c3146a2be9932353388d34d15862056f952631fb..f24503aade42275e8b96f6958535fc8439da356c 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -64,11 +64,11 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr T operator,(const TinyVector& v) const
+  constexpr friend T dot(const TinyVector& u, const TinyVector& v)
   {
-    T t = m_values[0] * v.m_values[0];
+    T t = u.m_values[0] * v.m_values[0];
     for (size_t i = 1; i < N; ++i) {
-      t += m_values[i] * v.m_values[i];
+      t += u.m_values[i] * v.m_values[i];
     }
     return t;
   }
@@ -239,7 +239,7 @@ l2Norm(const TinyVector<N, T>& x)
 {
   static_assert(std::is_arithmetic<T>(), "Cannot compute L2 norm for non-arithmetic types");
   static_assert(std::is_floating_point<T>::value, "L2 norm is defined for floating point types only");
-  return std::sqrt((x, x));
+  return std::sqrt(dot(x, x));
 }
 
 // Cross product is only defined for K^3 vectors
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index efe64b41a8940be9fadf7066dd5dc27bff90a633..089efa633161f4d3e83d16a5fc666ad5f9fc1696 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -33,18 +33,18 @@ class Vector   // LCOV_EXCL_LINE
     return x_copy;
   }
 
-  friend Vector operator*(const DataType& a, const Vector& x)
+  friend Vector
+  operator*(const DataType& a, const Vector& x)
   {
     Vector<std::remove_const_t<DataType>> y = copy(x);
     return y *= a;
   }
 
   template <typename DataType2>
-  PUGS_INLINE
-  auto
-  operator,(const Vector<DataType2>& y)
+  PUGS_INLINE friend auto
+  dot(const Vector& x, const Vector<DataType2>& y)
   {
-    Assert(this->size() == y.size());
+    Assert(x.size() == y.size());
     // Quite ugly, TODO: fix me in C++20
     auto promoted = [] {
       DataType a{0};
@@ -55,8 +55,8 @@ class Vector   // LCOV_EXCL_LINE
     decltype(promoted) sum = 0;
 
     // Does not use parallel_for to preserve sum order
-    for (index_type i = 0; i < this->size(); ++i) {
-      sum += m_values[i] * y.m_values[i];
+    for (index_type i = 0; i < x.size(); ++i) {
+      sum += x.m_values[i] * y.m_values[i];
     }
 
     return sum;
@@ -130,7 +130,8 @@ class Vector   // LCOV_EXCL_LINE
   }
 
   PUGS_INLINE
-  DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+  DataType&
+  operator[](index_type i) const noexcept(NO_ASSERT)
   {
     return m_values[i];
   }
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index 68b22b552f932a1fb48baeb96d6742a74a4717c2..7617463648fdee030760ef65ed9dec50bdb21232 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -147,4 +147,155 @@ 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); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a) -> double { return min(a); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         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 min(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(double, std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](double a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("min",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, double)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a,
+                              double b) -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a) -> double { return max(a); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         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 max(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(double, std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](double a, std::shared_ptr<const IDiscreteFunction> b)
+                             -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("max",
+                         std::make_shared<BuiltinFunctionEmbedder<
+                           std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, double)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a,
+                              double b) -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("sum_to_R",
+                         std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>(
+                           [](std::shared_ptr<const IDiscreteFunction> a) -> double { return sum_to<double>(a); }));
+
+  scheme_module._addBuiltinFunction("sum_to_R1",
+                                    std::make_shared<
+                                      BuiltinFunctionEmbedder<TinyVector<1>(std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<1> {
+                                        return sum_to<TinyVector<1>>(a);
+                                      }));
+
+  scheme_module._addBuiltinFunction("sum_to_R2",
+                                    std::make_shared<
+                                      BuiltinFunctionEmbedder<TinyVector<2>(std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<2> {
+                                        return sum_to<TinyVector<2>>(a);
+                                      }));
+
+  scheme_module._addBuiltinFunction("sum_to_R3",
+                                    std::make_shared<
+                                      BuiltinFunctionEmbedder<TinyVector<3>(std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<3> {
+                                        return sum_to<TinyVector<3>>(a);
+                                      }));
+
+  scheme_module._addBuiltinFunction("sum_to_R1x1",
+                                    std::make_shared<
+                                      BuiltinFunctionEmbedder<TinyMatrix<1>(std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<1> {
+                                        return sum_to<TinyMatrix<1>>(a);
+                                      }));
+
+  scheme_module._addBuiltinFunction("sum_to_R2x2",
+                                    std::make_shared<
+                                      BuiltinFunctionEmbedder<TinyMatrix<2>(std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<2> {
+                                        return sum_to<TinyMatrix<2>>(a);
+                                      }));
+
+  scheme_module._addBuiltinFunction("sum_to_R3x3",
+                                    std::make_shared<
+                                      BuiltinFunctionEmbedder<TinyMatrix<3>(std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<3> {
+                                        return sum_to<TinyMatrix<3>>(a);
+                                      }));
 }
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index aa961e4a0ae250581fad42b644f7d670ef4fdb32..69f6890fbc5e8c3e3515901acd5b5102100e8c35 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -69,6 +69,19 @@ MathModule::MathModule()
 
   this->_addBuiltinFunction("round", std::make_shared<BuiltinFunctionEmbedder<int64_t(double)>>(
                                        [](double x) -> int64_t { return std::lround(x); }));
+
+  this->_addBuiltinFunction("dot",
+                            std::make_shared<BuiltinFunctionEmbedder<double(const TinyVector<1>, const TinyVector<1>)>>(
+                              [](const TinyVector<1> x, const TinyVector<1> y) -> double { return dot(x, y); }));
+
+  this->_addBuiltinFunction("dot",
+                            std::make_shared<BuiltinFunctionEmbedder<double(const TinyVector<2>, const TinyVector<2>)>>(
+                              [](const TinyVector<2> x, const TinyVector<2> y) -> double { return dot(x, y); }));
+
+  this
+    ->_addBuiltinFunction("dot",
+                          std::make_shared<BuiltinFunctionEmbedder<double(const TinyVector<3>&, const TinyVector<3>&)>>(
+                            [](const TinyVector<3>& x, const TinyVector<3>& y) -> double { return dot(x, y); }));
 }
 
 void
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index d6f9227f4b66e1f1abc442ca36f2cc21a8b9b04b..b4022d660475cf93eb9f9e1cb291622de22dba22 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -6,12 +6,16 @@
 #include <language/utils/BinaryOperatorProcessorBuilder.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/TypeDescriptor.hpp>
+#include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
 #include <scheme/AcousticSolver.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
 #include <scheme/DiscreteFunctionInterpoler.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
@@ -25,8 +29,6 @@
 #include <scheme/NumberedBoundaryDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 
-#include <scheme/PerfectGas.hpp>
-
 #include <memory>
 
 SchemeModule::SchemeModule()
@@ -150,82 +152,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this
-    ->_addBuiltinFunction("perfect_gas_epsilon_from_rho_p_gamma",
-                          std::make_shared<BuiltinFunctionEmbedder<
-                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
-
-                            [](const std::shared_ptr<const IDiscreteFunction>& rho,
-                               const std::shared_ptr<const IDiscreteFunction>& p,
-                               const std::shared_ptr<const IDiscreteFunction>& gamma)
-                              -> std::shared_ptr<const IDiscreteFunction> {
-                              return perfect_gas::epsilonFromRhoPGamma(rho, p, gamma);
-                            }
-
-                            ));
-
-  this
-    ->_addBuiltinFunction("perfect_gas_p_from_rho_epsilon_gamma",
-                          std::make_shared<BuiltinFunctionEmbedder<
-                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
-
-                            [](const std::shared_ptr<const IDiscreteFunction>& rho,
-                               const std::shared_ptr<const IDiscreteFunction>& epsilon,
-                               const std::shared_ptr<const IDiscreteFunction>& gamma)
-                              -> std::shared_ptr<const IDiscreteFunction> {
-                              return perfect_gas::pFromRhoEpsilonGamma(rho, epsilon, gamma);
-                            }
-
-                            ));
-
-  this
-    ->_addBuiltinFunction("perfect_gas_sound_speed",
-                          std::make_shared<BuiltinFunctionEmbedder<
-                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
-
-                            [](const std::shared_ptr<const IDiscreteFunction>& rho,
-                               const std::shared_ptr<const IDiscreteFunction>& p,
-                               const std::shared_ptr<const IDiscreteFunction>& gamma)
-                              -> std::shared_ptr<const IDiscreteFunction> {
-                              return perfect_gas::soundSpeed(rho, p, gamma);
-                            }
-
-                            ));
-
-  this
-    ->_addBuiltinFunction("E_from_epsilon_u",
-                          std::make_shared<BuiltinFunctionEmbedder<
-                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
-
-                            [](const std::shared_ptr<const IDiscreteFunction>& epsilon,
-                               const std::shared_ptr<const IDiscreteFunction>& u)
-                              -> std::shared_ptr<const IDiscreteFunction> {
-                              return perfect_gas::totalEnergyFromEpsilonU(epsilon, u);
-                            }
-
-                            ));
-
-  this
-    ->_addBuiltinFunction("epsilon_from_E_u",
-                          std::make_shared<BuiltinFunctionEmbedder<
-                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
-                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
-
-                            [](const std::shared_ptr<const IDiscreteFunction>& E,
-                               const std::shared_ptr<const IDiscreteFunction>& u)
-                              -> std::shared_ptr<const IDiscreteFunction> {
-                              return perfect_gas::epsilonFromTotalEnergyU(E, u);
-                            }
-
-                            ));
-
   this->_addBuiltinFunction("glace_solver",
                             std::make_shared<BuiltinFunctionEmbedder<std::tuple<
                               std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
@@ -318,6 +244,48 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this
+    ->_addBuiltinFunction("cell_volume",
+                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                            const std::shared_ptr<const IMesh>&)>>(
+
+                            [](const std::shared_ptr<const IMesh>& i_mesh) -> std::shared_ptr<const IDiscreteFunction> {
+                              switch (i_mesh->dimension()) {
+                              case 1: {
+                                constexpr size_t Dimension = 1;
+                                using MeshType             = Mesh<Connectivity<Dimension>>;
+                                std::shared_ptr<const MeshType> mesh =
+                                  std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh);
+
+                                return std::make_shared<const DiscreteFunctionP0<
+                                  Dimension, double>>(mesh, copy(MeshDataManager::instance().getMeshData(*mesh).Vj()));
+                              }
+                              case 2: {
+                                constexpr size_t Dimension = 2;
+                                using MeshType             = Mesh<Connectivity<Dimension>>;
+                                std::shared_ptr<const MeshType> mesh =
+                                  std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh);
+
+                                return std::make_shared<const DiscreteFunctionP0<
+                                  Dimension, double>>(mesh, copy(MeshDataManager::instance().getMeshData(*mesh).Vj()));
+                              }
+                              case 3: {
+                                constexpr size_t Dimension = 3;
+                                using MeshType             = Mesh<Connectivity<Dimension>>;
+                                std::shared_ptr<const MeshType> mesh =
+                                  std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh);
+
+                                return std::make_shared<const DiscreteFunctionP0<
+                                  Dimension, double>>(mesh, copy(MeshDataManager::instance().getMeshData(*mesh).Vj()));
+                              }
+                              default: {
+                                throw UnexpectedError("invalid mesh dimension");
+                              }
+                              }
+                            }
+
+                            ));
+
   MathFunctionRegisterForVh{*this};
 }
 
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index 2939e011af950fc414cde730e2c1206bf8f4ae01..0c5f23fb63210a8f72ddda25688fc758f743c13c 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -9,7 +9,7 @@
 
 #define DISCRETE_FUNCTION_CALL(FUNCTION, ARG)                                                                         \
   if (ARG->dataType() == ASTNodeDataType::double_t and ARG->descriptor().type() == DiscreteFunctionType::P0) {        \
-    switch (f->mesh()->dimension()) {                                                                                 \
+    switch (ARG->mesh()->dimension()) {                                                                               \
     case 1: {                                                                                                         \
       using DiscreteFunctionType = DiscreteFunctionP0<1, double>;                                                     \
       return std::make_shared<const DiscreteFunctionType>(FUNCTION(dynamic_cast<const DiscreteFunctionType&>(*ARG))); \
@@ -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) {
@@ -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 operand types " << 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 operand types " << 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,455 @@ 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>&);
+
+double
+min(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return min(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return min(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return min(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw NormalError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+min(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
+      (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
+    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: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*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());
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+min(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*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());
   }
 }
+
+std::shared_ptr<const IDiscreteFunction>
+min(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*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());
+  }
+}
+
+double
+max(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return max(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return max(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return max(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw NormalError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+max(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
+      (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) {
+    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: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*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());
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+max(const double a, const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*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());
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+max(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
+{
+  if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*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 <typename ValueT>
+ValueT
+sum_to(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->dataType() == ast_node_data_type_from<ValueT> and f->descriptor().type() == DiscreteFunctionType::P0) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      using DiscreteFunctionType = DiscreteFunctionP0<1, ValueT>;
+      return sum(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, ValueT>;
+      return sum(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, ValueT>;
+      return sum(dynamic_cast<const DiscreteFunctionType&>(*f));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw NormalError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+template double sum_to<double>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyVector<1> sum_to<TinyVector<1>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyVector<2> sum_to<TinyVector<2>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyVector<3> sum_to<TinyVector<3>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyMatrix<1> sum_to<TinyMatrix<1>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyMatrix<2> sum_to<TinyMatrix<2>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyMatrix<3> sum_to<TinyMatrix<3>>(const std::shared_ptr<const IDiscreteFunction>&);
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
index 0030f46759dfe683a09228131140e98ec205fc8b..6eef58101016b95411613677cef2d27e5e52eac1 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -54,4 +54,36 @@ 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>&);
+
+double min(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> min(const std::shared_ptr<const IDiscreteFunction>&,
+                                             const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> min(const double, const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> min(const std::shared_ptr<const IDiscreteFunction>&, const double);
+
+double max(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> max(const std::shared_ptr<const IDiscreteFunction>&,
+                                             const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> max(const double, const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> max(const std::shared_ptr<const IDiscreteFunction>&, const double);
+
+template <typename ValueT>
+ValueT sum_to(const std::shared_ptr<const IDiscreteFunction>&);
+
 #endif   // EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 9866c9ad8389c509565ba26c8e94915c3f597931..3a83ef5c3f49d068fbfd8e62d4ff002df50a8845 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -1233,12 +1233,12 @@ GmshReader::__proceedData()
   dimension3_mask[12]                 = 1;
   dimension3_mask[13]                 = 1;
 
-  std::cout << "- dimension 0 entities: " << (dimension0_mask, elementNumber) << '\n';
-  std::cout << "- dimension 1 entities: " << (dimension1_mask, elementNumber) << '\n';
-  std::cout << "- dimension 2 entities: " << (dimension2_mask, elementNumber) << '\n';
-  std::cout << "- dimension 3 entities: " << (dimension3_mask, elementNumber) << '\n';
-  if ((dimension3_mask, elementNumber) > 0) {
-    const size_t nb_cells = (dimension3_mask, elementNumber);
+  std::cout << "- dimension 0 entities: " << dot(dimension0_mask, elementNumber) << '\n';
+  std::cout << "- dimension 1 entities: " << dot(dimension1_mask, elementNumber) << '\n';
+  std::cout << "- dimension 2 entities: " << dot(dimension2_mask, elementNumber) << '\n';
+  std::cout << "- dimension 3 entities: " << dot(dimension3_mask, elementNumber) << '\n';
+  if (dot(dimension3_mask, elementNumber) > 0) {
+    const size_t nb_cells = dot(dimension3_mask, elementNumber);
 
     GmshConnectivityBuilder<3> connectivity_builder(m_mesh_data, nb_cells);
 
@@ -1257,8 +1257,8 @@ GmshReader::__proceedData()
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
 
-  } else if ((dimension2_mask, elementNumber) > 0) {
-    const size_t nb_cells = (dimension2_mask, elementNumber);
+  } else if (dot(dimension2_mask, elementNumber) > 0) {
+    const size_t nb_cells = dot(dimension2_mask, elementNumber);
 
     GmshConnectivityBuilder<2> connectivity_builder(m_mesh_data, nb_cells);
 
@@ -1276,8 +1276,8 @@ GmshReader::__proceedData()
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
 
-  } else if ((dimension1_mask, elementNumber) > 0) {
-    const size_t nb_cells = (dimension1_mask, elementNumber);
+  } else if (dot(dimension1_mask, elementNumber) > 0) {
+    const size_t nb_cells = dot(dimension1_mask, elementNumber);
 
     GmshConnectivityBuilder<1> connectivity_builder(m_mesh_data, nb_cells);
 
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index f891f5ccb40e92cfbd0aa9e4cf3a600f352a8790..c15bb1a66ddacc34fe392a08d133e6f683b0e757 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -265,7 +265,7 @@ class MeshData : public IMeshData
                 const Rd xjxr0 = xr0 - xj;
                 const Rd xjxr1 = xr1 - xj;
 
-                const double Vjlr = (crossProduct(xjxr0, xjxr1), xjxl);
+                const double Vjlr = dot(crossProduct(xjxr0, xjxr1), xjxl);
 
                 sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1);
               }
@@ -297,7 +297,7 @@ class MeshData : public IMeshData
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         for (size_t R = 0; R < cell_nodes.size(); ++R) {
-          sum_cjr_xr += (xr[cell_nodes[R]], Cjr(j, R));
+          sum_cjr_xr += dot(xr[cell_nodes[R]], Cjr(j, R));
         }
         Vj[j] = inv_Dimension * sum_cjr_xr;
       });
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index a49db8b4b358f52e7dd199a64a55a43b1123de2e..92d354fc52ddfd3e35ab69a8d892edb32f6432de 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -157,7 +157,7 @@ MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
   parallel_for(
     m_node_list.size(), PUGS_LAMBDA(size_t r) {
       const R2& x = xr[m_node_list[r]];
-      if ((x - origin, normal) > 1E-13 * length) {
+      if (dot(x - origin, normal) > 1E-13 * length) {
         throw NormalError("this FlatBoundary is not flat!");
       }
     });
@@ -330,13 +330,13 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
   const R3 w = zmax - zmin;
 
   const R3 uv        = crossProduct(u, v);
-  const double uv_l2 = (uv, uv);
+  const double uv_l2 = dot(uv, uv);
 
   R3 normal        = uv;
   double normal_l2 = uv_l2;
 
   const R3 uw        = crossProduct(u, w);
-  const double uw_l2 = (uw, uw);
+  const double uw_l2 = dot(uw, uw);
 
   if (uw_l2 > uv_l2) {
     normal    = uw;
@@ -344,7 +344,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
   }
 
   const R3 vw        = crossProduct(v, w);
-  const double vw_l2 = (vw, vw);
+  const double vw_l2 = dot(vw, vw);
 
   if (vw_l2 > normal_l2) {
     normal    = vw;
@@ -385,7 +385,7 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType
     const auto& j0_nodes = cell_to_node_matrix[j0];
 
     for (size_t r = 0; r < j0_nodes.size(); ++r) {
-      const double height = (xr[j0_nodes[r]] - xr[r0], normal);
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
       if (std::abs(height) > std::abs(max_height)) {
         max_height = height;
       }
@@ -429,7 +429,7 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType
     const CellId j0      = node_to_cell_matrix[r0][0];
     const auto& j0_nodes = cell_to_node_matrix[j0];
     for (size_t r = 0; r < j0_nodes.size(); ++r) {
-      const double height = (xr[j0_nodes[r]] - xr[r0], normal);
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
       if (std::abs(height) > std::abs(max_height)) {
         max_height = height;
       }
@@ -474,7 +474,7 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType
     const auto& j0_nodes = cell_to_node_matrix[j0];
 
     for (size_t r = 0; r < j0_nodes.size(); ++r) {
-      const double height = (xr[j0_nodes[r]] - xr[r0], normal);
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
       if (std::abs(height) > std::abs(max_height)) {
         max_height = height;
       }
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 6bac5a6b3d920f9f3b5654e826bf9b974a9fcd0f..6973037571b23f861cec5c70ed28fbe67c036fcf 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -453,7 +453,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
         for (size_t R = 0; R < cell_nodes.size(); ++R) {
           const NodeId r = cell_nodes[R];
           momentum_fluxes += m_Fjr(j, R);
-          energy_fluxes += (m_Fjr(j, R), m_ur[r]);
+          energy_fluxes += dot(m_Fjr(j, R), m_ur[r]);
         }
         const double dt_over_Mj = dt / (rho[j] * Vj[j]);
         new_u[j] -= dt_over_Mj * momentum_fluxes;
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index b9a24259f7fcf800f4215f7b291115c0e110f41b..2131f2fec0748407ea7bb21a3375d398e1f6242a 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -5,5 +5,4 @@ add_library(
   AcousticSolver.cpp
   DiscreteFunctionInterpoler.cpp
   DiscreteFunctionUtils.cpp
-  DiscreteFunctionVectorInterpoler.cpp
-  PerfectGas.cpp)
+  DiscreteFunctionVectorInterpoler.cpp)
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 0b77aaea481451c6662dd450756affc02ce0226d..ccedb957d9e1fdb9b0c6643d66618709a9717cad 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -4,6 +4,7 @@
 #include <scheme/IDiscreteFunction.hpp>
 
 #include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -523,6 +524,141 @@ 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;
+  }
+
+  PUGS_INLINE friend double
+  min(const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+
+    return min(f.m_cell_values);
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  min(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
+    Assert(f.m_mesh == g.m_mesh);
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(f[cell_id], g[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  min(const double a, const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(a, f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  min(const DiscreteFunctionP0& f, const double a)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(f[cell_id], a); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend double
+  max(const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+
+    return max(f.m_cell_values);
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  max(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
+    Assert(f.m_mesh == g.m_mesh);
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(f[cell_id], g[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  max(const double a, const DiscreteFunctionP0& f)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(a, f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  max(const DiscreteFunctionP0& f, const double a)
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(f.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(f[cell_id], a); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DataType
+  sum(const DiscreteFunctionP0& f)
+  {
+    Assert(f.m_cell_values.isBuilt());
+    return sum(f.m_cell_values);
+  }
+
   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)
diff --git a/src/scheme/PerfectGas.cpp b/src/scheme/PerfectGas.cpp
deleted file mode 100644
index b371dea3c9e4aac8ee867dd0f938b72bcd02791e..0000000000000000000000000000000000000000
--- a/src/scheme/PerfectGas.cpp
+++ /dev/null
@@ -1,248 +0,0 @@
-#include <scheme/PerfectGas.hpp>
-
-#include <scheme/DiscreteFunctionP0.hpp>
-#include <scheme/DiscreteFunctionUtils.hpp>
-
-namespace perfect_gas
-{
-template <size_t Dimension, typename DataType>
-std::shared_ptr<IDiscreteFunction>
-soundSpeed(const DiscreteFunctionP0<Dimension, DataType>& rho,
-           const DiscreteFunctionP0<Dimension, DataType>& p,
-           const DiscreteFunctionP0<Dimension, DataType>& gamma)
-{
-  using MeshType       = Mesh<Connectivity<Dimension>>;
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
-
-  DiscreteFunctionP0<Dimension, DataType> sound_speed{mesh};
-
-  parallel_for(
-    mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { sound_speed[cell_id] = std::sqrt(gamma[cell_id] * p[cell_id] / rho[cell_id]); });
-
-  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(sound_speed);
-}
-
-std::shared_ptr<IDiscreteFunction>
-soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho,
-           const std::shared_ptr<const IDiscreteFunction>& p,
-           const std::shared_ptr<const IDiscreteFunction>& gamma)
-{
-  std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma});
-  if (not mesh.use_count()) {
-    throw NormalError("rho, p and gamma are not defined on the same mesh");
-  }
-
-  switch (mesh->dimension()) {
-  case 1: {
-    return soundSpeed(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
-                      dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p),
-                      dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
-  }
-  case 2: {
-    return soundSpeed(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
-                      dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p),
-                      dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
-  }
-  case 3: {
-    return soundSpeed(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
-                      dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p),
-                      dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-template <size_t Dimension, typename DataType>
-std::shared_ptr<IDiscreteFunction>
-epsilonFromRhoPGamma(const DiscreteFunctionP0<Dimension, DataType>& rho,
-                     const DiscreteFunctionP0<Dimension, DataType>& p,
-                     const DiscreteFunctionP0<Dimension, DataType>& gamma)
-{
-  using MeshType       = Mesh<Connectivity<Dimension>>;
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
-
-  DiscreteFunctionP0<Dimension, DataType> epsilon{mesh};
-
-  parallel_for(
-    mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = p[cell_id] / ((gamma[cell_id] - 1) * rho[cell_id]); });
-
-  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
-}
-
-std::shared_ptr<IDiscreteFunction>
-epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
-                     const std::shared_ptr<const IDiscreteFunction>& p,
-                     const std::shared_ptr<const IDiscreteFunction>& gamma)
-{
-  std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma});
-  if (not mesh.use_count()) {
-    throw NormalError("rho, p and gamma are not defined on the same mesh");
-  }
-
-  switch (mesh->dimension()) {
-  case 1: {
-    return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
-                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p),
-                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
-  }
-  case 2: {
-    return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
-                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p),
-                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
-  }
-  case 3: {
-    return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
-                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p),
-                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-template <size_t Dimension, typename DataType>
-std::shared_ptr<IDiscreteFunction>
-pFromRhoEpsilonGamma(const DiscreteFunctionP0<Dimension, DataType>& rho,
-                     const DiscreteFunctionP0<Dimension, DataType>& epsilon,
-                     const DiscreteFunctionP0<Dimension, DataType>& gamma)
-{
-  using MeshType       = Mesh<Connectivity<Dimension>>;
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
-
-  DiscreteFunctionP0<Dimension, DataType> p{mesh};
-
-  parallel_for(
-    mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { p[cell_id] = ((gamma[cell_id] - 1) * rho[cell_id]) * epsilon[cell_id]; });
-
-  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(p);
-}
-
-std::shared_ptr<IDiscreteFunction>
-pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
-                     const std::shared_ptr<const IDiscreteFunction>& epsilon,
-                     const std::shared_ptr<const IDiscreteFunction>& gamma)
-{
-  std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, epsilon, gamma});
-  if (not mesh.use_count()) {
-    throw NormalError("rho, epsilon and gamma are not defined on the same mesh");
-  }
-
-  switch (mesh->dimension()) {
-  case 1: {
-    return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
-                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon),
-                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
-  }
-  case 2: {
-    return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
-                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon),
-                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
-  }
-  case 3: {
-    return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
-                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon),
-                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-template <size_t Dimension, typename DataType>
-std::shared_ptr<IDiscreteFunction>
-totalEnergyFromEpsilonU(const DiscreteFunctionP0<Dimension, DataType>& epsilon,
-                        const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u)
-{
-  using MeshType       = Mesh<Connectivity<Dimension>>;
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(epsilon.mesh());
-
-  DiscreteFunctionP0<Dimension, DataType> E{mesh};
-
-  parallel_for(
-    mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { E[cell_id] = epsilon[cell_id] + 0.5 * (u[cell_id], u[cell_id]); });
-
-  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
-}
-
-std::shared_ptr<IDiscreteFunction>
-totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon,
-                        const std::shared_ptr<const IDiscreteFunction>& u)
-{
-  std::shared_ptr<const IMesh> mesh = getCommonMesh({epsilon, u});
-  if (not mesh.use_count()) {
-    throw NormalError("epsilon and u are not defined on the same mesh");
-  }
-
-  switch (mesh->dimension()) {
-  case 1: {
-    return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon),
-                                   dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u));
-  }
-  case 2: {
-    return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon),
-                                   dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u));
-  }
-  case 3: {
-    return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon),
-                                   dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u));
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-template <size_t Dimension, typename DataType>
-std::shared_ptr<IDiscreteFunction>
-epsilonFromTotalEnergyU(const DiscreteFunctionP0<Dimension, DataType>& E,
-                        const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u)
-{
-  using MeshType       = Mesh<Connectivity<Dimension>>;
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(u.mesh());
-
-  DiscreteFunctionP0<Dimension, DataType> epsilon{mesh};
-
-  parallel_for(
-    mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = E[cell_id] - 0.5 * (u[cell_id], u[cell_id]); });
-
-  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
-}
-
-std::shared_ptr<IDiscreteFunction>
-epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E,
-                        const std::shared_ptr<const IDiscreteFunction>& u)
-{
-  std::shared_ptr<const IMesh> mesh = getCommonMesh({E, u});
-  if (not mesh.use_count()) {
-    throw NormalError("E and u are not defined on the same mesh");
-  }
-
-  switch (mesh->dimension()) {
-  case 1: {
-    return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*E),
-                                   dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u));
-  }
-  case 2: {
-    return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*E),
-                                   dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u));
-  }
-  case 3: {
-    return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*E),
-                                   dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u));
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-}   // namespace perfect_gas
diff --git a/src/scheme/PerfectGas.hpp b/src/scheme/PerfectGas.hpp
deleted file mode 100644
index fb7b9d212b4681f365faa04e4565dfd78c316f66..0000000000000000000000000000000000000000
--- a/src/scheme/PerfectGas.hpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef PERFECT_GAS_HPP
-#define PERFECT_GAS_HPP
-
-#include <scheme/IDiscreteFunction.hpp>
-
-namespace perfect_gas
-{
-std::shared_ptr<IDiscreteFunction> soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho,
-                                              const std::shared_ptr<const IDiscreteFunction>& p,
-                                              const std::shared_ptr<const IDiscreteFunction>& gamma);
-
-std::shared_ptr<IDiscreteFunction> epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
-                                                        const std::shared_ptr<const IDiscreteFunction>& p,
-                                                        const std::shared_ptr<const IDiscreteFunction>& gamma);
-
-std::shared_ptr<IDiscreteFunction> pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
-                                                        const std::shared_ptr<const IDiscreteFunction>& epsilon,
-                                                        const std::shared_ptr<const IDiscreteFunction>& gamma);
-
-// This is a temporary function
-// Todo: TO REMOVE
-std::shared_ptr<IDiscreteFunction> totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon,
-                                                           const std::shared_ptr<const IDiscreteFunction>& u);
-
-// This is a temporary function
-// Todo: TO REMOVE
-std::shared_ptr<IDiscreteFunction> epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E,
-                                                           const std::shared_ptr<const IDiscreteFunction>& u);
-
-}   // namespace perfect_gas
-
-#endif   // PERFECT_GAS_HPP
diff --git a/src/utils/Messenger.hpp b/src/utils/Messenger.hpp
index 99de952b4d9e084a7981e0e0f43fe5feb4f8b4ab..7e8b30c1a8e92c05fc30178e325110617c73c823 100644
--- a/src/utils/Messenger.hpp
+++ b/src/utils/Messenger.hpp
@@ -510,7 +510,6 @@ class Messenger
   allReduceSum(const DataType& data) const
   {
     static_assert(not std::is_const_v<DataType>);
-    static_assert(std::is_arithmetic_v<DataType>);
     static_assert(not std::is_same_v<DataType, bool>);
 
 #ifdef PUGS_HAS_MPI
diff --git a/src/utils/PugsUtils.hpp b/src/utils/PugsUtils.hpp
index 7429a75b3be20ec3d62028cc9180a2017f7b4f15..91267a2af92928084ce07e1021baf6336e8a8f98 100644
--- a/src/utils/PugsUtils.hpp
+++ b/src/utils/PugsUtils.hpp
@@ -4,7 +4,6 @@
 #include <Kokkos_Core.hpp>
 #include <utils/PugsMacros.hpp>
 
-#include <functional>
 #include <string>
 
 template <typename FunctionType>
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index ac6b547a82d532e0edf2e5c57d68c41ae97afeee..e1be8069885b80a3e26ede233dac3941b7525521 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -48,7 +48,7 @@ TEST_CASE("BiCGStab", "[algebra]")
 
     BiCGStab{A, x, b, 1e-12, 10, false};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
   }
 
   SECTION("no preconditionner non-converged")
@@ -91,6 +91,6 @@ TEST_CASE("BiCGStab", "[algebra]")
 
     BiCGStab{A, x, b, 1e-12, 1, true};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) > 1E-5 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) > 1E-5 * std::sqrt(dot(x, x)));
   }
 }
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index dda3becc1897f2e659ed105256bc0ec09271e46d..a1fb3fe82b982bc50d2b5247c1158941da24ce0d 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -278,6 +278,39 @@ let z:Z, z = round(-1.2);
       CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "z", int64_t{-1});
     }
 
+    {   // dot
+      tested_function_set.insert("dot:R^1*R^1");
+      std::string_view data = R"(
+import math;
+let x:R^1, x = -2;
+let y:R^1, y = 4;
+let s:R, s = dot(x,y);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", double{-2 * 4});
+    }
+
+    {   // dot
+      tested_function_set.insert("dot:R^2*R^2");
+      std::string_view data = R"(
+import math;
+let x:R^2, x = (-2, 3);
+let y:R^2, y = (4, 3);
+let s:R, s = dot(x,y);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", dot(TinyVector<2>{-2, 3}, TinyVector<2>{4, 3}));
+    }
+
+    {   // dot
+      tested_function_set.insert("dot:R^3*R^3");
+      std::string_view data = R"(
+import math;
+let x:R^3, x = (-2, 3, 4);
+let y:R^3, y = (4, 3, 5);
+let s:R, s = dot(x,y);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", dot(TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
+    }
+
     MathModule math_module;
 
     bool missing_test = false;
diff --git a/tests/test_CG.cpp b/tests/test_CG.cpp
index 00751efb3ccbdf5fe86da3a242c0b6dad92b097b..05919858f3ada681eb371883c6fcf6fd7d15a724 100644
--- a/tests/test_CG.cpp
+++ b/tests/test_CG.cpp
@@ -48,7 +48,7 @@ TEST_CASE("CG", "[algebra]")
 
     CG{A, x, b, 1e-12, 10, true};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
   }
 
   SECTION("singular matrix with RHS in its image")
@@ -64,7 +64,7 @@ TEST_CASE("CG", "[algebra]")
     x = 0;
 
     CG{A, x, b, 1e-12, 10};
-    REQUIRE(std::sqrt((x, x)) == 0);
+    REQUIRE(std::sqrt(dot(x, x)) == 0);
   }
 
   SECTION("no preconditionner non-converged")
@@ -107,6 +107,6 @@ TEST_CASE("CG", "[algebra]")
 
     CG{A, x, b, 1e-12, 1, false};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) > 1E-10 * std::sqrt(dot(x, x)));
   }
 }
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 9ced9e9424a7d1c80a41cd9dd5a807476adea9d7..5504ae390fd259eaa700a911b158bf6a6e46a625 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -194,7 +194,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
       }
 
@@ -221,7 +221,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("CG Diagonal")
@@ -235,7 +235,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("CG ICholeski")
@@ -249,7 +249,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("CG AMG")
@@ -263,7 +263,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -281,7 +281,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
 
 #else    // PUGS_HAS_PETSC
@@ -349,7 +349,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
       }
 
@@ -376,7 +376,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("BICGStab Diagonal")
@@ -390,7 +390,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("BICGStab ILU")
@@ -404,7 +404,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -426,7 +426,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("BICGStab2 Diagonal")
@@ -440,7 +440,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -462,7 +462,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("GMRES Diagonal")
@@ -476,7 +476,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("GMRES ILU")
@@ -490,7 +490,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -508,7 +508,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
 
 #else    // PUGS_HAS_PETSC
diff --git a/tests/test_MathModule.cpp b/tests/test_MathModule.cpp
index e899190bf486ff605c1bd4f8af2ccc232cc806a8..7c2a9bc5438367b7f3cf70d14d3b8acc58a29ae7 100644
--- a/tests/test_MathModule.cpp
+++ b/tests/test_MathModule.cpp
@@ -4,6 +4,8 @@
 #include <language/modules/MathModule.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 
+#include <set>
+
 // clazy:excludeall=non-pod-global-static
 
 TEST_CASE("MathModule", "[language]")
@@ -13,7 +15,7 @@ TEST_CASE("MathModule", "[language]")
   MathModule math_module;
   const auto& name_builtin_function = math_module.getNameBuiltinFunctionMap();
 
-  REQUIRE(name_builtin_function.size() == 22);
+  REQUIRE(name_builtin_function.size() == 25);
 
   SECTION("double -> double")
   {
@@ -322,4 +324,61 @@ TEST_CASE("MathModule", "[language]")
       REQUIRE(std::get<decltype(result)>(result_variant) == Catch::Approx(result));
     }
   }
+
+  SECTION("(R^d, R^d) -> double")
+  {
+    SECTION("dot:R^1*R^1")
+    {
+      TinyVector<1> arg0 = 3;
+      TinyVector<1> arg1 = 2;
+
+      DataVariant arg0_variant = arg0;
+      DataVariant arg1_variant = arg1;
+
+      auto i_function = name_builtin_function.find("dot:R^1*R^1");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg0_variant, arg1_variant});
+
+      const double result = dot(arg0, arg1);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+
+    SECTION("dot:R^2*R^2")
+    {
+      TinyVector<2> arg0{3, 2};
+      TinyVector<2> arg1{-2, 5};
+
+      DataVariant arg0_variant = arg0;
+      DataVariant arg1_variant = arg1;
+
+      auto i_function = name_builtin_function.find("dot:R^2*R^2");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg0_variant, arg1_variant});
+
+      const double result = dot(arg0, arg1);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+
+    SECTION("dot:R^3*R^3")
+    {
+      TinyVector<3> arg0{3, 2, 4};
+      TinyVector<3> arg1{-2, 5, 2};
+
+      DataVariant arg0_variant = arg0;
+      DataVariant arg1_variant = arg1;
+
+      auto i_function = name_builtin_function.find("dot:R^3*R^3");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg0_variant, arg1_variant});
+
+      const double result = dot(arg0, arg1);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+  }
 }
diff --git a/tests/test_TinyVector.cpp b/tests/test_TinyVector.cpp
index 2775588af687aa274dbba5374af53a85cacce36a..499536f9c702d5c3a1f2fcf8c4da7b1ecb461038 100644
--- a/tests/test_TinyVector.cpp
+++ b/tests/test_TinyVector.cpp
@@ -26,7 +26,7 @@ TEST_CASE("TinyVector", "[algebra]")
   REQUIRE(Catch::Detail::stringify(v) == "(3,2,4)");
 
   TinyVector<3, int> w(1, 2, 6);
-  REQUIRE((v, w) == 31);
+  REQUIRE(dot(v, w) == 31);
 
   w = 2 * v;
   REQUIRE(w == TinyVector<3, int>(6, 4, 8));
diff --git a/tests/test_Vector.cpp b/tests/test_Vector.cpp
index d02481859e1c41c49509f973f09006ba50373cc3..89fedcec868e026cfdf569eb8667b57b730a6e16 100644
--- a/tests/test_Vector.cpp
+++ b/tests/test_Vector.cpp
@@ -145,8 +145,8 @@ TEST_CASE("Vector", "[algebra]")
     y[3] = 2;
     y[4] = 8;
 
-    const int dot = (x, y);
-    REQUIRE(dot == 49);
+    const int s = dot(x, y);
+    REQUIRE(s == 49);
   }
 
   SECTION("self scalar division")