diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 40808c6ab3165f00540025735a7bba6863886ea3..78ad84a0414adaa3d6f68cc303611d4430477974 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -611,6 +611,7 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& aT,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
                                  const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -618,10 +619,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, sigmad, p})}
+                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigmad, p, mu})}
                                   .solver()
                                   .apply(HypoelasticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, sigmad, aL, aT,
-                                         p, bc_descriptor_list);
+                                         p, mu, bc_descriptor_list);
                               }
 
                               ));
@@ -658,6 +659,7 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& aT,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
                                  const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -665,10 +667,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, sigmad, p})}
+                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigmad, mu, p})}
                                   .solver()
                                   .apply(HypoelasticSolverHandler::SolverType::Glace, dt, rho, u, E, sigmad, aL, aT, p,
-                                         bc_descriptor_list);
+                                         mu, bc_descriptor_list);
                               }
 
                               ));
@@ -680,7 +682,7 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& u,        //
                                  const std::shared_ptr<const DiscreteFunctionVariant>& E,        //
                                  const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,   //
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& mu,       //
                                  const std::shared_ptr<const ItemValueVariant>& ur,              //
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
                                  const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -688,9 +690,9 @@ SchemeModule::SchemeModule()
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, sigmad, p})}   //
+                                return HypoelasticSolverHandler{getCommonMesh({rho, u, E, sigmad, mu})}   //
                                   .solver()
-                                  .apply_fluxes(dt, rho, u, E, sigmad, p, ur, Fjr);
+                                  .apply_fluxes(dt, rho, u, E, sigmad, mu, ur, Fjr);
                               }
 
                               ));
diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp
index f39f214f0915b5cfc0634ee7310c5a61be01f1cb..0bea644c94677c55a2c33fe7cda0d1175758eee7 100644
--- a/src/scheme/HypoelasticSolver.cpp
+++ b/src/scheme/HypoelasticSolver.cpp
@@ -434,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, sigmad);
+    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -451,8 +451,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
                const DiscreteVectorFunction& u,
                const DiscreteScalarFunction& E,
                const DiscreteTensorFunction& sigmad,
-               // p probably not needed here
-               const DiscreteScalarFunction& p,
+               const DiscreteScalarFunction& mu,
                const NodeValue<const Rd>& ur,
                const NodeValuePerCell<const Rd>& Fjr) const
   {
@@ -491,14 +490,16 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
           energy_fluxes += dot(Fjr(j, R), ur[r]);
         }
         // TODO
-        // const Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
+        const Rdxd I            = identity;
+        const Rdxd D            = 0.5 * (gradv + transpose(gradv));
+        const Rdxd Dd           = D - 1 / 3 * trace(D) * I;
+        const Rdxd gradva       = 0.5 * (gradv - transpose(gradv));
+        const Rdxd jaumann      = sigmad[j] * gradva - gradva * sigmad[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_CG[j] += dt_over_Vj * cauchy_green_fluxes;
-        // new_CG[j] += transpose(new_CG[j]);
-        // new_CG[j] *= 0.5;
+        new_sigmad[j] += dt_over_Vj * (2 * mu[j] * Dd - jaumann);
       });
 
     CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
@@ -522,7 +523,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
                const std::shared_ptr<const DiscreteFunctionVariant>& u,
                const std::shared_ptr<const DiscreteFunctionVariant>& E,
                const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
-               const std::shared_ptr<const DiscreteFunctionVariant>& p,
+               const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                const std::shared_ptr<const ItemValueVariant>& ur,
                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
   {
@@ -541,7 +542,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
                               u->get<DiscreteVectorFunction>(),         //
                               E->get<DiscreteScalarFunction>(),         //
                               sigmad->get<DiscreteTensorFunction>(),    //
-                              p->get<DiscreteScalarFunction>(),         //
+                              mu->get<DiscreteScalarFunction>(),        //
                               ur->get<NodeValue<const Rd>>(),           //
                               Fjr->get<NodeValuePerCell<const Rd>>());
   }
@@ -560,10 +561,11 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
         const std::shared_ptr<const DiscreteFunctionVariant>& aL,
         const std::shared_ptr<const DiscreteFunctionVariant>& aT,
         const std::shared_ptr<const DiscreteFunctionVariant>& p,
+        const std::shared_ptr<const DiscreteFunctionVariant>& mu,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     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);
+    return apply_fluxes(dt, rho, u, E, sigmad, mu, ur, Fjr);
   }
 
   HypoelasticSolver()                    = default;
diff --git a/src/scheme/HypoelasticSolver.hpp b/src/scheme/HypoelasticSolver.hpp
index 5c21a83ac60adf5bb98c22a27e14ee62242e3c6e..f0239743f063d5752de0b30670ae72f07f722c59 100644
--- a/src/scheme/HypoelasticSolver.hpp
+++ b/src/scheme/HypoelasticSolver.hpp
@@ -47,7 +47,7 @@ class HypoelasticSolverHandler
                  const std::shared_ptr<const DiscreteFunctionVariant>& u,
                  const std::shared_ptr<const DiscreteFunctionVariant>& E,
                  const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                  const std::shared_ptr<const ItemValueVariant>& ur,
                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const = 0;
 
@@ -65,6 +65,7 @@ class HypoelasticSolverHandler
           const std::shared_ptr<const DiscreteFunctionVariant>& aL,
           const std::shared_ptr<const DiscreteFunctionVariant>& aT,
           const std::shared_ptr<const DiscreteFunctionVariant>& p,
+          const std::shared_ptr<const DiscreteFunctionVariant>& mu,
           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
     IHypoelasticSolver()                                = default;