diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 4d8f15775bc356f8ce7d85ce8f3a27bd35c45e22..f7f1baeafe34261a38033dd3d2db828df334ce25 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -5,6 +5,7 @@ add_library(PugsLanguageModules
   BuiltinModule.cpp
   CoreModule.cpp
   LinearSolverModule.cpp
+  MathFunctionRegisterForVh.cpp
   MathModule.cpp
   MeshModule.cpp
   ModuleRepository.cpp
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..68b22b552f932a1fb48baeb96d6742a74a4717c2
--- /dev/null
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -0,0 +1,150 @@
+#include <language/modules/MathFunctionRegisterForVh.hpp>
+
+#include <language/modules/SchemeModule.hpp>
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+#include <language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module)
+{
+  scheme_module._addBuiltinFunction("sqrt",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return sqrt(a); }));
+
+  scheme_module._addBuiltinFunction("abs",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return abs(a); }));
+
+  scheme_module._addBuiltinFunction("sin",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return sin(a); }));
+
+  scheme_module._addBuiltinFunction("cos",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return cos(a); }));
+
+  scheme_module._addBuiltinFunction("tan",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return tan(a); }));
+
+  scheme_module._addBuiltinFunction("asin",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return asin(a); }));
+
+  scheme_module._addBuiltinFunction("acos",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return acos(a); }));
+
+  scheme_module._addBuiltinFunction("atan",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return atan(a); }));
+
+  scheme_module
+    ._addBuiltinFunction("atan2",
+                         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 atan2(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("atan2",
+                         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 atan2(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("atan2",
+                         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 atan2(a, b); }));
+
+  scheme_module._addBuiltinFunction("sinh",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return sinh(a); }));
+
+  scheme_module._addBuiltinFunction("cosh",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return cosh(a); }));
+
+  scheme_module._addBuiltinFunction("tanh",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return tanh(a); }));
+
+  scheme_module._addBuiltinFunction("asinh",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return asinh(a); }));
+
+  scheme_module._addBuiltinFunction("acosh",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return acosh(a); }));
+
+  scheme_module._addBuiltinFunction("atanh",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return atanh(a); }));
+
+  scheme_module._addBuiltinFunction("exp",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return exp(a); }));
+
+  scheme_module._addBuiltinFunction("log",
+                                    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
+                                      std::shared_ptr<const IDiscreteFunction>)>>(
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return log(a); }));
+
+  scheme_module
+    ._addBuiltinFunction("pow",
+                         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 pow(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("pow",
+                         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 pow(a, b); }));
+
+  scheme_module
+    ._addBuiltinFunction("pow",
+                         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 pow(a, b); }));
+}
diff --git a/src/language/modules/MathFunctionRegisterForVh.hpp b/src/language/modules/MathFunctionRegisterForVh.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..95ee461bc8e999fc75e273c31c0aea4505edda0e
--- /dev/null
+++ b/src/language/modules/MathFunctionRegisterForVh.hpp
@@ -0,0 +1,12 @@
+#ifndef MATH_FUNCTION_REGISTER_FOR_VH_HPP
+#define MATH_FUNCTION_REGISTER_FOR_VH_HPP
+
+class SchemeModule;
+
+class MathFunctionRegisterForVh
+{
+ public:
+  MathFunctionRegisterForVh(SchemeModule& module);
+};
+
+#endif   // MATH_FUNCTION_REGISTER_FOR_VH_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index db1988abac81bb0613e473bfa7429022d0f8f90d..898f2a71c15b07e8bbc8946795f2450a9f315957 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,6 +1,7 @@
 #include <language/modules/SchemeModule.hpp>
 
 #include <language/modules/BinaryOperatorRegisterForVh.hpp>
+#include <language/modules/MathFunctionRegisterForVh.hpp>
 #include <language/modules/UnaryOperatorRegisterForVh.hpp>
 #include <language/utils/BinaryOperatorProcessorBuilder.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
@@ -296,6 +297,8 @@ SchemeModule::SchemeModule()
                               [](const std::shared_ptr<const IDiscreteFunction>& c) -> double { return acoustic_dt(c); }
 
                               ));
+
+  MathFunctionRegisterForVh{*this};
 }
 
 void
diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp
index 8e5ee95d8f6797e705cba0436829a565c7714e6c..758b5f9c56fcdc928adfc4678a0a02a8dcfa5bd8 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -27,6 +27,8 @@ inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IDiscreteFu
 
 class SchemeModule : public BuiltinModule
 {
+  friend class MathFunctionRegisterForVh;
+
  public:
   std::string_view
   name() const final
diff --git a/src/language/utils/CMakeLists.txt b/src/language/utils/CMakeLists.txt
index bed3da329d5534e8aa4f48ef48cc0751dfa00f3b..216a6be6a15971fe93ca2f4341f720e8b6e5cde3 100644
--- a/src/language/utils/CMakeLists.txt
+++ b/src/language/utils/CMakeLists.txt
@@ -24,6 +24,7 @@ add_library(PugsLanguageUtils
   DataVariant.cpp
   EmbeddedData.cpp
   EmbeddedIDiscreteFunctionOperators.cpp
+  EmbeddedIDiscreteFunctionMathFunctions.cpp
   FunctionSymbolId.cpp
   IncDecOperatorRegisterForN.cpp
   IncDecOperatorRegisterForR.cpp
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2939e011af950fc414cde730e2c1206bf8f4ae01
--- /dev/null
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -0,0 +1,307 @@
+#include <language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp>
+
+#include <language/utils/EmbeddedIDiscreteFunctionUtils.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+#define DISCRETE_FUNCTION_CALL(FUNCTION, ARG)                                                                         \
+  if (ARG->dataType() == ASTNodeDataType::double_t and ARG->descriptor().type() == DiscreteFunctionType::P0) {        \
+    switch (f->mesh()->dimension()) {                                                                                 \
+    case 1: {                                                                                                         \
+      using DiscreteFunctionType = DiscreteFunctionP0<1, double>;                                                     \
+      return std::make_shared<const DiscreteFunctionType>(FUNCTION(dynamic_cast<const DiscreteFunctionType&>(*ARG))); \
+    }                                                                                                                 \
+    case 2: {                                                                                                         \
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;                                                     \
+      return std::make_shared<const DiscreteFunctionType>(FUNCTION(dynamic_cast<const DiscreteFunctionType&>(*ARG))); \
+    }                                                                                                                 \
+    case 3: {                                                                                                         \
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;                                                     \
+      return std::make_shared<const DiscreteFunctionType>(FUNCTION(dynamic_cast<const DiscreteFunctionType&>(*ARG))); \
+    }                                                                                                                 \
+    default: {                                                                                                        \
+      throw UnexpectedError("invalid mesh dimension");                                                                \
+    }                                                                                                                 \
+    }                                                                                                                 \
+  } else {                                                                                                            \
+    throw UnexpectedError("invalid operand type " + operand_type_name(ARG));                                          \
+  }
+
+std::shared_ptr<const IDiscreteFunction>
+sqrt(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(sqrt, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+abs(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(abs, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+sin(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(sin, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+cos(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(cos, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+tan(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(tan, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+asin(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(asin, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+acos(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(acos, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+atan(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(atan, f);
+}
+
+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)) {
+    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>(
+        atan2(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        atan2(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        atan2(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>
+atan2(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>(atan2(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(atan2(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(atan2(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+atan2(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>(atan2(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(atan2(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(atan2(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+sinh(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(sinh, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+cosh(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(cosh, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+tanh(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(tanh, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+asinh(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(asinh, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+acosh(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(acosh, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+atanh(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(atanh, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+exp(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(exp, f);
+}
+
+std::shared_ptr<const IDiscreteFunction>
+log(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  DISCRETE_FUNCTION_CALL(log, f);
+}
+
+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)) {
+    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>(
+        pow(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        pow(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(
+        pow(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>
+pow(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>(pow(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(pow(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(pow(a, dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+pow(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>(pow(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 2: {
+      using DiscreteFunctionType = DiscreteFunctionP0<2, double>;
+      return std::make_shared<const DiscreteFunctionType>(pow(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    case 3: {
+      using DiscreteFunctionType = DiscreteFunctionP0<3, double>;
+      return std::make_shared<const DiscreteFunctionType>(pow(dynamic_cast<const DiscreteFunctionType&>(*f), a));
+    }
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+    }
+  } else {
+    throw UnexpectedError("invalid operand type " + operand_type_name(f));
+  }
+}
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0030f46759dfe683a09228131140e98ec205fc8b
--- /dev/null
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -0,0 +1,57 @@
+#ifndef EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP
+#define EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <memory>
+
+class IDiscreteFunction;
+
+std::shared_ptr<const IDiscreteFunction> sqrt(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> abs(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> sin(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> cos(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> tan(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> asin(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> acos(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> atan(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> atan2(const std::shared_ptr<const IDiscreteFunction>&,
+                                               const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> atan2(const double, const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> atan2(const std::shared_ptr<const IDiscreteFunction>&, const double);
+
+std::shared_ptr<const IDiscreteFunction> sinh(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> cosh(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> tanh(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> asinh(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> acosh(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> atanh(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> exp(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> log(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> pow(const std::shared_ptr<const IDiscreteFunction>&,
+                                             const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> pow(const double, const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> pow(const std::shared_ptr<const IDiscreteFunction>&, const double);
+
+#endif   // EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
index a9691615ba86788e978bdc7696ba705afd20d236..815436bb0d90d8ded91b66ee4e5f3d0e3345c558 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
@@ -2,58 +2,13 @@
 
 #include <language/node_processor/BinaryExpressionProcessor.hpp>
 #include <language/node_processor/UnaryExpressionProcessor.hpp>
+#include <language/utils/EmbeddedIDiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename T>
-PUGS_INLINE std::string
-operand_type_name(const T& t)
-{
-  if constexpr (is_shared_ptr_v<T>) {
-    Assert(t.use_count() > 0);
-    return operand_type_name(*t);
-  } else if constexpr (std::is_base_of_v<IDiscreteFunction, std::decay_t<T>>) {
-    return "Vh(" + name(t.descriptor().type()) + ':' + dataTypeName(t.dataType()) + ')';
-  } else {
-    return dataTypeName(ast_node_data_type_from<T>);
-  }
-}
-
-PUGS_INLINE
-bool
-isSameDiscretization(const IDiscreteFunction& f, const IDiscreteFunction& g)
-{
-  if ((f.dataType() == g.dataType()) and (f.descriptor().type() == g.descriptor().type())) {
-    switch (f.dataType()) {
-    case ASTNodeDataType::double_t: {
-      return true;
-    }
-    case ASTNodeDataType::vector_t: {
-      return f.dataType().dimension() == g.dataType().dimension();
-    }
-    case ASTNodeDataType::matrix_t: {
-      return (f.dataType().nbRows() == g.dataType().nbRows()) and
-             (f.dataType().nbColumns() == g.dataType().nbColumns());
-    }
-    default: {
-      throw UnexpectedError("invalid data type " + operand_type_name(f));
-    }
-    }
-  } else {
-    return false;
-  }
-}
-
-PUGS_INLINE
-bool
-isSameDiscretization(const std::shared_ptr<const IDiscreteFunction>& f,
-                     const std::shared_ptr<const IDiscreteFunction>& g)
-{
-  return isSameDiscretization(*f, *g);
-}
-
 template <typename LHS_T, typename RHS_T>
 PUGS_INLINE std::string
 invalid_operands(const LHS_T& f, const RHS_T& g)
@@ -286,14 +241,16 @@ std::shared_ptr<const IDiscreteFunction>
 innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
                     const std::shared_ptr<const IDiscreteFunction>& g)
 {
-  if (f->mesh() != g->mesh()) {
-    throw NormalError("discrete functions defined on different meshes");
+  std::shared_ptr mesh = getCommonMesh({f, g});
+  if (mesh.use_count() == 0) {
+    throw NormalError("operands are defined on different meshes");
   }
+
   if (not isSameDiscretization(f, g)) {
     throw NormalError(invalid_operands(f, g));
   }
 
-  switch (f->mesh()->dimension()) {
+  switch (mesh->dimension()) {
   case 1: {
     return innerCompositionLaw<BinOperatorT, 1>(f, g);
   }
@@ -473,13 +430,14 @@ std::shared_ptr<const IDiscreteFunction>
 applyBinaryOperation(const std::shared_ptr<const IDiscreteFunction>& f,
                      const std::shared_ptr<const IDiscreteFunction>& g)
 {
-  if (f->mesh() != g->mesh()) {
-    throw NormalError("functions defined on different meshes");
+  std::shared_ptr mesh = getCommonMesh({f, g});
+  if (mesh.use_count() == 0) {
+    throw NormalError("operands are defined on different meshes");
   }
 
   Assert(not isSameDiscretization(f, g), "should call inner composition instead");
 
-  switch (f->mesh()->dimension()) {
+  switch (mesh->dimension()) {
   case 1: {
     return applyBinaryOperation<BinOperatorT, 1>(f, g);
   }
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp b/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9321c3273d7955e77934f9f4233ae3603cf14990
--- /dev/null
+++ b/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
@@ -0,0 +1,58 @@
+#ifndef EMBEDDED_I_DISCRETE_FUNCTION_UTILS_HPP
+#define EMBEDDED_I_DISCRETE_FUNCTION_UTILS_HPP
+
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <sstream>
+#include <string>
+
+template <typename T>
+PUGS_INLINE std::string
+operand_type_name(const T& t)
+{
+  if constexpr (is_shared_ptr_v<T>) {
+    Assert(t.use_count() > 0);
+    return operand_type_name(*t);
+  } else if constexpr (std::is_base_of_v<IDiscreteFunction, std::decay_t<T>>) {
+    return "Vh(" + name(t.descriptor().type()) + ':' + dataTypeName(t.dataType()) + ')';
+  } else {
+    return dataTypeName(ast_node_data_type_from<T>);
+  }
+}
+
+PUGS_INLINE
+bool
+isSameDiscretization(const IDiscreteFunction& f, const IDiscreteFunction& g)
+{
+  if ((f.dataType() == g.dataType()) and (f.descriptor().type() == g.descriptor().type())) {
+    switch (f.dataType()) {
+    case ASTNodeDataType::double_t: {
+      return true;
+    }
+    case ASTNodeDataType::vector_t: {
+      return f.dataType().dimension() == g.dataType().dimension();
+    }
+    case ASTNodeDataType::matrix_t: {
+      return (f.dataType().nbRows() == g.dataType().nbRows()) and
+             (f.dataType().nbColumns() == g.dataType().nbColumns());
+    }
+    default: {
+      throw UnexpectedError("invalid data type " + operand_type_name(f));
+    }
+    }
+  } else {
+    return false;
+  }
+}
+
+PUGS_INLINE
+bool
+isSameDiscretization(const std::shared_ptr<const IDiscreteFunction>& f,
+                     const std::shared_ptr<const IDiscreteFunction>& g)
+{
+  return isSameDiscretization(*f, *g);
+}
+
+#endif   // EMBEDDED_I_DISCRETE_FUNCTION_UTILS_HPP
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 2c90637fbf9069d4e940d8d7edf20302a55f7613..0b77aaea481451c6662dd450756affc02ce0226d 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -114,6 +114,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator+(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
     Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
     DiscreteFunctionP0<Dimension, decltype(DataType{} + DataType2T{})> sum(f.m_mesh);
     parallel_for(
@@ -126,6 +127,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator+(const LHSDataType& a, const DiscreteFunctionP0& g)
   {
     using SumDataType = decltype(LHSDataType{} + DataType{});
+    Assert(g.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, SumDataType> sum(g.m_mesh);
     parallel_for(
       g.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = a + g[cell_id]; });
@@ -137,6 +139,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator+(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
     using SumDataType = decltype(DataType{} + RHSDataType{});
+    Assert(f.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, SumDataType> sum(f.m_mesh);
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = f[cell_id] + b; });
@@ -148,6 +151,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator-(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
     Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
     DiscreteFunctionP0<Dimension, decltype(DataType{} - DataType2T{})> difference(f.m_mesh);
     parallel_for(
@@ -160,6 +164,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator-(const LHSDataType& a, const DiscreteFunctionP0& g)
   {
     using DifferenceDataType = decltype(LHSDataType{} - DataType{});
+    Assert(g.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, DifferenceDataType> difference(g.m_mesh);
     parallel_for(
       g.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = a - g[cell_id]; });
@@ -171,6 +176,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator-(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
     using DifferenceDataType = decltype(DataType{} - RHSDataType{});
+    Assert(f.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, DifferenceDataType> difference(f.m_mesh);
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = f[cell_id] - b; });
@@ -182,6 +188,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator*(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
+    Assert(f.m_cell_values.isBuilt());
     Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
     DiscreteFunctionP0<Dimension, decltype(DataType{} * DataType2T{})> product(f.m_mesh);
     parallel_for(
@@ -194,6 +201,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator*(const LHSDataType& a, const DiscreteFunctionP0& f)
   {
     using ProductDataType = decltype(LHSDataType{} * DataType{});
+    Assert(f.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, ProductDataType> product(f.m_mesh);
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = a * f[cell_id]; });
@@ -205,6 +213,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator*(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
     using ProductDataType = decltype(DataType{} * RHSDataType{});
+    Assert(f.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, ProductDataType> product(f.m_mesh);
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = f[cell_id] * b; });
@@ -216,6 +225,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator/(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
+    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
     Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
     DiscreteFunctionP0<Dimension, decltype(DataType{} / DataType2T{})> ratio(f.m_mesh);
     parallel_for(
@@ -228,6 +238,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator/(const LHSDataType& a, const DiscreteFunctionP0& f)
   {
     using RatioDataType = decltype(LHSDataType{} / DataType{});
+    Assert(f.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, RatioDataType> ratio(f.m_mesh);
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = a / f[cell_id]; });
@@ -239,12 +250,279 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator/(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
     using RatioDataType = decltype(DataType{} / RHSDataType{});
+    Assert(f.m_cell_values.isBuilt());
     DiscreteFunctionP0<Dimension, RatioDataType> ratio(f.m_mesh);
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = f[cell_id] / b; });
     return ratio;
   }
 
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  sqrt(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::sqrt(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  abs(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::abs(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  sin(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::sin(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  cos(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::cos(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  tan(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::tan(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  asin(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::asin(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  acos(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::acos(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  atan(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::atan(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  sinh(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::sinh(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  cosh(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::cosh(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  tanh(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::tanh(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  asinh(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::asinh(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  acosh(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::acosh(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  atanh(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::atanh(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  exp(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::exp(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  log(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::log(f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  atan2(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::atan2(f[cell_id], g[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  atan2(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::atan2(a, f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  atan2(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::atan2(f[cell_id], a); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  pow(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::pow(f[cell_id], g[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  pow(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::pow(a, f[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  pow(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::pow(f[cell_id], a); });
+
+    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)