diff --git a/doc/userdoc.org b/doc/userdoc.org
index 01a88b0af5800e6565b3dc01d99bdfb4de569aba..ba7cd798a6a7991d6451c405a5b9c81dfa2c0362 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -2587,6 +2587,61 @@ $\mathbb{R}^1$, $\mathbb{R}^2$ and $\mathbb{R}^3$.
 The output is
 #+RESULTS: dot-examples
 
+A set of ~det~ functions is defined to get the determinant of matrices
+of $\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and
+$\mathbb{R}^{3\times3}$.
+#+NAME: det-examples
+#+BEGIN_SRC pugs :exports both :results output
+   import math;
+   cout << "det([[1.2]]) = " << det([[1.2]]) << "\n";
+   cout << "det([[1,2],[3,4]]) = " << det([[1,2],[3,4]]) << "\n";
+   cout << "det([[1,2,3],[4,5,6],[7,8,9]]) = "
+        << det([[1,2,3],[4,5,6],[7,8,9]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: det-examples
+
+The ~trace~ functions compute the trace of matrices of
+$\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$.
+#+NAME: trace-examples
+#+BEGIN_SRC pugs :exports both :results output
+   import math;
+   cout << "trace([[1.2]]) = " << trace([[1.2]]) << "\n";
+   cout << "trace([[1,2],[3,4]]) = " << trace([[1,2],[3,4]]) << "\n";
+   cout << "trace([[1,2,3],[4,5,6],[7,8,9]]) = "
+        << trace([[1,2,3],[4,5,6],[7,8,9]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: trace-examples
+
+Also, one can compute inverses of $\mathbb{R}^{1\times1}$,
+$\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ matrices using the
+~inverse~ function set.
+#+NAME: inverse-examples
+#+BEGIN_SRC pugs :exports both :results output
+   import math;
+   cout << "inverse([[1.2]]) = " << inverse([[1.2]]) << "\n";
+   cout << "inverse([[1,2],[3,4]]) = " << inverse([[1,2],[3,4]]) << "\n";
+   cout << "inverse([[3,2,1],[5,6,4],[7,8,9]]) = "
+        << inverse([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: inverse-examples
+
+Transpose of matrices of $\mathbb{R}^{1\times1}$,
+$\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ are obtained using the
+~transpose~ functions.
+#+NAME: transpose-examples
+#+BEGIN_SRC pugs :exports both :results output
+   import math;
+   cout << "transpose([[1.2]]) = " << transpose([[1.2]]) << "\n";
+   cout << "transpose([[1,2],[3,4]]) = " << transpose([[1,2],[3,4]]) << "\n";
+   cout << "transpose([[3,2,1],[5,6,4],[7,8,9]]) = "
+        << transpose([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: transpose-examples
+
 #+BEGIN_note
 Observe that the use of a proper rounding or truncation function is
 the right way to convert a real value to an integer one. Available
@@ -3195,11 +3250,12 @@ description.
 
 ****** ~Vh -> Vh~
 
-These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
-return value is also a $\mathbb{P}_0(\mathbb{R})$ function. The value
-is simply the application of the function to the cell values.
+The majority of these functions are defined for
+$\mathbb{P}_0(\mathbb{R})$ data and the return value is also a
+$\mathbb{P}_0(\mathbb{R})$ function. The value is simply the
+application of the function to the cell values.
 
-Here is the list of the functions
+Here is the list of these functions
 - ~abs: Vh -> Vh~
 - ~acos: Vh -> Vh~
 - ~acosh: Vh -> Vh~
@@ -3217,6 +3273,21 @@ Here is the list of the functions
 - ~tan: Vh -> Vh~
 - ~tanh: Vh -> Vh~
 
+A few functions are defined for $\mathbb{P}_0(\mathbb{R^{}}^{1\times1})$,
+$\mathbb{P}_0(\mathbb{R^{}}^{2\times2})$ and
+$\mathbb{P}_0(\mathbb{R^{}}^{3\times3})$ data and the return value is a
+$\mathbb{P}_0(\mathbb{R})$ function. The value is simply the
+application of the function to the cell values.
+- ~det: Vh -> Vh~
+- ~trace: Vh -> Vh~
+
+Also, functions are defined for $\mathbb{P}_0(\mathbb{R}^{d\times d})$
+data and the return value is a $\mathbb{P}_0(\mathbb{R}^{d\times d})$
+function, with $d\in\{1,2,3\}$. The value is simply the application of
+the function to the cell values.
+- ~inverse: Vh -> Vh~
+- ~transpose: Vh -> Vh~
+
 ******  ~Vh*Vh -> Vh~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index e34f98ed8c25df7a77ba90bbd7cd05305b5eb115..1717b48c78deaedc7cabd9f60feb215e29a80ebe 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -21,16 +21,21 @@ class [[nodiscard]] TinyMatrix
   using data_type = T;
 
  private:
+  static_assert((M > 0), "TinyMatrix number of rows must be strictly positive");
+  static_assert((N > 0), "TinyMatrix number of columns must be strictly positive");
+
   T m_values[M * N];
 
   PUGS_FORCEINLINE
-  constexpr size_t _index(size_t i, size_t j) const noexcept   // LCOV_EXCL_LINE (due to forced inline)
+  constexpr size_t
+  _index(size_t i, size_t j) const noexcept   // LCOV_EXCL_LINE (due to forced inline)
   {
     return i * N + j;
   }
 
   template <typename... Args>
-  PUGS_FORCEINLINE constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
+  PUGS_FORCEINLINE constexpr void
+  _unpackVariadicInput(const T& t, Args&&... args) noexcept
   {
     m_values[M * N - 1 - sizeof...(args)] = t;
     if constexpr (sizeof...(args) > 0) {
@@ -40,13 +45,15 @@ class [[nodiscard]] TinyMatrix
 
  public:
   PUGS_INLINE
-  constexpr bool isSquare() const noexcept
+  constexpr bool
+  isSquare() const noexcept
   {
     return M == N;
   }
 
   PUGS_INLINE
-  constexpr friend TinyMatrix<N, M, T> transpose(const TinyMatrix& A)
+  constexpr friend TinyMatrix<N, M, T>
+  transpose(const TinyMatrix& A)
   {
     TinyMatrix<N, M, T> tA;
     for (size_t i = 0; i < M; ++i) {
@@ -58,31 +65,36 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr size_t dimension() const
+  constexpr size_t
+  dimension() const
   {
     return M * N;
   }
 
   PUGS_INLINE
-  constexpr size_t numberOfValues() const
+  constexpr size_t
+  numberOfValues() const
   {
     return this->dimension();
   }
 
   PUGS_INLINE
-  constexpr size_t numberOfRows() const
+  constexpr size_t
+  numberOfRows() const
   {
     return M;
   }
 
   PUGS_INLINE
-  constexpr size_t numberOfColumns() const
+  constexpr size_t
+  numberOfColumns() const
   {
     return N;
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix operator-() const
+  constexpr TinyMatrix
+  operator-() const
   {
     TinyMatrix opposite;
     for (size_t i = 0; i < M * N; ++i) {
@@ -92,20 +104,23 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr friend TinyMatrix operator*(const T& t, const TinyMatrix& A)
+  constexpr friend TinyMatrix
+  operator*(const T& t, const TinyMatrix& A)
   {
     TinyMatrix B = A;
     return B *= t;
   }
 
   PUGS_INLINE
-  constexpr friend TinyMatrix operator*(const T& t, TinyMatrix&& A)
+  constexpr friend TinyMatrix
+  operator*(const T& t, TinyMatrix&& A)
   {
     return std::move(A *= t);
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix& operator*=(const T& t)
+  constexpr TinyMatrix&
+  operator*=(const T& t)
   {
     for (size_t i = 0; i < M * N; ++i) {
       m_values[i] *= t;
@@ -114,7 +129,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   template <size_t P>
-  PUGS_INLINE constexpr TinyMatrix<M, P, T> operator*(const TinyMatrix<N, P, T>& B) const
+  PUGS_INLINE constexpr TinyMatrix<M, P, T>
+  operator*(const TinyMatrix<N, P, T>& B) const
   {
     const TinyMatrix& A = *this;
     TinyMatrix<M, P, T> AB;
@@ -131,7 +147,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr TinyVector<M, T> operator*(const TinyVector<N, T>& x) const
+  constexpr TinyVector<M, T>
+  operator*(const TinyVector<N, T>& x) const
   {
     const TinyMatrix& A = *this;
     TinyVector<M, T> Ax;
@@ -146,7 +163,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
+  constexpr friend std::ostream&
+  operator<<(std::ostream& os, const TinyMatrix& A)
   {
     os << '[';
     for (size_t i = 0; i < M; ++i) {
@@ -165,7 +183,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr bool operator==(const TinyMatrix& A) const
+  constexpr bool
+  operator==(const TinyMatrix& A) const
   {
     for (size_t i = 0; i < M * N; ++i) {
       if (m_values[i] != A.m_values[i])
@@ -175,13 +194,15 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr bool operator!=(const TinyMatrix& A) const
+  constexpr bool
+  operator!=(const TinyMatrix& A) const
   {
     return not this->operator==(A);
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix operator+(const TinyMatrix& A) const
+  constexpr TinyMatrix
+  operator+(const TinyMatrix& A) const
   {
     TinyMatrix sum;
     for (size_t i = 0; i < M * N; ++i) {
@@ -191,14 +212,16 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix operator+(TinyMatrix&& A) const
+  constexpr TinyMatrix
+  operator+(TinyMatrix&& A) const
   {
     A += *this;
     return std::move(A);
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix operator-(const TinyMatrix& A) const
+  constexpr TinyMatrix
+  operator-(const TinyMatrix& A) const
   {
     TinyMatrix difference;
     for (size_t i = 0; i < M * N; ++i) {
@@ -208,7 +231,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix operator-(TinyMatrix&& A) const
+  constexpr TinyMatrix
+  operator-(TinyMatrix&& A) const
   {
     for (size_t i = 0; i < M * N; ++i) {
       A.m_values[i] = m_values[i] - A.m_values[i];
@@ -217,7 +241,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix& operator+=(const TinyMatrix& A)
+  constexpr TinyMatrix&
+  operator+=(const TinyMatrix& A)
   {
     for (size_t i = 0; i < M * N; ++i) {
       m_values[i] += A.m_values[i];
@@ -226,7 +251,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr void operator+=(const volatile TinyMatrix& A) volatile
+  constexpr void
+  operator+=(const volatile TinyMatrix& A) volatile
   {
     for (size_t i = 0; i < M * N; ++i) {
       m_values[i] += A.m_values[i];
@@ -234,7 +260,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix& operator-=(const TinyMatrix& A)
+  constexpr TinyMatrix&
+  operator-=(const TinyMatrix& A)
   {
     for (size_t i = 0; i < M * N; ++i) {
       m_values[i] -= A.m_values[i];
@@ -243,21 +270,24 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr T& operator()(size_t i, size_t j) noexcept(NO_ASSERT)
+  constexpr T&
+  operator()(size_t i, size_t j) noexcept(NO_ASSERT)
   {
     Assert((i < M) and (j < N));
     return m_values[_index(i, j)];
   }
 
   PUGS_INLINE
-  constexpr const T& operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
+  constexpr const T&
+  operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
   {
     Assert((i < M) and (j < N));
     return m_values[_index(i, j)];
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix& operator=(ZeroType) noexcept
+  constexpr TinyMatrix&
+  operator=(ZeroType) noexcept
   {
     static_assert(std::is_arithmetic<T>(), "Cannot assign 'zero' value for non-arithmetic types");
     for (size_t i = 0; i < M * N; ++i) {
@@ -267,7 +297,8 @@ class [[nodiscard]] TinyMatrix
   }
 
   PUGS_INLINE
-  constexpr TinyMatrix& operator=(IdentityType) noexcept
+  constexpr TinyMatrix&
+  operator=(IdentityType) noexcept
   {
     static_assert(std::is_arithmetic<T>(), "Cannot assign 'identity' value for non-arithmetic types");
     for (size_t i = 0; i < M; ++i) {
@@ -329,7 +360,7 @@ class [[nodiscard]] TinyMatrix
   constexpr TinyMatrix(const TinyMatrix&) noexcept = default;
 
   PUGS_INLINE
-  TinyMatrix(TinyMatrix && A) noexcept = default;
+  TinyMatrix(TinyMatrix&& A) noexcept = default;
 
   PUGS_INLINE
   ~TinyMatrix() = default;
@@ -450,6 +481,19 @@ getMinor(const TinyMatrix<M, N, T>& A, size_t I, size_t J)
   return m;
 }
 
+template <size_t N, typename T>
+PUGS_INLINE T
+trace(const TinyMatrix<N, N, T>& A)
+{
+  static_assert(std::is_arithmetic<T>::value, "trace is not defined for non-arithmetic types");
+
+  T t = A(0, 0);
+  for (size_t i = 1; i < N; ++i) {
+    t += A(i, i);
+  }
+  return t;
+}
+
 template <size_t N, typename T>
 PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
 
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index fcf3b31a73da2c47bb59e27bc169195255effb30..05d2ea703cfbce9d233fc285f2b38ebbb805ae14 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -209,6 +209,35 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
 
                                              ));
 
+  scheme_module._addBuiltinFunction("det", std::function(
+
+                                             [](std::shared_ptr<const IDiscreteFunction> A)
+                                               -> std::shared_ptr<const IDiscreteFunction> { return det(A); }
+
+                                             ));
+
+  scheme_module._addBuiltinFunction("trace", std::function(
+
+                                               [](std::shared_ptr<const IDiscreteFunction> A)
+                                                 -> std::shared_ptr<const IDiscreteFunction> { return trace(A); }
+
+                                               ));
+
+  scheme_module._addBuiltinFunction("inverse", std::function(
+
+                                                 [](std::shared_ptr<const IDiscreteFunction> A)
+                                                   -> std::shared_ptr<const IDiscreteFunction> { return inverse(A); }
+
+                                                 ));
+
+  scheme_module._addBuiltinFunction("transpose",
+                                    std::function(
+
+                                      [](std::shared_ptr<const IDiscreteFunction> A)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return transpose(A); }
+
+                                      ));
+
   scheme_module._addBuiltinFunction("min", std::function(
 
                                              [](std::shared_ptr<const IDiscreteFunction> a) -> double { return min(a); }
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index 4e8e10d1534f55f8c02cccb2f3f48974fdc944d2..4a1e7b35b76a4d9b8b05ebb1355f3cdca4f751d9 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -70,6 +70,36 @@ MathModule::MathModule()
   this->_addBuiltinFunction("dot", std::function([](const TinyVector<3>& x, const TinyVector<3>& y) -> double {
                               return dot(x, y);
                             }));
+
+  this->_addBuiltinFunction("det", std::function([](const TinyMatrix<1> A) -> double { return det(A); }));
+
+  this->_addBuiltinFunction("det", std::function([](const TinyMatrix<2>& A) -> double { return det(A); }));
+
+  this->_addBuiltinFunction("det", std::function([](const TinyMatrix<3>& A) -> double { return det(A); }));
+
+  this->_addBuiltinFunction("trace", std::function([](const TinyMatrix<1> A) -> double { return trace(A); }));
+
+  this->_addBuiltinFunction("trace", std::function([](const TinyMatrix<2>& A) -> double { return trace(A); }));
+
+  this->_addBuiltinFunction("trace", std::function([](const TinyMatrix<3>& A) -> double { return trace(A); }));
+
+  this->_addBuiltinFunction("inverse",
+                            std::function([](const TinyMatrix<1> A) -> TinyMatrix<1> { return inverse(A); }));
+
+  this->_addBuiltinFunction("inverse",
+                            std::function([](const TinyMatrix<2>& A) -> TinyMatrix<2> { return inverse(A); }));
+
+  this->_addBuiltinFunction("inverse",
+                            std::function([](const TinyMatrix<3>& A) -> TinyMatrix<3> { return inverse(A); }));
+
+  this->_addBuiltinFunction("transpose",
+                            std::function([](const TinyMatrix<1> A) -> TinyMatrix<1> { return transpose(A); }));
+
+  this->_addBuiltinFunction("transpose",
+                            std::function([](const TinyMatrix<2>& A) -> TinyMatrix<2> { return transpose(A); }));
+
+  this->_addBuiltinFunction("transpose",
+                            std::function([](const TinyMatrix<3>& A) -> TinyMatrix<3> { return transpose(A); }));
 }
 
 void
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index f83bd418592c29b5cded8b77bda6c717ff3af804..67622f4a7c2ebbc7284f7305be98aeb3beb8768f 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -503,6 +503,210 @@ template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<2>&,
 template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<3>&,
                                                       const std::shared_ptr<const IDiscreteFunction>&);
 
+template <typename MatrixType>
+std::shared_ptr<const IDiscreteFunction>
+det(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  switch (A->mesh()->dimension()) {
+  case 1: {
+    using DiscreteFunctionType = DiscreteFunctionP0<1, MatrixType>;
+    return std::make_shared<const DiscreteFunctionP0<1, double>>(det(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 2: {
+    using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+    return std::make_shared<const DiscreteFunctionP0<2, double>>(det(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 3: {
+    using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+    return std::make_shared<const DiscreteFunctionP0<3, double>>(det(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+det(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  if (A->dataType() == ASTNodeDataType::matrix_t and A->descriptor().type() == DiscreteFunctionType::P0) {
+    if (A->dataType().numberOfRows() != A->dataType().numberOfColumns()) {
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+    }
+    switch (A->dataType().numberOfRows()) {
+    case 1: {
+      return det<TinyMatrix<1>>(A);
+    }
+    case 2: {
+      return det<TinyMatrix<2>>(A);
+    }
+    case 3: {
+      return det<TinyMatrix<3>>(A);
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("invalid matrix type");
+    }
+      // LCOV_EXCL_STOP
+    }
+  } else {
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+  }
+}
+
+template <typename MatrixType>
+std::shared_ptr<const IDiscreteFunction>
+trace(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  switch (A->mesh()->dimension()) {
+  case 1: {
+    using DiscreteFunctionType = DiscreteFunctionP0<1, MatrixType>;
+    return std::make_shared<const DiscreteFunctionP0<1, double>>(trace(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 2: {
+    using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+    return std::make_shared<const DiscreteFunctionP0<2, double>>(trace(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 3: {
+    using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+    return std::make_shared<const DiscreteFunctionP0<3, double>>(trace(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+trace(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  if (A->dataType() == ASTNodeDataType::matrix_t and A->descriptor().type() == DiscreteFunctionType::P0) {
+    if (A->dataType().numberOfRows() != A->dataType().numberOfColumns()) {
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+    }
+    switch (A->dataType().numberOfRows()) {
+    case 1: {
+      return trace<TinyMatrix<1>>(A);
+    }
+    case 2: {
+      return trace<TinyMatrix<2>>(A);
+    }
+    case 3: {
+      return trace<TinyMatrix<3>>(A);
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("invalid matrix type");
+    }
+      // LCOV_EXCL_STOP
+    }
+  } else {
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+  }
+}
+
+template <typename MatrixType>
+std::shared_ptr<const IDiscreteFunction>
+inverse(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  switch (A->mesh()->dimension()) {
+  case 1: {
+    using DiscreteFunctionType = DiscreteFunctionP0<1, MatrixType>;
+    return std::make_shared<const DiscreteFunctionType>(inverse(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 2: {
+    using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+    return std::make_shared<const DiscreteFunctionType>(inverse(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 3: {
+    using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+    return std::make_shared<const DiscreteFunctionType>(inverse(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+inverse(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  if (A->dataType() == ASTNodeDataType::matrix_t and A->descriptor().type() == DiscreteFunctionType::P0) {
+    if (A->dataType().numberOfRows() != A->dataType().numberOfColumns()) {
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+    }
+    switch (A->dataType().numberOfRows()) {
+    case 1: {
+      return inverse<TinyMatrix<1>>(A);
+    }
+    case 2: {
+      return inverse<TinyMatrix<2>>(A);
+    }
+    case 3: {
+      return inverse<TinyMatrix<3>>(A);
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("invalid matrix type");
+    }
+      // LCOV_EXCL_STOP
+    }
+  } else {
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+  }
+}
+
+template <typename MatrixType>
+std::shared_ptr<const IDiscreteFunction>
+transpose(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  switch (A->mesh()->dimension()) {
+  case 1: {
+    using DiscreteFunctionType = DiscreteFunctionP0<1, MatrixType>;
+    return std::make_shared<const DiscreteFunctionType>(transpose(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 2: {
+    using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+    return std::make_shared<const DiscreteFunctionType>(transpose(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  case 3: {
+    using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+    return std::make_shared<const DiscreteFunctionType>(transpose(dynamic_cast<const DiscreteFunctionType&>(*A)));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+transpose(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+  if (A->dataType() == ASTNodeDataType::matrix_t and A->descriptor().type() == DiscreteFunctionType::P0) {
+    if (A->dataType().numberOfRows() != A->dataType().numberOfColumns()) {
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+    }
+    switch (A->dataType().numberOfRows()) {
+    case 1: {
+      return transpose<TinyMatrix<1>>(A);
+    }
+    case 2: {
+      return transpose<TinyMatrix<2>>(A);
+    }
+    case 3: {
+      return transpose<TinyMatrix<3>>(A);
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("invalid matrix type");
+    }
+      // LCOV_EXCL_STOP
+    }
+  } else {
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+  }
+}
+
 double
 min(const std::shared_ptr<const IDiscreteFunction>& f)
 {
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
index 2369429a6267a8fdbc803d75b4213703dd970a48..24994c743428d98befb8d34cc18c2cbbf12a2ba4 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -65,6 +65,14 @@ template <size_t VectorDimension>
 std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&,
                                              const std::shared_ptr<const IDiscreteFunction>&);
 
+std::shared_ptr<const IDiscreteFunction> det(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> trace(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> inverse(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> transpose(const std::shared_ptr<const IDiscreteFunction>&);
+
 std::shared_ptr<const IDiscreteFunction> sum_of_Vh_components(const std::shared_ptr<const IDiscreteFunction>&);
 
 std::shared_ptr<const IDiscreteFunction> vectorize(
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index bb6b484a6bc10bd803e16b8bb54a7cba25a349de..d1bd72c3394346db10898cd83f1c609974f17697 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -524,11 +524,62 @@ class DiscreteFunctionP0 : public IDiscreteFunction
     return result;
   }
 
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  det(const DiscreteFunctionP0& A)
+  {
+    static_assert(is_tiny_matrix_v<DataType>);
+    Assert(A.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, double> result{A.m_mesh};
+    parallel_for(
+      A.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = det(A[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  trace(const DiscreteFunctionP0& A)
+  {
+    static_assert(is_tiny_matrix_v<DataType>);
+    Assert(A.m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, double> result{A.m_mesh};
+    parallel_for(
+      A.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = trace(A[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  inverse(const DiscreteFunctionP0& A)
+  {
+    static_assert(is_tiny_matrix_v<DataType>);
+    static_assert(DataType::NumberOfRows == DataType::NumberOfColumns, "cannot compute inverse of non square matrices");
+    Assert(A.m_cell_values.isBuilt());
+    DiscreteFunctionP0 result{A.m_mesh};
+    parallel_for(
+      A.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = inverse(A[cell_id]); });
+
+    return result;
+  }
+
+  PUGS_INLINE friend DiscreteFunctionP0
+  transpose(const DiscreteFunctionP0& A)
+  {
+    static_assert(is_tiny_matrix_v<DataType>);
+    static_assert(DataType::NumberOfRows == DataType::NumberOfColumns, "cannot compute inverse of non square matrices");
+    Assert(A.m_cell_values.isBuilt());
+    DiscreteFunctionP0 result{A.m_mesh};
+    parallel_for(
+      A.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = transpose(A[cell_id]); });
+
+    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());
+    Assert(f.m_mesh == g.m_mesh);
     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]); });
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index 60a4bd491f392b7cc0c303d94531edcf8634edbc..8de9ec2c2ed9e61479aa2d7165bf6e506fd6caf6 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -358,6 +358,158 @@ 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}));
     }
 
+    {   // det
+      tested_function_set.insert("det:R^1x1");
+      std::string_view data = R"(
+import math;
+let A:R^1x1, A = -2;
+let d:R, d = det(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "d", det(TinyMatrix<1>{-2}));
+    }
+
+    {   // det
+      tested_function_set.insert("det:R^2x2");
+      std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+                  [ 3.2, 2.3]];
+let d:R, d = det(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "d",
+                                               det(TinyMatrix<2>{-2.5, 3.2,   //
+                                                                 +3.2, 2.3}));
+    }
+
+    {   // det
+      tested_function_set.insert("det:R^3x3");
+      std::string_view data = R"(
+import math;
+let A:R^3x3, A = [[-2, 2,-1],
+                  [ 3, 2, 2],
+                  [-2, 5,-3]];
+let d:R, d = det(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "d",
+                                               det(TinyMatrix<3>{-2, 2, -1,   //
+                                                                 +3, 2, +2,   //
+                                                                 -2, 5, -3}));
+    }
+
+    {   // trace
+      tested_function_set.insert("trace:R^1x1");
+      std::string_view data = R"(
+import math;
+let A:R^1x1, A = -2;
+let t:R, t = trace(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "t", trace(TinyMatrix<1>{-2}));
+    }
+
+    {   // trace
+      tested_function_set.insert("trace:R^2x2");
+      std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+                  [ 3.2, 2.3]];
+let t:R, t = trace(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "t",
+                                               trace(TinyMatrix<2>{-2.5, 3.2,   //
+                                                                   +3.2, 2.3}));
+    }
+
+    {   // trace
+      tested_function_set.insert("trace:R^3x3");
+      std::string_view data = R"(
+import math;
+let A:R^3x3, A = [[-2.5, 2.9,-1.3],
+                  [ 3.2, 2.3, 2.7],
+                  [-2.6, 5.2,-3.5]];
+let t:R, t = trace(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "t",
+                                               trace(TinyMatrix<3>{-2.5, 2.9, -1.3,   //
+                                                                   +3.2, 2.3, +2.7,   //
+                                                                   -2.6, 5.2, -3.5}));
+    }
+
+    {   // inverse
+      tested_function_set.insert("inverse:R^1x1");
+      std::string_view data = R"(
+import math;
+let A:R^1x1, A = -2;
+let invA:R^1x1, invA = inverse(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "invA", inverse(TinyMatrix<1>{-2}));
+    }
+
+    {   // inverse
+      tested_function_set.insert("inverse:R^2x2");
+      std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+                  [ 3.2, 2.3]];
+let invA:R^2x2, invA = inverse(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "invA",
+                                               inverse(TinyMatrix<2>{-2.5, 3.2,   //
+                                                                     +3.2, 2.3}));
+    }
+
+    {   // inverse
+      tested_function_set.insert("inverse:R^3x3");
+      std::string_view data = R"(
+import math;
+let A:R^3x3, A = [[-2, 2,-1],
+                  [ 3, 2, 2],
+                  [-2, 5,-3]];
+let invA:R^3x3, invA = inverse(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "invA",
+                                               inverse(TinyMatrix<3>{-2, 2, -1,   //
+                                                                     +3, 2, +2,   //
+                                                                     -2, 5, -3}));
+    }
+
+    {   // transpose
+      tested_function_set.insert("transpose:R^1x1");
+      std::string_view data = R"(
+import math;
+let A:R^1x1, A = -2;
+let tA:R^1x1, tA = transpose(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "tA", transpose(TinyMatrix<1>{-2}));
+    }
+
+    {   // transpose
+      tested_function_set.insert("transpose:R^2x2");
+      std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+                  [ 3.2, 2.3]];
+let tA:R^2x2, tA = transpose(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "tA",
+                                               transpose(TinyMatrix<2>{-2.5, 3.2,   //
+                                                                       +3.2, 2.3}));
+    }
+
+    {   // transpose
+      tested_function_set.insert("transpose:R^3x3");
+      std::string_view data = R"(
+import math;
+let A:R^3x3, A = [[-2.5, 2.9,-1.3],
+                  [ 3.2, 2.3, 2.7],
+                  [-2.6, 5.2,-3.5]];
+let tA:R^3x3, tA = transpose(A);
+)";
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "tA",
+                                               transpose(TinyMatrix<3>{-2.5, 2.9, -1.3,   //
+                                                                       +3.2, 2.3, +2.7,   //
+                                                                       -2.6, 5.2, -3.5}));
+    }
+
     MathModule math_module;
 
     bool missing_test = false;
diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index 9d83de5368f81f1dc46570234c2a7ed19ebc2da4..48c107e4ad2f07391a49270a04c9a3dc20d2db4a 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -3428,6 +3428,94 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
             CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, dot);
           }
 
+          SECTION("det(Ah)")
+          {
+            DiscreteFunctionP0<Dimension, TinyMatrix<2>> Ah{mesh};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                const double x = xj[cell_id][0];
+                Ah[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3,   //
+                                               -0.2 * x - 1, 2 + x};
+              });
+
+            {
+              DiscreteFunctionP0 result = det(Ah);
+              bool is_same              = true;
+              parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+                if (result[cell_id] != det(Ah[cell_id])) {
+                  is_same = false;
+                }
+              });
+              REQUIRE(is_same);
+            }
+          }
+
+          SECTION("trace(Ah)")
+          {
+            DiscreteFunctionP0<Dimension, TinyMatrix<2>> Ah{mesh};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                const double x = xj[cell_id][0];
+                Ah[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3,   //
+                                               -0.2 * x - 1, 2 + x};
+              });
+
+            {
+              DiscreteFunctionP0 result = trace(Ah);
+              bool is_same              = true;
+              parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+                if (result[cell_id] != trace(Ah[cell_id])) {
+                  is_same = false;
+                }
+              });
+              REQUIRE(is_same);
+            }
+          }
+
+          SECTION("inverse(Ah)")
+          {
+            DiscreteFunctionP0<Dimension, TinyMatrix<2>> Ah{mesh};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                const double x = xj[cell_id][0];
+                Ah[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3,   //
+                                               -0.2 * x - 1, 2 + x};
+              });
+
+            {
+              DiscreteFunctionP0 result = inverse(Ah);
+              bool is_same              = true;
+              parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+                if (result[cell_id] != inverse(Ah[cell_id])) {
+                  is_same = false;
+                }
+              });
+              REQUIRE(is_same);
+            }
+          }
+
+          SECTION("transpose(Ah)")
+          {
+            DiscreteFunctionP0<Dimension, TinyMatrix<2>> Ah{mesh};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                const double x = xj[cell_id][0];
+                Ah[cell_id]    = TinyMatrix<2>{x + 1, 2 * x - 3,   //
+                                               -0.2 * x - 1, 2 + x};
+              });
+
+            {
+              DiscreteFunctionP0 result = transpose(Ah);
+              bool is_same              = true;
+              parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+                if (result[cell_id] != transpose(Ah[cell_id])) {
+                  is_same = false;
+                }
+              });
+              REQUIRE(is_same);
+            }
+          }
+
           SECTION("scalar sum")
           {
             const CellValue<const double> cell_value = positive_function.cellValues();
diff --git a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
index 9508afc6987a8bf67a2a090892ba1bee211cd176..6d00782f96379548ebb66927044782def35612e8 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -538,6 +538,302 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
         }
 
+        SECTION("det Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = det(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = det(p_R2x2_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = det(p_R3x3_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(det(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(det(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(det(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("trace Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = trace(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = trace(p_R2x2_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = trace(p_R3x3_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(trace(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(trace(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(trace(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(trace(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("inverse Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = inverse(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id](0, 0) != inverse(values[cell_id])(0, 0)) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = inverse(p_R2x2_u);
+            constexpr size_t nb_row = 2;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = inverse(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = inverse(p_R3x3_u);
+            constexpr size_t nb_row = 3;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = inverse(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(inverse(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(inverse(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(inverse(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("transpose Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = transpose(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id](0, 0) != transpose(values[cell_id])(0, 0)) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = transpose(p_R2x2_u);
+            constexpr size_t nb_row = 2;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = transpose(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = transpose(p_R3x3_u);
+            constexpr size_t nb_row = 3;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = transpose(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(transpose(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(transpose(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(transpose(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(transpose(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
         SECTION("sum_of_Vh Vh -> Vh")
         {
           {
@@ -1110,6 +1406,302 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
         }
 
+        SECTION("det Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = det(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = det(p_R2x2_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = det(p_R3x3_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(det(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(det(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(det(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("trace Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = trace(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = trace(p_R2x2_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = trace(p_R3x3_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(trace(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(trace(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(trace(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(trace(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("inverse Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = inverse(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id](0, 0) != inverse(values[cell_id])(0, 0)) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = inverse(p_R2x2_u);
+            constexpr size_t nb_row = 2;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = inverse(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = inverse(p_R3x3_u);
+            constexpr size_t nb_row = 3;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = inverse(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(inverse(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(inverse(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(inverse(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("transpose Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = transpose(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id](0, 0) != transpose(values[cell_id])(0, 0)) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = transpose(p_R2x2_u);
+            constexpr size_t nb_row = 2;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = transpose(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = transpose(p_R3x3_u);
+            constexpr size_t nb_row = 3;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              const TinyMatrix invA          = transpose(values[cell_id]);
+              const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+              for (size_t i = 0; i < nb_row; ++i) {
+                for (size_t j = 0; j < nb_row; ++j) {
+                  if (fu_i(i, j) != invA(i, j)) {
+                    is_same = false;
+                  }
+                }
+              }
+              if (not is_same) {
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(transpose(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(transpose(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(transpose(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(transpose(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
         SECTION("sum_of_Vh Vh -> Vh")
         {
           {
@@ -1680,6 +2272,270 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
         }
 
+        SECTION("det Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = det(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = det(p_R2x2_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = det(p_R3x3_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != det(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(det(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(det(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(det(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("trace Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = trace(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = trace(p_R2x2_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu = trace(p_R3x3_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != trace(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(trace(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(trace(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(trace(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(trace(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("inverse Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = inverse(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != inverse(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = inverse(p_R2x2_u);
+            constexpr size_t nb_row = 2;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != inverse(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = inverse(p_R3x3_u);
+            constexpr size_t nb_row = 3;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != inverse(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(inverse(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(inverse(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(inverse(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("transpose Vh -> Vh")
+        {
+          {
+            std::shared_ptr p_fu = transpose(p_R1x1_u);
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+            auto values  = p_R1x1_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != transpose(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = transpose(p_R2x2_u);
+            constexpr size_t nb_row = 2;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R2x2_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != transpose(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          {
+            std::shared_ptr p_fu    = transpose(p_R3x3_u);
+            constexpr size_t nb_row = 3;
+
+            REQUIRE(p_fu.use_count() > 0);
+            REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+            const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+            auto values  = p_R3x3_u->cellValues();
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+              if (fu[cell_id] != transpose(values[cell_id])) {
+                is_same = false;
+                break;
+              }
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(transpose(p_u), "error: invalid operand type Vh(P0:R)");
+          REQUIRE_THROWS_WITH(transpose(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+          REQUIRE_THROWS_WITH(transpose(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+          REQUIRE_THROWS_WITH(transpose(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
         SECTION("sum_of_Vh Vh -> Vh")
         {
           {
diff --git a/tests/test_MathModule.cpp b/tests/test_MathModule.cpp
index 99966d4bacebb938fce5bc686661b5c53e8d708e..2cb3481b1e6cde7561733d8ea03decd46153db9d 100644
--- a/tests/test_MathModule.cpp
+++ b/tests/test_MathModule.cpp
@@ -13,7 +13,7 @@ TEST_CASE("MathModule", "[language]")
   MathModule math_module;
   const auto& name_builtin_function = math_module.getNameBuiltinFunctionMap();
 
-  REQUIRE(name_builtin_function.size() == 30);
+  REQUIRE(name_builtin_function.size() == 42);
 
   SECTION("Z -> N")
   {
@@ -457,4 +457,210 @@ TEST_CASE("MathModule", "[language]")
       REQUIRE(std::get<double>(result_variant) == result);
     }
   }
+
+  SECTION("R^dxd -> double")
+  {
+    SECTION("det:R^1x1")
+    {
+      TinyMatrix<1> arg{3};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("det:R^1x1");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const double result = det(arg);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+
+    SECTION("det:R^2x2")
+    {
+      TinyMatrix<2> arg{3, 2, 1, -2};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("det:R^2x2");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const double result = det(arg);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+
+    SECTION("det:R^3x3")
+    {
+      TinyMatrix<3> arg{3,  2,  4,   //
+                        -1, 4,  1,   //
+                        -4, -1, 5};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("det:R^3x3");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const double result = det(arg);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+
+    SECTION("trace:R^1x1")
+    {
+      TinyMatrix<1> arg{3};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("trace:R^1x1");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const double result = trace(arg);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+
+    SECTION("trace:R^2x2")
+    {
+      TinyMatrix<2> arg{3, 2, 1, -2};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("trace:R^2x2");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const double result = trace(arg);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+
+    SECTION("trace:R^3x3")
+    {
+      TinyMatrix<3> arg{3,  2,  4,   //
+                        -1, 4,  1,   //
+                        -4, -1, 5};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("trace:R^3x3");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const double result = trace(arg);
+      REQUIRE(std::get<double>(result_variant) == result);
+    }
+  }
+
+  SECTION("R^dxd -> R^dxd")
+  {
+    SECTION("inverse:R^1x1")
+    {
+      TinyMatrix<1> arg{3};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("inverse:R^1x1");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const TinyMatrix<1> result = inverse(arg);
+      REQUIRE(std::get<TinyMatrix<1>>(result_variant) == result);
+    }
+
+    SECTION("inverse:R^2x2")
+    {
+      TinyMatrix<2> arg{3, 2, 1, -2};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("inverse:R^2x2");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const TinyMatrix<2> result = inverse(arg);
+      REQUIRE(std::get<TinyMatrix<2>>(result_variant) == result);
+    }
+
+    SECTION("inverse:R^3x3")
+    {
+      TinyMatrix<3> arg{3,  2,  4,   //
+                        -1, 4,  1,   //
+                        -4, -1, 5};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("inverse:R^3x3");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const TinyMatrix<3> result = inverse(arg);
+      REQUIRE(std::get<TinyMatrix<3>>(result_variant) == result);
+    }
+
+    SECTION("transpose:R^1x1")
+    {
+      TinyMatrix<1> arg{3};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("transpose:R^1x1");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const TinyMatrix<1> result = transpose(arg);
+      REQUIRE(std::get<TinyMatrix<1>>(result_variant) == result);
+    }
+
+    SECTION("transpose:R^2x2")
+    {
+      TinyMatrix<2> arg{3, 2, 1, -2};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("transpose:R^2x2");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const TinyMatrix<2> result = transpose(arg);
+      REQUIRE(std::get<TinyMatrix<2>>(result_variant) == result);
+    }
+
+    SECTION("transpose:R^3x3")
+    {
+      TinyMatrix<3> arg{3,  2,  4,   //
+                        -1, 4,  1,   //
+                        -4, -1, 5};
+
+      DataVariant arg_variant = arg;
+
+      auto i_function = name_builtin_function.find("transpose:R^3x3");
+      REQUIRE(i_function != name_builtin_function.end());
+
+      IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+      DataVariant result_variant                  = function_embedder.apply({arg_variant});
+
+      const TinyMatrix<3> result = transpose(arg);
+      REQUIRE(std::get<TinyMatrix<3>>(result_variant) == result);
+    }
+  }
 }
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 2f4c09d4ad9ce1dbb5e3d856d8737d27344459f2..f6f09cf55dd2ad1ffef55c9275b989c203598f7e 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -194,6 +194,17 @@ TEST_CASE("TinyMatrix", "[algebra]")
     REQUIRE(det(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
   }
 
+  SECTION("checking for trace calculations")
+  {
+    REQUIRE(trace(TinyMatrix<1, 1, int>(6)) == 6);
+    REQUIRE(trace(TinyMatrix<2, 2, int>(5, 1, -3, 6)) == 5 + 6);
+    REQUIRE(trace(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1 + 2 + 3);
+    REQUIRE(trace(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == 6 + 34 + 7);
+    REQUIRE(trace(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
+            Catch::Approx(1 + 4 + 2 + 17.5).epsilon(1E-14));
+    REQUIRE(trace(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 1 + 0 + 1 + 2);
+  }
+
   SECTION("checking for inverse calculations")
   {
     {