diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 3e9d9ce9986e595bd69ec753d1d2b753631b824d..40808c6ab3165f00540025735a7bba6863886ea3 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -586,16 +586,16 @@ SchemeModule::SchemeModule()
                                  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::shared_ptr<const DiscreteFunctionVariant>& sigmad,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list)
                                 -> std::tuple<std::shared_ptr<const ItemValueVariant>,
                                               std::shared_ptr<const SubItemValuePerItemVariant>> {
-                                return HypoelasticSolverHandler{getCommonMesh({rho, aL, aT, u, sigma, p})}
+                                return HypoelasticSolverHandler{getCommonMesh({rho, aL, aT, u, sigmad, p})}
                                   .solver()
-                                  .compute_fluxes(HypoelasticSolverHandler::SolverType::Eucclhyd, rho, aL, aT, u, sigma,
-                                                  p, bc_descriptor_list);
+                                  .compute_fluxes(HypoelasticSolverHandler::SolverType::Eucclhyd, rho, aL, aT, u,
+                                                  sigmad, p, bc_descriptor_list);
                               }
 
                               ));
@@ -609,7 +609,7 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& CG,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& aL,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& aT,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
@@ -617,12 +617,11 @@ SchemeModule::SchemeModule()
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma, p})}
+                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigmad, p})}
                                   .solver()
-                                  .apply(HypoelasticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, CG, aL, aT,
-                                         sigma, p, bc_descriptor_list);
+                                  .apply(HypoelasticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, sigmad, aL, aT,
+                                         p, bc_descriptor_list);
                               }
 
                               ));
@@ -634,15 +633,16 @@ SchemeModule::SchemeModule()
                                  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::shared_ptr<const DiscreteFunctionVariant>& sigmad,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list)
                                 -> std::tuple<std::shared_ptr<const ItemValueVariant>,
                                               std::shared_ptr<const SubItemValuePerItemVariant>> {
-                                return HypoelasticSolverHandler{getCommonMesh({rho, aL, aT, u, sigma})}
+                                return HypoelasticSolverHandler{getCommonMesh({rho, aL, aT, u, sigmad, p})}
                                   .solver()
                                   .compute_fluxes(HypoelasticSolverHandler::SolverType::Glace,   //
-                                                  rho, aL, aT, u, sigma, bc_descriptor_list);
+                                                  rho, aL, aT, u, sigmad, p, bc_descriptor_list);
                               }
 
                               ));
@@ -656,7 +656,8 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& CG,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& aL,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& aT,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
                                  const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -664,10 +665,10 @@ SchemeModule::SchemeModule()
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
+                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigmad, p})}
                                   .solver()
-                                  .apply(HypoelasticSolverHandler::SolverType::Glace, dt, rho, u, E, CG, aL, aT, sigma,
-                                         p, bc_descriptor_list);
+                                  .apply(HypoelasticSolverHandler::SolverType::Glace, dt, rho, u, E, sigmad, aL, aT, p,
+                                         bc_descriptor_list);
                               }
 
                               ));
@@ -678,7 +679,8 @@ SchemeModule::SchemeModule()
                               [](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>& CG,       //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,   //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,        //
                                  const std::shared_ptr<const ItemValueVariant>& ur,              //
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
                                  const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -686,9 +688,9 @@ SchemeModule::SchemeModule()
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, CG})}   //
+                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, sigmad, p})}   //
                                   .solver()
-                                  .apply_fluxes(dt, rho, u, E, CG, ur, Fjr);
+                                  .apply_fluxes(dt, rho, u, E, sigmad, p, ur, Fjr);
                               }
 
                               ));
diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp
index 044c8978c5e43fd32ad606f6f70a7d5f740729fe..f39f214f0915b5cfc0634ee7310c5a61be01f1cb 100644
--- a/src/scheme/HypoelasticSolver.cpp
+++ b/src/scheme/HypoelasticSolver.cpp
@@ -399,25 +399,28 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
                  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::shared_ptr<const DiscreteFunctionVariant>& sigmad_v,
                  const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+    std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigmad_v});
     if (not i_mesh) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
-    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigmad_v}, DiscreteFunctionType::P0)) {
       throw NormalError("hypoelastic solver expects P0 functions");
     }
 
-    const MeshType& mesh                = dynamic_cast<const MeshType&>(*i_mesh);
-    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 DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+    const MeshType& mesh                 = dynamic_cast<const MeshType&>(*i_mesh);
+    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 DiscreteTensorFunction& sigmad = sigmad_v->get<DiscreteTensorFunction>();
+    const DiscreteScalarFunction& p      = p_v->get<DiscreteScalarFunction>();
+    const Rdxd I                         = identity;
+    const DiscreteTensorFunction& sigma  = -p * I + sigmad;
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
@@ -431,7 +434,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
     synchronize(br);
 
     NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
+    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigmad);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -447,7 +450,8 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
                const DiscreteScalarFunction& rho,
                const DiscreteVectorFunction& u,
                const DiscreteScalarFunction& E,
-               const DiscreteTensorFunction& sigma,
+               const DiscreteTensorFunction& sigmad,
+               // p probably not needed here
                const DiscreteScalarFunction& p,
                const NodeValue<const Rd>& ur,
                const NodeValuePerCell<const Rd>& Fjr) const
@@ -468,10 +472,10 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
     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<Rdxd> new_sigma = copy(sigma.cellValues());
+    CellValue<double> new_rho  = copy(rho.cellValues());
+    CellValue<Rd> new_u        = copy(u.cellValues());
+    CellValue<double> new_E    = copy(E.cellValues());
+    CellValue<Rdxd> new_sigmad = copy(sigmad.cellValues());
 
     parallel_for(
       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
@@ -505,7 +509,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
     return {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>(DiscreteTensorFunction(new_mesh, new_sigma))};
+            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_sigmad))};
   }
 
   std::tuple<std::shared_ptr<const IMesh>,
@@ -517,7 +521,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
                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>& sigma,
+               const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
                const std::shared_ptr<const DiscreteFunctionVariant>& p,
                const std::shared_ptr<const ItemValueVariant>& ur,
                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
@@ -536,7 +540,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
                               rho->get<DiscreteScalarFunction>(),       //
                               u->get<DiscreteVectorFunction>(),         //
                               E->get<DiscreteScalarFunction>(),         //
-                              sigma->get<DiscreteTensorFunction>(),     //
+                              sigmad->get<DiscreteTensorFunction>(),    //
                               p->get<DiscreteScalarFunction>(),         //
                               ur->get<NodeValue<const Rd>>(),           //
                               Fjr->get<NodeValuePerCell<const Rd>>());
@@ -552,14 +556,14 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
         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>& sigma,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
         const std::shared_ptr<const DiscreteFunctionVariant>& aL,
         const std::shared_ptr<const DiscreteFunctionVariant>& aT,
         const std::shared_ptr<const DiscreteFunctionVariant>& p,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
-    auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, p, bc_descriptor_list);
-    return apply_fluxes(dt, rho, u, E, sigma, p, ur, Fjr);
+    auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigmad, p, bc_descriptor_list);
+    return apply_fluxes(dt, rho, u, E, sigmad, p, ur, Fjr);
   }
 
   HypoelasticSolver()                    = default;
diff --git a/src/scheme/HypoelasticSolver.hpp b/src/scheme/HypoelasticSolver.hpp
index aa2501f3a2168bf27b28cda6a899a635837ff7b6..5c21a83ac60adf5bb98c22a27e14ee62242e3c6e 100644
--- a/src/scheme/HypoelasticSolver.hpp
+++ b/src/scheme/HypoelasticSolver.hpp
@@ -33,7 +33,7 @@ class HypoelasticSolverHandler
       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::shared_ptr<const DiscreteFunctionVariant>& sigmad,
       const std::shared_ptr<const DiscreteFunctionVariant>& p,
       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
@@ -46,7 +46,7 @@ class HypoelasticSolverHandler
                  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>& sigma,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
                  const std::shared_ptr<const ItemValueVariant>& ur,
                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const = 0;
@@ -61,7 +61,7 @@ class HypoelasticSolverHandler
           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>& sigma,
+          const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
           const std::shared_ptr<const DiscreteFunctionVariant>& aL,
           const std::shared_ptr<const DiscreteFunctionVariant>& aT,
           const std::shared_ptr<const DiscreteFunctionVariant>& p,