diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 30691b1da34424f349ef9c47a70848f21ebccc75..fc57d41ae78776aa6c977ae6f0895e1cd9fc0388 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -132,6 +132,16 @@ class [[nodiscard]] TinyMatrix
     return *this;
   }
 
+  [[nodiscard]] PUGS_INLINE constexpr friend T
+  doubleDot(const TinyMatrix& A, const TinyMatrix& B)
+  {
+    T t = A.m_values[0] * B.m_values[0];
+    for (size_t i = 1; i < M * N; ++i) {
+      t += A.m_values[i] * B.m_values[i];
+    }
+    return t;
+  }
+
   template <size_t P>
   [[nodiscard]] PUGS_INLINE constexpr TinyMatrix<M, P, T>
   operator*(const TinyMatrix<N, P, T>& B) const noexcept
@@ -201,18 +211,6 @@ class [[nodiscard]] TinyMatrix
     return not this->operator==(A);
   }
 
-  [[nodiscard]] PUGS_INLINE constexpr friend T
-  dot(const TinyMatrix& A, const TinyMatrix& B) noexcept
-  {
-    T sum = A.m_values[0] * B.m_values[0];
-
-    for (size_t i = 1; i < M * N; ++i) {
-      sum += A.m_values[i] * B.m_values[i];
-    }
-
-    return sum;
-  }
-
   [[nodiscard]] PUGS_INLINE constexpr TinyMatrix
   operator+(const TinyMatrix& A) const noexcept
   {
@@ -269,24 +267,21 @@ class [[nodiscard]] TinyMatrix
     return *this;
   }
 
-  PUGS_INLINE
-  constexpr T&
+  [[nodiscard]] PUGS_INLINE 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&
+  [[nodiscard]] PUGS_INLINE 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&
+  PUGS_INLINE constexpr TinyMatrix&
   operator=(ZeroType) noexcept
   {
     static_assert(std::is_arithmetic<T>(), "Cannot assign 'zero' value for non-arithmetic types");
@@ -386,16 +381,9 @@ tensorProduct(const TinyVector<M, T>& x, const TinyVector<N, T>& y) noexcept
   return A;
 }
 
-template <size_t M, size_t N, typename T>
-[[nodiscard]] PUGS_INLINE constexpr auto
-frobeniusNorm(const TinyMatrix<M, N, T>& A) noexcept
-{
-  static_assert(std::is_arithmetic<T>(), "Cannot compute frobeniusNorm value for non-arithmetic types");
-  return std::sqrt(dot(A, A));
-}
-
 template <size_t N, typename T>
 [[nodiscard]] PUGS_INLINE constexpr T
+
 det(const TinyMatrix<N, N, T>& A) noexcept
 {
   static_assert(std::is_arithmetic<T>::value, "determinant is not defined for non-arithmetic types");
@@ -509,23 +497,17 @@ trace(const TinyMatrix<N, N, T>& A) noexcept
   return t;
 }
 
-template <size_t N, typename T>
-PUGS_INLINE T
-frobeniusNorm(const TinyMatrix<N, N, T>& A)
+template <size_t M, size_t N, typename T>
+[[nodiscard]] PUGS_INLINE constexpr decltype(std::sqrt(std::declval<T>()))
+frobeniusNorm(const TinyMatrix<M, N, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "norm is not defined for non-arithmetic types");
+  static_assert(std::is_floating_point<T>::value, "Frobenius norm is defined for floating point types only");
 
-  T t = 0;
-  for (size_t i = 0; i < N; ++i) {
-    for (size_t j = 0; j < N; ++j) {
-      t += A(i, j) * A(i, j);
-    }
-  }
-  return std::sqrt(t);
+  return std::sqrt(doubleDot(A, A));
 }
 
 template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
 [[nodiscard]] PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
 
 template <typename T>
@@ -579,73 +561,4 @@ inverse(const TinyMatrix<3, 3, T>& A) noexcept(NO_ASSERT)
   return A_cofactors_T *= 1. / determinent;
 }
 
-template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N, N, T>
-inverse(const TinyMatrix<N, N, T>& A)
-{
-  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
-  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
-  TinyVector<N, size_t> indexc(zero);
-  TinyVector<N, size_t> indexr(zero);
-  TinyVector<N, size_t> indpiv(zero);
-  TinyMatrix<N> inv_A(A);
-  double big    = 0;
-  double pivinv = 0;
-  double dum    = 0;
-  size_t irow   = 0;
-  size_t icol   = 0;
-
-  for (size_t i = 0; i < N; ++i) {
-    big = 0;
-    for (size_t j = 0; j < N; ++j) {
-      if (indpiv[j] != 1) {
-        for (size_t k = 0; k < N; k++) {
-          if (indpiv[k] == 0) {
-            if (std::abs(inv_A(j, k)) >= big) {
-              big  = std::abs(inv_A(j, k));
-              irow = j;
-              icol = k;
-            }
-          } else {
-            Assert(indpiv[k] <= 1, "Gauss pivoting: singular matrix");
-          }
-        }
-      }
-    }
-    ++indpiv[icol];
-    // We found the pivot A(irow,icol), now we swap columns to put the pivot on the diagonal
-    if (irow != icol) {
-      for (size_t l = 0; l < N; ++l) {
-        std::swap(inv_A(irow, l), inv_A(icol, l));
-      }
-    }
-    // we save the swap to invert it at the end
-    indexr[i] = irow;
-    indexc[i] = icol;
-    Assert(inv_A(icol, icol) != 0, "Gauss pivoting: singular matrix");
-    pivinv            = 1 / inv_A(icol, icol);
-    inv_A(icol, icol) = 1;
-    for (size_t l = 0; l < N; ++l) {
-      inv_A(icol, l) *= pivinv;
-    }
-    for (size_t ll = 0; ll < N; ++ll) {
-      if (ll != icol) {
-        dum             = inv_A(ll, icol);
-        inv_A(ll, icol) = 0;
-        for (size_t l = 0; l < N; ++l) {
-          inv_A(ll, l) -= inv_A(icol, l) * dum;
-        }
-      }
-    }
-  }
-  for (size_t l = N; l > 0; --l) {
-    if (indexr[l - 1] != indexc[l - 1]) {
-      for (size_t k = 0; k < N; ++k) {
-        std::swap(inv_A(k, indexr[l - 1]), inv_A(k, indexc[l - 1]));
-      }
-    }
-  }
-  return inv_A;
-}
-
 #endif   // TINYMATRIX_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b4adfde7a9c29785856ae2a8415639459bd0e1ed..e179415db7a0bb10fa5d182a2f8d4ba9279b804a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -37,6 +37,7 @@
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
 #include <scheme/HyperelasticSolver.hpp>
 #include <scheme/HyperplasticSolver.hpp>
+#include <scheme/HyperplasticSolverO2.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/InflowBoundaryConditionDescriptor.hpp>
@@ -1009,6 +1010,375 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("hyperplastic_eucclhyd_fluxes_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list)
+                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
+                                return HyperplasticSolverO2Handler{getCommonMesh({rho, aL, aT, u, sigma})}
+                                  .solver()
+                                  .compute_fluxes(HyperplasticSolverO2Handler::SolverType::Eucclhyd, rho, aL, aT, u,
+                                                  sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_eucclhyd_solver_implicit_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{
+                                  getCommonMesh({rho, u, E, Fe, zeta, yield, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(HyperplasticSolverO2Handler::SolverType::Eucclhyd,
+                                         HyperplasticSolverO2Handler::RelaxationType::Implicit, dt, rho, u, E, Fe, zeta,
+                                         yield, aL, aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_eucclhyd_solver_exponential_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{
+                                  getCommonMesh({rho, u, E, Fe, zeta, yield, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(HyperplasticSolverO2Handler::SolverType::Eucclhyd,
+                                         HyperplasticSolverO2Handler::RelaxationType::Exponential, dt, rho, u, E, Fe,
+                                         zeta, yield, aL, aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_eucclhyd_solver_cauchygreen_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{
+                                  getCommonMesh({rho, u, E, Fe, zeta, yield, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(HyperplasticSolverO2Handler::SolverType::Eucclhyd,
+                                         HyperplasticSolverO2Handler::RelaxationType::CauchyGreen, dt, rho, u, E, Fe,
+                                         zeta, yield, aL, aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_plastic_step_implicit_o2",
+                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& sigman,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& sigmanp1,
+                                             const double& dt) -> std::shared_ptr<const DiscreteFunctionVariant> {
+                              return HyperplasticSolverO2Handler{getCommonMesh({Fe, zeta, yield, sigman, sigmanp1})}
+                                .solver()
+                                .apply_plastic_relaxation(HyperplasticSolverO2Handler::RelaxationType::Implicit, dt, Fe,
+                                                          zeta, yield, sigman, sigmanp1);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("hyperplastic_plastic_step_exponential_o2",
+                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& sigman,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& sigmanp1,
+                                             const double& dt) -> std::shared_ptr<const DiscreteFunctionVariant> {
+                              return HyperplasticSolverO2Handler{getCommonMesh({Fe, zeta, yield, sigman, sigmanp1})}
+                                .solver()
+                                .apply_plastic_relaxation(HyperplasticSolverO2Handler::RelaxationType::Exponential, dt,
+                                                          Fe, zeta, yield, sigman, sigmanp1);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("hyperplastic_elastic_step_eucclhyd_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe, aL, aT, sigma})}
+                                  .solver()
+                                  .apply_elastic(HyperplasticSolverO2Handler::SolverType::Eucclhyd, dt, rho, u, E, Fe,
+                                                 aL, aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_glace_fluxes_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list)
+                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
+                                return HyperplasticSolverO2Handler{getCommonMesh({rho, aL, aT, u, sigma})}
+                                  .solver()
+                                  .compute_fluxes(HyperplasticSolverO2Handler::SolverType::Glace,   //
+                                                  rho, aL, aT, u, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_elastic_step_glace_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe, aL, aT, sigma})}
+                                  .solver()
+                                  .apply_elastic(HyperplasticSolverO2Handler::SolverType::Glace, dt, rho, u, E, Fe, aL,
+                                                 aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_glace_solver_implicit_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{
+                                  getCommonMesh({rho, u, E, Fe, zeta, yield, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(HyperplasticSolverO2Handler::SolverType::Glace,
+                                         HyperplasticSolverO2Handler::RelaxationType::Implicit, dt, rho, u, E, Fe, zeta,
+                                         yield, aL, aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_glace_solver_exponential_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{
+                                  getCommonMesh({rho, u, E, Fe, zeta, yield, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(HyperplasticSolverO2Handler::SolverType::Glace,
+                                         HyperplasticSolverO2Handler::RelaxationType::Exponential, dt, rho, u, E, Fe,
+                                         zeta, yield, aL, aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperplastic_glace_solver_cauchygreen_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{
+                                  getCommonMesh({rho, u, E, Fe, zeta, yield, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(HyperplasticSolverO2Handler::SolverType::Glace,
+                                         HyperplasticSolverO2Handler::RelaxationType::CauchyGreen, dt, rho, u, E, Fe,
+                                         zeta, yield, aL, aT, sigma, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("apply_hyperplastic_fluxes_implicit_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,      //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,       //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,     //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,    //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,    //
+                                 const std::shared_ptr<const ItemValueVariant>& ur,              //
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})}   //
+                                  .solver()
+                                  .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr,
+                                                HyperplasticSolverO2Handler::RelaxationType::Implicit);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("apply_hyperplastic_fluxes_cauchygreen_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,      //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,       //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,     //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,    //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,    //
+                                 const std::shared_ptr<const ItemValueVariant>& ur,              //
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})}   //
+                                  .solver()
+                                  .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr,
+                                                HyperplasticSolverO2Handler::RelaxationType::CauchyGreen);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("apply_hyperplastic_fluxes_exponential_o2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,      //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,       //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,     //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,    //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,    //
+                                 const std::shared_ptr<const ItemValueVariant>& ur,              //
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})}   //
+                                  .solver()
+                                  .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr,
+                                                HyperplasticSolverO2Handler::RelaxationType::Exponential);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("lagrangian",
                             std::function(
 
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 812504f72000fef674a7938f707d707fd1d982b4..931047e36b957471f41232cb69228b0bf7b2480a 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -5,6 +5,7 @@ add_library(
   AcousticSolver.cpp
   HyperelasticSolver.cpp
   HyperplasticSolver.cpp
+  HyperplasticSolverO2.cpp
   DiscreteFunctionIntegrator.cpp
   DiscreteFunctionInterpoler.cpp
   DiscreteFunctionUtils.cpp
diff --git a/src/scheme/HyperplasticSolverO2.cpp b/src/scheme/HyperplasticSolverO2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c99d0100cb902f8d274696acaa2ff179bb996565
--- /dev/null
+++ b/src/scheme/HyperplasticSolverO2.cpp
@@ -0,0 +1,1707 @@
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
+#include <algebra/EigenvalueSolver.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <analysis/Polynomial.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/StencilArray.hpp>
+#include <mesh/StencilManager.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/HyperplasticSolverO2.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
+#include <variant>
+#include <vector>
+
+template <MeshConcept MeshType>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
+  : public HyperplasticSolverO2Handler::IHyperplasticSolverO2
+{
+ private:
+  static constexpr size_t Dimension = MeshType::Dimension;
+  using Rdxd                        = TinyMatrix<Dimension>;
+  using R3x3                        = TinyMatrix<3>;
+  using Rd                          = TinyVector<Dimension>;
+  using R3                          = TinyVector<3>;
+
+  using MeshDataType = MeshData<MeshType>;
+
+  using DiscreteScalarFunction   = DiscreteFunctionP0<const double>;
+  using DiscreteVectorFunction   = DiscreteFunctionP0<const Rd>;
+  using DiscreteTensorFunction   = DiscreteFunctionP0<const Rdxd>;
+  using DiscreteTensorFunction3d = DiscreteFunctionP0<const R3x3>;
+
+  class FixedBoundaryCondition;
+  class PressureBoundaryCondition;
+  class NormalStressBoundaryCondition;
+  class SymmetryBoundaryCondition;
+  class VelocityBoundaryCondition;
+
+  using BoundaryCondition = std::variant<FixedBoundaryCondition,
+                                         PressureBoundaryCondition,
+                                         NormalStressBoundaryCondition,
+                                         SymmetryBoundaryCondition,
+                                         VelocityBoundaryCondition>;
+
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+  R3x3
+  to3D(const Rdxd tensor) const
+  {
+    R3x3 tensor3d = zero;
+    for (size_t i = 0; i < Dimension; i++) {
+      for (size_t j = 0; j < Dimension; j++) {
+        tensor3d(i, j) = tensor(i, j);
+      }
+    }
+    return tensor3d;
+  }
+
+  R3
+  to3D(const Rd vector) const
+  {
+    R3 vector3d = zero;
+    for (size_t i = 0; i < Dimension; i++) {
+      vector3d[i] = vector[i];
+    }
+    return vector3d;
+  }
+
+  CellValue<const R3x3>
+  to3D(const MeshType& mesh, const CellValue<Rdxd>& tensor) const
+  {
+    CellValue<R3x3> tensor3d{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { tensor3d[j] = to3D(tensor[j]); });
+    return tensor3d;
+  }
+
+  CellValue<const R3>
+  to3D(const MeshType& mesh, const CellValue<Rd> vector) const
+  {
+    CellValue<R3> vector3d{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { vector3d[j] = to3D(vector[j]); });
+    return vector3d;
+  }
+
+  Rdxd
+  toDimension(const R3x3 tensor3d) const
+  {
+    Rdxd tensor = zero;
+    for (size_t i = 0; i < Dimension; i++) {
+      for (size_t j = 0; j < Dimension; j++) {
+        tensor(i, j) = tensor3d(i, j);
+      }
+    }
+    return tensor;
+  }
+
+  Rd
+  toDimension(const R3 vector3d) const
+  {
+    Rd vector = zero;
+    for (size_t i = 0; i < Dimension; i++) {
+      vector[i] = vector3d[i];
+    }
+    return vector;
+  }
+
+  CellValue<const Rdxd>
+  toDimension(const MeshType& mesh, const CellValue<R3x3>& tensor3d) const
+  {
+    CellValue<Rdxd> tensor{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { tensor[j] = toDimension(tensor3d[j]); });
+    return tensor;
+  }
+
+  CellValue<const Rd>
+  toDimension(const MeshType& mesh, const CellValue<R3>& vector3d) const
+  {
+    CellValue<Rd> vector{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { vector[j] = to3D(vector3d[j]); });
+    return vector;
+  }
+
+  double
+  _compute_chi(const R3x3 S, const double yield) const
+  {
+    const R3x3 Id       = identity;
+    const R3x3 Sd       = S - 1. / 3 * trace(S) * Id;
+    double chi          = 0;
+    const double normS2 = frobeniusNorm(Sd) * frobeniusNorm(Sd);
+    double limit2       = 2. / 3 * yield * yield;
+    chi                 = std::max(0., std::sqrt(normS2) - std::sqrt(limit2));
+    return chi;
+  }
+
+  double
+  _compute_chi2(const R3x3 S, const double yield) const
+  {
+    const R3x3 Id       = identity;
+    const R3x3 Sd       = S - 1. / 3 * trace(S) * Id;
+    double chi          = 0;
+    const double normS2 = frobeniusNorm(Sd) * frobeniusNorm(Sd);
+    double limit2       = 2. / 3 * yield * yield;
+    double alpha        = (normS2 - limit2);
+    double sign         = alpha / std::abs(alpha);
+    if (sign < 0.) {
+      return chi;
+    }
+    int power    = 10;
+    double ratio = normS2 / limit2;
+    chi          = ratio;
+    for (int i = 0; i < power; i++) {
+      chi *= ratio;
+    }
+    return chi;
+  }
+
+  DiscreteFunctionDPk<Dimension, Rd>
+  _limit(const MeshType& mesh,
+         const DiscreteFunctionP0<const Rd>& f,
+         const DiscreteFunctionDPk<Dimension, const Rd>& DPk_fh,
+         const double eps) const
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+    DiscreteFunctionDPk<Dimension, Rd> DPk_fh_lim = copy(DPk_fh);
+    const auto xj                                 = mesh_data.xj();
+
+    const double eps2 = 1E-14;
+    double cmin       = 1. - eps;
+    double cmax       = 1. + eps;
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const Rd fj             = f[cell_id];
+        const auto cell_stencil = stencil[cell_id];
+        Rd coef1;
+        Rd coef2;
+        for (size_t dim = 0; dim < Dimension; ++dim) {
+          coef1[dim] = 1.;
+          coef2[dim] = 1.;
+        }
+        for (size_t dim = 0; dim < Dimension; ++dim) {
+          double f_min = std::min(cmin * fj[dim], cmax * fj[dim]);
+          double f_max = std::max(cmin * fj[dim], cmax * fj[dim]);
+
+          for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+            f_min = std::min(f_min, std::min(cmin * f[cell_stencil[i_cell]][dim], cmax * f[cell_stencil[i_cell]][dim]));
+            f_max = std::max(f_max, std::max(cmin * f[cell_stencil[i_cell]][dim], cmax * f[cell_stencil[i_cell]][dim]));
+          }
+          double f_bar_min = fj[dim];
+          double f_bar_max = fj[dim];
+
+          for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+            const CellId cell_k_id = cell_stencil[i_cell];
+            const double f_xk      = DPk_fh[cell_id](xj[cell_k_id])[dim];
+
+            f_bar_min = std::min(f_bar_min, f_xk);
+            f_bar_max = std::max(f_bar_max, f_xk);
+          }
+          if (std::abs(f_bar_max - fj[dim]) > eps2) {
+            coef1[dim] = (f_max - fj[dim]) / ((f_bar_max - fj[dim]));
+          }
+
+          if (std::abs(f_bar_min - fj[dim]) > eps2) {
+            coef2[dim] = (f_min - fj[dim]) / ((f_bar_min - fj[dim]));
+          }
+        }
+        double min_coef1 = min(coef1);
+        double min_coef2 = min(coef2);
+
+        const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
+
+        auto coefficients = DPk_fh_lim.coefficients(cell_id);
+
+        coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+        for (size_t i = 1; i < coefficients.size(); ++i) {
+          coefficients[i] *= lambda;
+        }
+      });
+    synchronize(DPk_fh_lim.cellArrays());
+    return DPk_fh_lim;
+  }
+
+  double
+  _min(const Rdxd M) const
+  {
+    double min = M(0, 0);
+    for (size_t i = 0; i < Dimension; ++i) {
+      for (size_t j = 0; j < Dimension; ++j) {
+        min = std::min(min, M(i, j));
+      }
+    }
+    return min;
+  }
+
+  double
+  _max(const Rdxd M) const
+  {
+    double max = M(0, 0);
+    for (size_t i = 0; i < Dimension; ++i) {
+      for (size_t j = 0; j < Dimension; ++j) {
+        max = std::max(max, M(i, j));
+      }
+    }
+    return max;
+  }
+
+  DiscreteFunctionDPk<Dimension, Rdxd>
+  _limit(const MeshType& mesh,
+         const DiscreteFunctionP0<const Rdxd>& f,
+         const DiscreteFunctionDPk<Dimension, const Rdxd>& DPk_fh,
+         const double eps) const
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+    DiscreteFunctionDPk<Dimension, Rdxd> DPk_fh_lim = copy(DPk_fh);
+    const auto xj                                   = mesh_data.xj();
+
+    const double eps2 = 1E-14;
+    double cmin       = 1. - eps;
+    double cmax       = 1. + eps;
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const Rdxd fj           = f[cell_id];
+        const auto cell_stencil = stencil[cell_id];
+        Rdxd coef1;
+        Rdxd coef2;
+        for (size_t i = 0; i < Dimension; ++i) {
+          for (size_t j = 0; j < Dimension; ++j) {
+            coef1(i, j) = 1.;
+            coef2(i, j) = 1.;
+          }
+        }
+        for (size_t i = 0; i < Dimension; ++i) {
+          for (size_t j = 0; j < Dimension; ++j) {
+            double f_min = std::min(cmin * fj(i, j), cmax * fj(i, j));
+            double f_max = std::max(cmin * fj(i, j), cmax * fj(i, j));
+
+            for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+              f_min =
+                std::min(f_min, std::min(cmin * f[cell_stencil[i_cell]](i, j), cmax * f[cell_stencil[i_cell]](i, j)));
+              f_max =
+                std::max(f_max, std::max(cmin * f[cell_stencil[i_cell]](i, j), cmax * f[cell_stencil[i_cell]](i, j)));
+            }
+            double f_bar_min = fj(i, j);
+            double f_bar_max = fj(i, j);
+
+            for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+              const CellId cell_k_id = cell_stencil[i_cell];
+              const double f_xk      = DPk_fh[cell_id](xj[cell_k_id])(i, j);
+
+              f_bar_min = std::min(f_bar_min, f_xk);
+              f_bar_max = std::max(f_bar_max, f_xk);
+            }
+            if (std::abs(f_bar_max - fj(i, j)) > eps2) {
+              coef1(i, j) = (f_max - fj(i, j)) / ((f_bar_max - fj(i, j)));
+            }
+
+            if (std::abs(f_bar_min - fj(i, j)) > eps2) {
+              coef2(i, j) = (f_min - fj(i, j)) / ((f_bar_min - fj(i, j)));
+            }
+          }
+        }
+        double min_coef1 = _min(coef1);
+        double min_coef2 = _min(coef2);
+
+        const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
+
+        auto coefficients = DPk_fh_lim.coefficients(cell_id);
+
+        coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+        for (size_t i = 1; i < coefficients.size(); ++i) {
+          coefficients[i] *= lambda;
+        }
+      });
+
+    synchronize(DPk_fh_lim.cellArrays());
+    return DPk_fh_lim;
+  }
+
+  NodeValuePerCell<const Rdxd>
+  _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) const
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+
+    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
+    const Rdxd I = identity;
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const size_t& nb_nodes = Ajr.numberOfSubValues(j);
+        const double& rhoaL_j  = rhoaL[j];
+        const double& rhoaT_j  = rhoaT[j];
+        for (size_t r = 0; r < nb_nodes; ++r) {
+          const Rdxd& M = tensorProduct(Cjr(j, r), njr(j, r));
+          Ajr(j, r)     = rhoaL_j * M + rhoaT_j * (l2Norm(Cjr(j, r)) * I - M);
+        }
+      });
+
+    return Ajr;
+  }
+
+  NodeValuePerCell<const Rdxd>
+  _computeEucclhydAjr(const MeshType& mesh,
+                      const DiscreteScalarFunction& rhoaL,
+                      const DiscreteScalarFunction& rhoaT) const
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+    const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
+
+    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
+
+    parallel_for(
+      Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
+    const Rdxd I = identity;
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        const auto& cell_faces = cell_to_face_matrix[j];
+
+        const double& rho_aL = rhoaL[j];
+        const double& rho_aT = rhoaT[j];
+
+        for (size_t L = 0; L < cell_faces.size(); ++L) {
+          const FaceId& l        = cell_faces[L];
+          const auto& face_nodes = face_to_node_matrix[l];
+
+          auto local_node_number_in_cell = [&](NodeId node_number) {
+            for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+              if (node_number == cell_nodes[i_node]) {
+                return i_node;
+              }
+            }
+            return std::numeric_limits<size_t>::max();
+          };
+
+          for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+            const size_t R = local_node_number_in_cell(face_nodes[rl]);
+            const Rdxd& M  = tensorProduct(Nlr(l, rl), nlr(l, rl));
+            Ajr(j, R) += rho_aL * M + rho_aT * (l2Norm(Nlr(l, rl)) * I - M);
+          }
+        }
+      });
+
+    return Ajr;
+  }
+
+  NodeValuePerCell<const Rdxd>
+  _computeAjr(const SolverType& solver_type,
+              const MeshType& mesh,
+              const DiscreteScalarFunction& rhoaL,
+              const DiscreteScalarFunction& rhoaT) const
+  {
+    if constexpr (Dimension == 1) {
+      return _computeGlaceAjr(mesh, rhoaL, rhoaT);
+    } else {
+      switch (solver_type) {
+      case SolverType::Glace: {
+        return _computeGlaceAjr(mesh, rhoaL, rhoaT);
+      }
+      case SolverType::Eucclhyd: {
+        return _computeEucclhydAjr(mesh, rhoaL, rhoaT);
+      }
+      default: {
+        throw UnexpectedError("invalid solver type");
+      }
+      }
+    }
+  }
+
+  NodeValue<Rdxd>
+  _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) const
+  {
+    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+    NodeValue<Rdxd> Ar{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        Rdxd sum                                   = zero;
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+        for (size_t j = 0; j < node_to_cell.size(); ++j) {
+          const CellId J       = node_to_cell[j];
+          const unsigned int R = node_local_number_in_its_cells[j];
+          sum += Ajr(J, R);
+        }
+        Ar[r] = sum;
+      });
+
+    return Ar;
+  }
+
+  NodeValue<Rd>
+  _computeBr(const Mesh<Dimension>& mesh,
+             const NodeValuePerCell<const Rdxd>& Ajr,
+             const DiscreteFunctionDPk<Dimension, const Rd>& DPk_u,
+             const DiscreteFunctionDPk<Dimension, const Rdxd>& DPk_sigma) const
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+
+    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+    NodeValue<Rd> b{mesh.connectivity()};
+    auto xr = mesh.xr();
+
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+        Rd br = zero;
+        for (size_t j = 0; j < node_to_cell.size(); ++j) {
+          const CellId J       = node_to_cell[j];
+          const unsigned int R = node_local_number_in_its_cells[j];
+          br += Ajr(J, R) * DPk_u[J](xr[r]) - DPk_sigma[J](xr[r]) * Cjr(J, R);
+        }
+
+        b[r] = br;
+      });
+
+    return b;
+  }
+
+  BoundaryConditionList
+  _getBCList(const MeshType& mesh,
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+  {
+    BoundaryConditionList bc_list;
+
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      bool is_valid_boundary_condition = true;
+
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        bc_list.emplace_back(
+          SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::fixed: {
+        bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::dirichlet: {
+        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+        if (dirichlet_bc_descriptor.name() == "velocity") {
+          MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+
+          Array<const Rd> value_list =
+            InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
+                                                                               mesh.xr(),
+                                                                               mesh_node_boundary.nodeList());
+
+          bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
+        } else if (dirichlet_bc_descriptor.name() == "pressure") {
+          const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+          if constexpr (Dimension == 1) {
+            MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            Array<const double> node_values =
+              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
+                                                                                     mesh_node_boundary.nodeList());
+
+            bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
+          } else {
+            MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+            Array<const double> face_values =
+              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
+                                                                                     mesh_face_boundary.faceList());
+            bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
+          }
+
+        } else if (dirichlet_bc_descriptor.name() == "normal-stress") {
+          const FunctionSymbolId normal_stress_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+          if constexpr (Dimension == 1) {
+            MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            Array<const Rdxd> node_values =
+              InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::node>(normal_stress_id, mesh.xr(),
+                                                                                   mesh_node_boundary.nodeList());
+
+            bc_list.emplace_back(NormalStressBoundaryCondition{mesh_node_boundary, node_values});
+          } else {
+            MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+            Array<const Rdxd> face_values =
+              InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::face>(normal_stress_id, mesh_data.xl(),
+                                                                                   mesh_face_boundary.faceList());
+            bc_list.emplace_back(NormalStressBoundaryCondition{mesh_face_boundary, face_values});
+          }
+
+        } else {
+          is_valid_boundary_condition = false;
+        }
+        break;
+      }
+      default: {
+        is_valid_boundary_condition = false;
+      }
+      }
+      if (not is_valid_boundary_condition) {
+        std::ostringstream error_msg;
+        error_msg << *bc_descriptor << " is an invalid boundary condition for hyperplastic solver";
+        throw NormalError(error_msg.str());
+      }
+    }
+
+    return bc_list;
+  }
+
+  void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
+  void _applyNormalStressBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
+  void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
+  void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
+  void
+  _applyBoundaryConditions(const BoundaryConditionList& bc_list,
+                           const MeshType& mesh,
+                           NodeValue<Rdxd>& Ar,
+                           NodeValue<Rd>& br) const
+  {
+    this->_applyPressureBC(bc_list, mesh, br);
+    this->_applyNormalStressBC(bc_list, mesh, br);
+    this->_applySymmetryBC(bc_list, Ar, br);
+    this->_applyVelocityBC(bc_list, Ar, br);
+  }
+
+  NodeValue<const Rd>
+  _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
+  {
+    NodeValue<Rd> u{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });
+
+    return u;
+  }
+
+  NodeValuePerCell<Rd>
+  _computeFjr(const MeshType& mesh,
+              const NodeValuePerCell<const Rdxd>& Ajr,
+              const NodeValue<const Rd>& ur,
+              const DiscreteFunctionDPk<Dimension, const Rd> DPk_uh,
+              const DiscreteFunctionDPk<Dimension, const Rdxd> DPk_sigmah) const
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const NodeValue<const Rd>& xr   = mesh.xr();
+
+    NodeValuePerCell<Rd> F{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        for (size_t r = 0; r < cell_nodes.size(); ++r) {
+          const NodeId R = cell_nodes[r];
+          F(j, r)        = -Ajr(j, r) * (DPk_uh[j](xr[R]) - ur[R]) + DPk_sigmah[j](xr[R]) * Cjr(j, r);
+        }
+      });
+
+    return F;
+  }
+
+ public:
+  std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+  compute_fluxes(const SolverType& solver_type,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+  {
+    std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+    if (not mesh_v) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperplastic solver expects P0 functions");
+    }
+
+    const MeshType& mesh                    = *mesh_v->get<MeshType>();
+    const DiscreteScalarFunction& rho       = rho_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u         = u_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& aL        = aL_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& aT        = aT_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction3d& sigma3d = sigma_v->get<DiscreteTensorFunction3d>();
+    const R3x3 I                            = identity;
+
+    //    std::shared_ptr<const MeshType> p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+    CellValue<R3x3> tensor_values3d{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { tensor_values3d[j] = sigma3d[j]; });
+    CellValue<const Rdxd> tensor_values = toDimension(mesh, tensor_values3d);
+    DiscreteTensorFunction sigma(mesh_v, tensor_values);
+
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+
+    for (auto&& bc_descriptor : bc_descriptor_list) {
+      if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+        symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+
+    PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list);
+
+    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u, sigma);
+    auto DPk_uh         = reconstruction[0]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto DPk_sigmah     = reconstruction[1]->template get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+
+    // DiscreteFunctionDPk<Dimension, Rd> u_lim       = copy(DPk_uh);
+    // DiscreteFunctionDPk<Dimension, Rdxd> sigma_lim = copy(DPk_sigmah);
+    double epsilon = 1.e-6;
+    auto u_lim     = _limit(mesh, u, DPk_uh, epsilon);
+    auto sigma_lim = _limit(mesh, sigma, DPk_sigmah, epsilon);
+
+    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
+
+    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
+
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+
+    synchronize(Ar);
+    synchronize(br);
+
+    NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
+    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
+
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes(const double& dt,
+               const MeshType& mesh,
+               const DiscreteScalarFunction& rho,
+               const DiscreteVectorFunction& u,
+               const DiscreteScalarFunction& E,
+               const DiscreteTensorFunction3d& Fe,
+               const DiscreteScalarFunction& zeta,
+               const DiscreteScalarFunction& yield,
+               const DiscreteTensorFunction3d& sigma,
+               const NodeValue<const Rd>& ur,
+               const NodeValuePerCell<const Rd>& Fjr) const
+  {
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
+      throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+    }
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
+    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+
+    CellValue<double> new_rho   = copy(rho.cellValues());
+    CellValue<Rd> new_u         = copy(u.cellValues());
+    CellValue<double> new_E     = copy(E.cellValues());
+    CellValue<R3x3> new_Fe      = copy(Fe.cellValues());
+    CellValue<double> new_zeta  = copy(zeta.cellValues());
+    CellValue<double> new_yield = copy(yield.cellValues());
+    CellValue<R3x3> sigma_s     = copy(sigma.cellValues());
+    CellValue<double> chi{mesh.connectivity()};
+    CellValue<double> rationorm{mesh.connectivity()};
+    chi.fill(0.);
+    rationorm.fill(0.);
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        Rdxd gradv           = zero;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += dot(Fjr(j, R), ur[r]);
+        }
+        R3x3 I                    = identity;
+        R3x3 gradv3d              = to3D(gradv);
+        R3x3 sigmas               = sigma[j] - 1. / 3 * trace(sigma[j]) * I;
+        const R3x3 elastic_fluxes = gradv3d * Fe[j];
+        const double dt_over_Mj   = dt / (rho[j] * Vj[j]);
+        const double dt_over_Vj   = dt / Vj[j];
+        new_u[j] += dt_over_Mj * momentum_fluxes;
+        new_E[j] += dt_over_Mj * energy_fluxes;
+        new_Fe[j] += dt_over_Vj * elastic_fluxes;
+        const double fenorm0 = frobeniusNorm(new_Fe[j]);
+        chi[j]               = _compute_chi(sigma[j], yield[j]);
+        const R3x3 M         = -dt_over_Vj * chi[j] * zeta[j] * sigmas;
+
+        new_Fe[j]            = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
+        const double fenorm1 = frobeniusNorm(new_Fe[j]);
+        rationorm[j]         = fenorm1 / fenorm0;
+      });
+    std::cout << "sum_chi " << sum(chi) << " min norm ratio " << min(rationorm) << " max norm ratio " << max(rationorm)
+              << "\n";
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {std::make_shared<MeshVariant>(new_mesh),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction3d(new_mesh, new_Fe))};
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes2(const double& dt,
+                const MeshType& mesh,
+                const DiscreteScalarFunction& rho,
+                const DiscreteVectorFunction& u,
+                const DiscreteScalarFunction& E,
+                const DiscreteTensorFunction3d& Fe,
+                const DiscreteScalarFunction& zeta,
+                const DiscreteScalarFunction& yield,
+                const DiscreteTensorFunction3d& sigma,
+                const NodeValue<const Rd>& ur,
+                const NodeValuePerCell<const Rd>& Fjr) const
+  {
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
+      throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+    }
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
+    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+
+    CellValue<double> new_rho   = copy(rho.cellValues());
+    CellValue<Rd> new_u         = copy(u.cellValues());
+    CellValue<double> new_E     = copy(E.cellValues());
+    CellValue<R3x3> new_Fe      = copy(Fe.cellValues());
+    CellValue<double> new_zeta  = copy(zeta.cellValues());
+    CellValue<double> new_yield = copy(yield.cellValues());
+    CellValue<R3x3> sigma_s     = copy(sigma.cellValues());
+    CellValue<double> chi{mesh.connectivity()};
+    CellValue<double> rationorm{mesh.connectivity()};
+    chi.fill(0.);
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        Rdxd gradv           = zero;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += dot(Fjr(j, R), ur[r]);
+        }
+        R3x3 I                    = identity;
+        R3x3 gradv3d              = to3D(gradv);
+        R3x3 sigmas               = sigma[j] - 1. / 3 * trace(sigma[j]) * I;
+        const R3x3 elastic_fluxes = gradv3d * Fe[j];
+        const double dt_over_Mj   = dt / (rho[j] * Vj[j]);
+        const double dt_over_Vj   = dt / Vj[j];
+        new_u[j] += dt_over_Mj * momentum_fluxes;
+        new_E[j] += dt_over_Mj * energy_fluxes;
+        new_Fe[j] += dt_over_Vj * elastic_fluxes;
+        chi[j]     = _compute_chi(sigma[j], yield[j]);
+        int sub_it = static_cast<int>(std::ceil(dt_over_Vj * chi[j] * zeta[j] * frobeniusNorm(sigmas)));
+        if (sub_it > 1) {
+          std::cout << sub_it << " dt_over_Vj " << dt_over_Vj << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
+                    << " frobeniusNorm(sigmas) " << frobeniusNorm(sigmas) << "\n";
+        }
+        for (int i = 0; i < sub_it; i++) {
+          const R3x3 M = Vj[j] * I + dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * sigmas;
+          new_Fe[j]    = Vj[j] * inverse(M) * new_Fe[j];
+        }
+      });
+    std::cout << "sum_chi " << sum(chi) << "\n";
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {std::make_shared<MeshVariant>(new_mesh),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction3d(new_mesh, new_Fe))};
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes_B(const double& dt,
+                 const MeshType& mesh,
+                 const DiscreteScalarFunction& rho,
+                 const DiscreteVectorFunction& u,
+                 const DiscreteScalarFunction& E,
+                 const DiscreteTensorFunction3d& Be,
+                 const DiscreteScalarFunction& zeta,
+                 const DiscreteScalarFunction& yield,
+                 const DiscreteTensorFunction3d& sigma,
+                 const NodeValue<const Rd>& ur,
+                 const NodeValuePerCell<const Rd>& Fjr) const
+  {
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
+      throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+    }
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
+    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+
+    CellValue<double> new_rho   = copy(rho.cellValues());
+    CellValue<Rd> new_u         = copy(u.cellValues());
+    CellValue<double> new_E     = copy(E.cellValues());
+    CellValue<R3x3> new_Be      = copy(Be.cellValues());
+    CellValue<double> new_zeta  = copy(zeta.cellValues());
+    CellValue<double> new_yield = copy(yield.cellValues());
+    CellValue<R3x3> sigma_s     = copy(sigma.cellValues());
+    CellValue<double> chi{mesh.connectivity()};
+    CellValue<double> rationorm{mesh.connectivity()};
+    chi.fill(0.);
+    // rationorm.fill(0.);
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+        Rd momentum_fluxes     = zero;
+        double energy_fluxes   = 0;
+        Rdxd gradv             = zero;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += dot(Fjr(j, R), ur[r]);
+        }
+        R3x3 I       = identity;
+        R3x3 gradv3d = to3D(gradv);
+        R3x3 sigmas  = sigma[j] - 1. / 3 * trace(sigma[j]) * I;
+        //        const R3x3 elastic_fluxes = gradv3d + transpose(gradv3d);
+        const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+        const double dt_over_Vj = dt / Vj[j];
+        new_u[j] += dt_over_Mj * momentum_fluxes;
+        new_E[j] += dt_over_Mj * energy_fluxes;
+        // new_Be[j] += dt_over_Vj * elastic_fluxes;
+        // const double fenorm0 = frobeniusNorm(new_Be[j]);
+        chi[j]           = _compute_chi(sigma[j], yield[j]);
+        const R3x3 M     = dt_over_Vj * (gradv3d - chi[j] * zeta[j] * sigmas);
+        const R3x3 expM  = EigenvalueSolver{}.computeExpForTinyMatrix(M);
+        const R3x3 expMT = EigenvalueSolver{}.computeExpForTinyMatrix(transpose(M));
+        new_Be[j]        = expM * Be[j] * expMT;
+        // const double fenorm1 = frobeniusNorm(new_Be[j]);
+        // rationorm[j]         = fenorm1 / fenorm0;
+      });
+    std::cout << "sum_chi " << sum(chi) << "\n";
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {std::make_shared<MeshVariant>(new_mesh),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction3d(new_mesh, new_Be))};
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes(const double& dt,
+               const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+               const std::shared_ptr<const DiscreteFunctionVariant>& u,
+               const std::shared_ptr<const DiscreteFunctionVariant>& E,
+               const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+               const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+               const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+               const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+               const std::shared_ptr<const ItemValueVariant>& ur,
+               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+               const RelaxationType& relaxation_type) const
+  {
+    std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
+    if (not mesh_v) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperplastic solver expects P0 functions");
+    }
+    switch (relaxation_type) {
+    case RelaxationType::Exponential: {
+      return this->apply_fluxes(dt,                                       //
+                                *mesh_v->get<MeshType>(),                 //
+                                rho->get<DiscreteScalarFunction>(),       //
+                                u->get<DiscreteVectorFunction>(),         //
+                                E->get<DiscreteScalarFunction>(),         //
+                                Fe->get<DiscreteTensorFunction3d>(),      //
+                                zeta->get<DiscreteScalarFunction>(),      //
+                                yield->get<DiscreteScalarFunction>(),     //
+                                sigma->get<DiscreteTensorFunction3d>(),   //
+                                ur->get<NodeValue<const Rd>>(),           //
+                                Fjr->get<NodeValuePerCell<const Rd>>());
+    }
+    case RelaxationType::Implicit: {
+      return this->apply_fluxes2(dt,                                       //
+                                 *mesh_v->get<MeshType>(),                 //
+                                 rho->get<DiscreteScalarFunction>(),       //
+                                 u->get<DiscreteVectorFunction>(),         //
+                                 E->get<DiscreteScalarFunction>(),         //
+                                 Fe->get<DiscreteTensorFunction3d>(),      //
+                                 zeta->get<DiscreteScalarFunction>(),      //
+                                 yield->get<DiscreteScalarFunction>(),     //
+                                 sigma->get<DiscreteTensorFunction3d>(),   //
+                                 ur->get<NodeValue<const Rd>>(),           //
+                                 Fjr->get<NodeValuePerCell<const Rd>>());
+    }
+    case RelaxationType::CauchyGreen: {
+      return this->apply_fluxes_B(dt,                                       //
+                                  *mesh_v->get<MeshType>(),                 //
+                                  rho->get<DiscreteScalarFunction>(),       //
+                                  u->get<DiscreteVectorFunction>(),         //
+                                  E->get<DiscreteScalarFunction>(),         //
+                                  Fe->get<DiscreteTensorFunction3d>(),      //
+                                  zeta->get<DiscreteScalarFunction>(),      //
+                                  yield->get<DiscreteScalarFunction>(),     //
+                                  sigma->get<DiscreteTensorFunction3d>(),   //
+                                  ur->get<NodeValue<const Rd>>(),           //
+                                  Fjr->get<NodeValuePerCell<const Rd>>());
+    }
+    default: {
+      throw UnexpectedError("invalid relaxation type");
+    }
+    }
+  }
+
+  std::shared_ptr<const DiscreteFunctionVariant>
+  apply_plastic_relaxation(const RelaxationType& relaxation_type,
+                           const double& dt,
+                           const std::shared_ptr<const MeshType>& mesh,
+                           const DiscreteTensorFunction3d& Fe,
+                           const DiscreteScalarFunction& zeta,
+                           const DiscreteScalarFunction& yield,
+                           const DiscreteTensorFunction3d& sigman,
+                           const DiscreteTensorFunction3d& sigmanp1) const
+  {
+    CellValue<R3x3> new_Fe = copy(Fe.cellValues());
+    CellValue<double> chi{mesh->connectivity()};
+    chi.fill(0.);
+    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+    R3x3 I                     = identity;
+    switch (relaxation_type) {
+    case RelaxationType::Exponential: {
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+          const double dt_over_Vj = dt / Vj[j];
+          R3x3 sigmasn            = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
+          R3x3 sigmasnp1          = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
+          chi[j]                  = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
+          const R3x3 M            = -dt_over_Vj * chi[j] * zeta[j] * 0.5 * (sigmasn + sigmasnp1);
+          new_Fe[j]               = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
+        });
+      break;
+    }
+    case RelaxationType::Implicit: {
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+          const double dt_over_Vj = dt / Vj[j];
+          R3x3 sigmasn            = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
+          R3x3 sigmasnp1          = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
+          chi[j]                  = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
+          int sub_it =
+            static_cast<int>(std::ceil(dt_over_Vj * chi[j] * zeta[j] * frobeniusNorm(0.5 * (sigmasn + sigmasnp1))));
+          if (sub_it > 1) {
+            std::cout << sub_it << " dt_over_Vj " << dt_over_Vj << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
+                      << " frobeniusNorm(sigmas) " << frobeniusNorm(0.5 * (sigmasn + sigmasnp1)) << "\n";
+          }
+          for (int i = 0; i < sub_it; i++) {
+            const R3x3 M =
+              Vj[j] * I + 0.5 * dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * (sigmasn + sigmasnp1);
+            new_Fe[j] = Vj[j] * inverse(M) * new_Fe[j];
+          }
+        });
+      break;
+    }
+    default: {
+      throw UnexpectedError("invalid relaxation type");
+    }
+    }
+
+    // for (CellId j = 0; j < mesh->numberOfCells(); j++) {
+    //   if ((std::abs(frobeniusNorm(new_Fe[j]) - std::sqrt(3.)) > 1.e-1) or (j == 6399)) {
+    //     std::cout << j << " new_Fe " << new_Fe[j] << " chi " << chi[j] << "\n";
+    //   }
+    // }
+    std::cout << "sum_chi " << sum(chi) << "\n";
+
+    return {std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction3d(mesh, new_Fe))};
+  }
+
+  std::shared_ptr<const DiscreteFunctionVariant>
+  apply_plastic_relaxation(const RelaxationType& relaxation_type,
+                           const double& dt,
+                           const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                           const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                           const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                           const std::shared_ptr<const DiscreteFunctionVariant>& sigman,
+                           const std::shared_ptr<const DiscreteFunctionVariant>& sigmanp1) const
+  {
+    std::shared_ptr mesh_v = getCommonMesh({sigman, sigmanp1});
+    if (not mesh_v) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({sigman, sigmanp1}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperplastic solver expects P0 functions");
+    }
+
+    return this->apply_plastic_relaxation(relaxation_type, dt,                    //
+                                          mesh_v->get<MeshType>(),                //
+                                          Fe->get<DiscreteTensorFunction3d>(),    //
+                                          zeta->get<DiscreteScalarFunction>(),    //
+                                          yield->get<DiscreteScalarFunction>(),   //
+                                          sigman->get<DiscreteTensorFunction3d>(),
+                                          sigmanp1->get<DiscreteTensorFunction3d>());
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_elastic_fluxes(const double& dt,
+                       const MeshType& mesh,
+                       const DiscreteScalarFunction& rho,
+                       const DiscreteVectorFunction& u,
+                       const DiscreteScalarFunction& E,
+                       const DiscreteTensorFunction3d& Fe,
+                       const NodeValue<const Rd>& ur,
+                       const NodeValuePerCell<const Rd>& Fjr) const
+  {
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
+      throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+    }
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
+    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+
+    CellValue<double> new_rho = copy(rho.cellValues());
+    CellValue<Rd> new_u       = copy(u.cellValues());
+    CellValue<double> new_E   = copy(E.cellValues());
+    CellValue<R3x3> new_Fe    = copy(Fe.cellValues());
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        Rdxd gradv           = zero;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += dot(Fjr(j, R), ur[r]);
+        }
+        R3x3 gradv3d              = to3D(gradv);
+        const R3x3 elastic_fluxes = gradv3d * Fe[j];
+        const double dt_over_Mj   = dt / (rho[j] * Vj[j]);
+        const double dt_over_Vj   = dt / Vj[j];
+        new_u[j] += dt_over_Mj * momentum_fluxes;
+        new_E[j] += dt_over_Mj * energy_fluxes;
+        new_Fe[j] += dt_over_Vj * elastic_fluxes;
+      });
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {std::make_shared<MeshVariant>(new_mesh),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction3d(new_mesh, new_Fe))};
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_elastic_fluxes(const double& dt,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                       const std::shared_ptr<const ItemValueVariant>& ur,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
+  {
+    std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
+    if (not mesh_v) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperplastic solver expects P0 functions");
+    }
+
+    return this->apply_elastic_fluxes(dt,                                    //
+                                      *mesh_v->get<MeshType>(),              //
+                                      rho->get<DiscreteScalarFunction>(),    //
+                                      u->get<DiscreteVectorFunction>(),      //
+                                      E->get<DiscreteScalarFunction>(),      //
+                                      Fe->get<DiscreteTensorFunction3d>(),   //
+                                      ur->get<NodeValue<const Rd>>(),        //
+                                      Fjr->get<NodeValuePerCell<const Rd>>());
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply(const SolverType& solver_type,
+        const RelaxationType& relaxation_type,
+        const double& dt,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E,
+        const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+        const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+        const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+  {
+    auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
+    return apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, relaxation_type);
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_elastic(const SolverType& solver_type,
+                const double& dt,
+                const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+  {
+    auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
+    return apply_elastic_fluxes(dt, rho, u, E, Fe, ur, Fjr);
+  }
+
+  HyperplasticSolverO2()                       = default;
+  HyperplasticSolverO2(HyperplasticSolverO2&&) = default;
+  ~HyperplasticSolverO2()                      = default;
+};
+
+template <MeshConcept MeshType>
+void
+HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
+                                                                              const MeshType& mesh,
+                                                                              NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+          if constexpr (Dimension == 1) {
+            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+            const auto& node_list  = bc.nodeList();
+            const auto& value_list = bc.valueList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                const NodeId node_id       = node_list[i_node];
+                const auto& node_cell_list = node_to_cell_matrix[node_id];
+                Assert(node_cell_list.size() == 1);
+
+                CellId node_cell_id              = node_cell_list[0];
+                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+                br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+              });
+          } else {
+            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
+            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
+            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();
+
+            const auto& face_list  = bc.faceList();
+            const auto& value_list = bc.valueList();
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id       = face_list[i_face];
+              const auto& face_cell_list = face_to_cell_matrix[face_id];
+              Assert(face_cell_list.size() == 1);
+
+              CellId face_cell_id              = face_cell_list[0];
+              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+              const auto& face_nodes = face_to_node_matrix[face_id];
+
+              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+                NodeId node_id = face_nodes[i_node];
+                br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+              }
+            }
+          }
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <MeshConcept MeshType>
+void
+HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::_applyNormalStressBC(const BoundaryConditionList& bc_list,
+                                                                                  const MeshType& mesh,
+                                                                                  NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<NormalStressBoundaryCondition, T>) {
+          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+          if constexpr (Dimension == 1) {
+            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+            const auto& node_list  = bc.nodeList();
+            const auto& value_list = bc.valueList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                const NodeId node_id       = node_list[i_node];
+                const auto& node_cell_list = node_to_cell_matrix[node_id];
+                Assert(node_cell_list.size() == 1);
+
+                CellId node_cell_id              = node_cell_list[0];
+                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+                br[node_id] += value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+              });
+          } else {
+            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
+            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
+            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();
+
+            const auto& face_list  = bc.faceList();
+            const auto& value_list = bc.valueList();
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id       = face_list[i_face];
+              const auto& face_cell_list = face_to_cell_matrix[face_id];
+              Assert(face_cell_list.size() == 1);
+
+              CellId face_cell_id              = face_cell_list[0];
+              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+              const auto& face_nodes = face_to_node_matrix[face_id];
+
+              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+                NodeId node_id = face_nodes[i_node];
+                br[node_id] += sign * value_list[i_face] * Nlr(face_id, i_node);
+              }
+            }
+          }
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <MeshConcept MeshType>
+void
+HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
+                                                                              NodeValue<Rdxd>& Ar,
+                                                                              NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+          const Rd& n = bc.outgoingNormal();
+
+          const Rdxd I   = identity;
+          const Rdxd nxn = tensorProduct(n, n);
+          const Rdxd P   = I - nxn;
+
+          const Array<const NodeId>& node_list = bc.nodeList();
+          parallel_for(
+            bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+              const NodeId r = node_list[r_number];
+
+              Ar[r] = P * Ar[r] * P + nxn;
+              br[r] = P * br[r];
+            });
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <MeshConcept MeshType>
+void
+HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
+                                                                              NodeValue<Rdxd>& Ar,
+                                                                              NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+          const auto& node_list  = bc.nodeList();
+          const auto& value_list = bc.valueList();
+
+          parallel_for(
+            node_list.size(), PUGS_LAMBDA(size_t i_node) {
+              NodeId node_id    = node_list[i_node];
+              const auto& value = value_list[i_node];
+
+              Ar[node_id] = identity;
+              br[node_id] = value;
+            });
+        } else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
+          const auto& node_list = bc.nodeList();
+          parallel_for(
+            node_list.size(), PUGS_LAMBDA(size_t i_node) {
+              NodeId node_id = node_list[i_node];
+
+              Ar[node_id] = identity;
+              br[node_id] = zero;
+            });
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <MeshConcept MeshType>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::FixedBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary m_mesh_node_boundary;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+
+  ~FixedBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::PressureBoundaryCondition
+{
+ private:
+  const MeshFaceBoundary m_mesh_face_boundary;
+  const Array<const double> m_value_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_mesh_face_boundary.faceList();
+  }
+
+  const Array<const double>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  PressureBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const double>& value_list)
+    : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
+  {}
+
+  ~PressureBoundaryCondition() = default;
+};
+
+template <>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2<Mesh<1>>::PressureBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary m_mesh_node_boundary;
+  const Array<const double> m_value_list;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  const Array<const double>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  PressureBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const double>& value_list)
+    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+  {}
+
+  ~PressureBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::NormalStressBoundaryCondition
+{
+ private:
+  const MeshFaceBoundary m_mesh_face_boundary;
+  const Array<const Rdxd> m_value_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_mesh_face_boundary.faceList();
+  }
+
+  const Array<const Rdxd>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  NormalStressBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const Rdxd>& value_list)
+    : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
+  {}
+
+  ~NormalStressBoundaryCondition() = default;
+};
+
+template <>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2<Mesh<1>>::NormalStressBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary m_mesh_node_boundary;
+  const Array<const Rdxd> m_value_list;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  const Array<const Rdxd>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  NormalStressBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const Rdxd>& value_list)
+    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+  {}
+
+  ~NormalStressBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::VelocityBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary m_mesh_node_boundary;
+
+  const Array<const TinyVector<Dimension>> m_value_list;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  const Array<const TinyVector<Dimension>>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  VelocityBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+                            const Array<const TinyVector<Dimension>>& value_list)
+    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+  {}
+
+  ~VelocityBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+
+ public:
+  const Rd&
+  outgoingNormal() const
+  {
+    return m_mesh_flat_node_boundary.outgoingNormal();
+  }
+
+  size_t
+  numberOfNodes() const
+  {
+    return m_mesh_flat_node_boundary.nodeList().size();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_flat_node_boundary.nodeList();
+  }
+
+  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
+    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  {
+    ;
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
+HyperplasticSolverO2Handler::HyperplasticSolverO2Handler(const std::shared_ptr<const MeshVariant>& mesh_v)
+{
+  if (not mesh_v) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::visit(
+    [&](auto&& mesh) {
+      using MeshType = mesh_type_t<decltype(mesh)>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        m_hyperplastic_solver = std::make_unique<HyperplasticSolverO2<MeshType>>();
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    mesh_v->variant());
+}
diff --git a/src/scheme/HyperplasticSolverO2.hpp b/src/scheme/HyperplasticSolverO2.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff89780d141f9aaea04221ca6be1ad3f46497085
--- /dev/null
+++ b/src/scheme/HyperplasticSolverO2.hpp
@@ -0,0 +1,131 @@
+#ifndef HYPERPLASTIC_SOLVER_O2_HPP
+#define HYPERPLASTIC_SOLVER_O2_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+class IBoundaryConditionDescriptor;
+class MeshVariant;
+class ItemValueVariant;
+class SubItemValuePerItemVariant;
+class DiscreteFunctionVariant;
+
+double hyperplastic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c);
+
+class HyperplasticSolverO2Handler
+{
+ public:
+  enum class SolverType
+  {
+    Glace,
+    Eucclhyd
+  };
+
+  enum class RelaxationType
+  {
+    Implicit,
+    Exponential,
+    CauchyGreen
+  };
+
+ private:
+  struct IHyperplasticSolverO2
+  {
+    virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>>
+    compute_fluxes(
+      const SolverType& solver_type,
+      const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+      const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+      const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+      const std::shared_ptr<const DiscreteFunctionVariant>& u,
+      const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+
+    virtual std::tuple<std::shared_ptr<const MeshVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>>
+    apply_fluxes(const double& dt,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                 const std::shared_ptr<const ItemValueVariant>& ur,
+                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+                 const RelaxationType& relaxation_type) const = 0;
+
+    virtual std::tuple<std::shared_ptr<const MeshVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>>
+    apply(const SolverType& solver_type,
+          const RelaxationType& relaxation_type,
+          const double& dt,
+          const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+          const std::shared_ptr<const DiscreteFunctionVariant>& u,
+          const std::shared_ptr<const DiscreteFunctionVariant>& E,
+          const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+          const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+          const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+          const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+          const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+          const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+
+    virtual std::tuple<std::shared_ptr<const MeshVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>>
+    apply_elastic(const SolverType& solver_type,
+                  const double& dt,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+
+    virtual std::shared_ptr<const DiscreteFunctionVariant> apply_plastic_relaxation(
+      const RelaxationType& relaxation_type,
+      const double& dt,
+      const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+      const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+      const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+      const std::shared_ptr<const DiscreteFunctionVariant>& sigman,
+      const std::shared_ptr<const DiscreteFunctionVariant>& sigmanp1) const = 0;
+
+    IHyperplasticSolverO2()                                   = default;
+    IHyperplasticSolverO2(IHyperplasticSolverO2&&)            = default;
+    IHyperplasticSolverO2& operator=(IHyperplasticSolverO2&&) = default;
+
+    virtual ~IHyperplasticSolverO2() = default;
+  };
+
+  template <MeshConcept MeshType>
+  class HyperplasticSolverO2;
+
+  std::unique_ptr<IHyperplasticSolverO2> m_hyperplastic_solver;
+
+ public:
+  const IHyperplasticSolverO2&
+  solver() const
+  {
+    return *m_hyperplastic_solver;
+  }
+
+  HyperplasticSolverO2Handler(const std::shared_ptr<const MeshVariant>& mesh);
+};
+
+#endif   // HYPERPLASTIC_SOLVER_O2_HPP
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index a5edec7477f9dadf71f2413da3e9c910d635d52c..2eaadd46a7cb0780d3d54fe8e7d06fc4e059a9b3 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -1,5 +1,6 @@
 #include <catch2/catch_approx.hpp>
 #include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
 
 #include <Kokkos_Core.hpp>
 
@@ -209,36 +210,26 @@ TEST_CASE("TinyMatrix", "[algebra]")
     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 dot product calculation")
+  SECTION("checking for doubleDot calculations")
   {
-    {
-      TinyMatrix<1, 1, int> M(6);
-      TinyMatrix<1, 1, int> N(7);
-      REQUIRE(dot(M, N) == trace(M * transpose(N)));
-    }
-
-    {
-      TinyMatrix<2, 3, int> M(6, 3, -2,   //
-                              5, -1, 4);
-      TinyMatrix<2, 3, int> N(7, 8, -6,   //
-                              -3, 4, 2);
-      REQUIRE(dot(M, N) == trace(M * transpose(N)));
-      REQUIRE(dot(M, N) == trace(N * transpose(M)));
-    }
+    REQUIRE(doubleDot(TinyMatrix<1, 1, int>(6), TinyMatrix<1, 1, int>(4)) == 24);
+    REQUIRE(doubleDot(TinyMatrix<2, 2>(5, 1, -3, 6), TinyMatrix<2, 2>(-2, 3, -5, 1)) ==
+            Catch::Approx(-10 + 3 + 15 + 6).epsilon(1E-14));
+    REQUIRE(doubleDot(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5),
+                      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 + 2.3 * 2.3 + 7 * 7 + 6.2 * 6.2 + 9 + 16 + 81 + 1 + 4.1 * 4.1 + 25 + 4 + 9 + 4 + 27 * 27 +
+                          9 + 17.5 * 17.5)
+              .epsilon(1E-14));
   }
 
-  SECTION("checking for Frobenius norm calculation")
+  SECTION("checking for norm calculations")
   {
-    {
-      TinyMatrix<1, 1, int> M(-6);
-      REQUIRE(frobeniusNorm(M) == 6);
-    }
-
-    {
-      TinyMatrix<2, 3, int> M(6, 3, -2,   //
-                              5, -1, 4);
-      REQUIRE(frobeniusNorm(M) == std::sqrt(dot(M, M)));
-    }
+    REQUIRE(frobeniusNorm(TinyMatrix<1, 1>(6)) == Catch::Approx(6));
+    REQUIRE(frobeniusNorm(TinyMatrix<2, 2>(5, 1, -3, 6)) == Catch::Approx(std::sqrt(25 + 1 + 9 + 36)).epsilon(1E-14));
+    REQUIRE(frobeniusNorm(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
+            Catch::Approx(std::sqrt(1 + 2.3 * 2.3 + 7 * 7 + 6.2 * 6.2 + 9 + 16 + 81 + 1 + 4.1 * 4.1 + 25 + 4 + 9 + 4 +
+                                    27 * 27 + 9 + 17.5 * 17.5))
+              .epsilon(1E-14));
   }
 
   SECTION("checking for inverse calculations")
@@ -272,27 +263,6 @@ TEST_CASE("TinyMatrix", "[algebra]")
       REQUIRE(I(2, 1) == Catch::Approx(0).margin(1E-14));
       REQUIRE(I(2, 2) == Catch::Approx(1).epsilon(1E-14));
     }
-    {
-      const TinyMatrix<4> A4(2, 3, 1, 4, -1, 5, -2, 3, 4, 1, 2, 3, 4, -1, -1, 0);
-      const TinyMatrix<4> I = inverse(A4) * A4;
-
-      REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
-      REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(0, 2) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(0, 3) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(1, 0) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(1, 1) == Catch::Approx(1).epsilon(1E-14));
-      REQUIRE(I(1, 2) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(1, 3) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(2, 0) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(2, 1) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(2, 2) == Catch::Approx(1).epsilon(1E-14));
-      REQUIRE(I(2, 3) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(3, 0) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(3, 1) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(3, 2) == Catch::Approx(0).margin(1E-14));
-      REQUIRE(I(3, 3) == Catch::Approx(1).epsilon(1E-14));
-    }
   }
 
   SECTION("checking for sizes")