diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 9915d4e8d9417bd5948c41cc1b32fddc154a97e4..0df2b95344a231490914a62adea82a55eb7dba32 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -410,9 +410,15 @@ LocalFSIModule::LocalFSIModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p2,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list2,
-                                 const double& lambda1, const double& mu1,
-                                 const double& lambda2, const double& mu2, 
-                                 const double& gamma1,  const double& gamma2 , const double& dt1)
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& lambda1, 
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& mu1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& lambda2, 
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& mu2, 
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& gamma1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2, 
+                                 const double& dt1)
                                 -> std::tuple<
                                   std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
                                   std::shared_ptr<const DiscreteFunctionVariant>,
@@ -429,7 +435,7 @@ LocalFSIModule::LocalFSIModule()
                                   .solver()
                                   .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
                                          u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, p1, p2, bc_descriptor_list1,
-                                         bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, 0, gamma2, 0);
+                                         bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, p_inf1, gamma2, p_inf2);
                               }
 
                               ));
@@ -457,10 +463,15 @@ LocalFSIModule::LocalFSIModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p2,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list2,
-                                 const double& lambda1, const double& mu1,
-                                 const double& lambda2, const double& mu2, 
-                                 const double& gamma1,  const double& p_inf1,
-                                 const double& gamma2,  const double& p_inf2, const double& dt1)
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& lambda1, 
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& mu1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& mu2, 
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& gamma1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2, 
+                                 const double& dt1)
                                 -> std::tuple<
                                   std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
                                   std::shared_ptr<const DiscreteFunctionVariant>,
diff --git a/src/scheme/LocalDtUtils.hpp b/src/scheme/LocalDtUtils.hpp
index cfd8ab7ff3264e1d1e888d5e1fe870583b3e4d7d..fc59315ac326ac7d65a7bf7d6d931252de465447 100644
--- a/src/scheme/LocalDtUtils.hpp
+++ b/src/scheme/LocalDtUtils.hpp
@@ -67,7 +67,7 @@ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
       aT_d = sqrt(mu_d / rho_d);
     } else if (law == 2){
       p_d = (gamma_d - 1) * rho_d * eps_d - gamma_d * p_inf_d;
-      sigma_d = mu_d / sqrt(det(CG_d)) * (CG_d - 1. / 3. * trace(CG_d) * I) - p_d * I;
+      sigma_d = 2*mu_d / sqrt(det(CG_d)) * (CG_d - 1. / 3. * trace(CG_d) * I) - p_d * I;
       aL_d = sqrt((2 * mu_d) / rho_d + (gamma_d * (p_d + p_inf_d)) / rho_d);
       aT_d = sqrt(mu_d / rho_d);
     } else {
@@ -97,8 +97,8 @@ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
                   std::vector<NodeId>& map1,
                   std::vector<NodeId>& map2)
   {
-    const MeshType mesh1 = *i_mesh1->get<MeshType>();
-    const MeshType mesh2 = *i_mesh2->get<MeshType>();
+    const MeshType& mesh1 = *i_mesh1->get<MeshType>();
+    const MeshType& mesh2 = *i_mesh2->get<MeshType>();
 
     NodeValue<TinyVector<MeshType::Dimension>> xr1     = copy(mesh1.xr());
     NodeValue<TinyVector<MeshType::Dimension>> new_xr2 = copy(mesh2.xr());
@@ -112,13 +112,14 @@ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
     }
 
     std::shared_ptr<const MeshType> new_mesh2 = std::make_shared<MeshType>(mesh2.shared_connectivity(), new_xr2);
+    const std::shared_ptr<const MeshVariant>& mesh_copy = std::make_shared<MeshVariant>(new_mesh2);
 
-    const std::shared_ptr<const DiscreteFunctionVariant>& new_rho = shallowCopy(new_mesh2, rho);
-    const std::shared_ptr<const DiscreteFunctionVariant>& new_u   = shallowCopy(new_mesh2, u);
-    const std::shared_ptr<const DiscreteFunctionVariant>& new_E   = shallowCopy(new_mesh2, E);
-    const std::shared_ptr<const DiscreteFunctionVariant>& new_CG  = shallowCopy(new_mesh2, CG);
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_rho = shallowCopy(mesh_copy, rho);
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_u   = shallowCopy(mesh_copy, u);
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_E   = shallowCopy(mesh_copy, E);
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_CG  = shallowCopy(mesh_copy, CG);
 
-    return std::make_tuple(new_mesh2,new_rho,new_u,new_E,new_CG);
+    return std::make_tuple(mesh_copy,new_rho,new_u,new_E,new_CG);
   }
 
 #endif //LOCAL_DT_UTILS_HPP
\ No newline at end of file
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index 0dfe808d2171ee0bc61cff2b929f8167fafb664b..2def47e16c7ba64386df3566558da343cfef6353 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -24,6 +24,8 @@
 #include <scheme/HyperelasticSolver.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/LocalDtUtils.hpp>
+#include <scheme/Order2Limiters.hpp>
 #include <scheme/PolynomialReconstruction.hpp>
 #include <scheme/PolynomialReconstructionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -1799,64 +1801,6 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                                      Fjr_o1->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>>
-  mesh_correction(const MeshType& mesh1,
-                  const MeshType& mesh2,
-                  const DiscreteScalarFunction& rho,
-                  const DiscreteVectorFunction& u,
-                  const DiscreteScalarFunction& E,
-                  const DiscreteTensorFunction& CG,
-                  std::vector<NodeId>& map1,
-                  std::vector<NodeId>& map2) const
-  {
-    NodeValue<Rd> xr1     = copy(mesh1.xr());
-    NodeValue<Rd> new_xr2 = copy(mesh2.xr());
-
-    size_t n = map1.size();
-
-    for (size_t i = 0; i < n; i++) {
-      for (size_t j = 0; j < Dimension; j++) {
-        new_xr2[map2[i]][j] = xr1[map1[i]][j];
-      }
-    }
-
-    std::shared_ptr<const MeshType> new_mesh2 = std::make_shared<MeshType>(mesh2.shared_connectivity(), new_xr2);
-
-    CellValue<double> new_rho = copy(rho.cellValues());
-    CellValue<Rd> new_u       = copy(u.cellValues());
-    CellValue<double> new_E   = copy(E.cellValues());
-    CellValue<Rdxd> new_CG    = copy(CG.cellValues());
-
-    return {std::make_shared<MeshVariant>(new_mesh2),
-            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_rho)),
-            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh2, new_u)),
-            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_E)),
-            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh2, new_CG))};
-  }
-
-  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>>
-  mesh_correction(const std::shared_ptr<const MeshVariant>& i_mesh1,
-                  const std::shared_ptr<const MeshVariant>& i_mesh2,
-                  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,
-                  std::vector<NodeId>& map1,
-                  std::vector<NodeId>& map2) const
-  {
-    return this->mesh_correction(*i_mesh1->get<MeshType>(), *i_mesh2->get<MeshType>(),
-                                 rho->get<DiscreteScalarFunction>(), u->get<DiscreteVectorFunction>(),
-                                 E->get<DiscreteScalarFunction>(), CG->get<DiscreteTensorFunction>(), map1, map2);
-  }
-
   std::tuple<std::shared_ptr<const MeshVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1874,21 +1818,25 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                          NodeValue<Rd> CR_ur,
                          NodeValuePerCell<Rd> CR_Fjr,
                          std::vector<NodeId> map2,
-                         const size_t& chi_solid,
-                         const double& lambda,
-                         const double& mu,
-                         const double& gamma,
-                         const double& p_inf) const
+                         const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
   {
-    std::shared_ptr sub_m = getCommonMesh({rho, u});
+    std::shared_ptr sub_m = getCommonMesh({rho, u, E, CG});
     std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
     std::shared_ptr<const DiscreteFunctionVariant> sub_u   = u;
     std::shared_ptr<const DiscreteFunctionVariant> sub_E   = E;
     std::shared_ptr<const DiscreteFunctionVariant> sub_CG  = CG;
 
+    std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = lambda;
+    std::shared_ptr<const DiscreteFunctionVariant> sub_mu     = mu;
+    std::shared_ptr<const DiscreteFunctionVariant> sub_gamma  = gamma;
+    std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf  = p_inf;
+
     double sum_dt = 0;
     while(sum_dt < dt){
-      const auto& [sigma, aL, aT, p] = _apply_state_law(law,sub_rho,sub_u,sub_E,sub_CG,chi_solid,lambda,mu,gamma,p_inf);
+      const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);
 
       const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
       const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
@@ -1905,6 +1853,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
       std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_ur2, sub_ur2_o1, sub_Fjr2, sub_Fjr2_o1);
 
+      sub_lambda = shallowCopy(sub_m,sub_lambda);
+      sub_mu     = shallowCopy(sub_m,sub_mu);
+      sub_gamma  = shallowCopy(sub_m,sub_gamma);
+      sub_p_inf  = shallowCopy(sub_m,sub_p_inf);
+
       sum_dt += sub_dt;
     }
 
@@ -1928,11 +1881,10 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                                          NodeValue<Rd> CR_ur,
                                          NodeValuePerCell<Rd> CR_Fjr,
                                          std::vector<NodeId> map2,
-                                         const size_t& chi_solid,
-                                         const double& lambda,
-                                         const double& mu,
-                                         const double& gamma,
-                                         const double& p_inf) const
+                                         const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
+                                         const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+                                         const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+                                         const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
   {
     std::shared_ptr sub_m = getCommonMesh({rho, u});
     std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
@@ -1942,7 +1894,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
     double sum_dt = 0;
     while(sum_dt < dt){
-      const auto& [sigma, aL, aT, p] = _apply_state_law(law,sub_rho,sub_u,sub_E,sub_CG,chi_solid,lambda,mu,gamma,p_inf);
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = shallowCopy(sub_m,lambda); 
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_mu     = shallowCopy(sub_m,mu); 
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_gamma  = shallowCopy(sub_m,gamma);
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf  = shallowCopy(sub_m,p_inf);       
+      const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);
 
       const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
       const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
@@ -1959,7 +1915,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
       auto [sub_m_app, sub_rho_app, sub_u_app, sub_E_app, sub_CG_app] = apply_fluxes(sub_dt/2., sub_rho, sub_u, sub_E, sub_CG, sub_ur, sub_ur_o1, sub_Fjr, sub_Fjr_o1);
 
-      const auto& [sigma_app, aL_app, aT_app, p_app] = _apply_state_law(law,sub_rho_app,sub_u_app,sub_E_app,sub_CG_app,chi_solid,lambda,mu,gamma,p_inf);
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_lambda_app = shallowCopy(sub_m_app,lambda); 
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_mu_app     = shallowCopy(sub_m_app,mu); 
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_gamma_app  = shallowCopy(sub_m_app,gamma);
+      const std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf_app  = shallowCopy(sub_m_app,p_inf);
+      const auto& [sigma_app, p_app, aL_app, aT_app] = apply_state_law<MeshType>(law,sub_rho_app,sub_u_app,sub_E_app,sub_CG_app,sub_lambda_app,sub_mu_app,sub_gamma_app,sub_p_inf_app);
 
       auto [sub_ur_app, sub_Fjr_app]       = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, p_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);
       auto [sub_ur_o1_app, sub_Fjr_o1_app] = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);
@@ -2099,43 +2059,46 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
         const std::shared_ptr<const DiscreteFunctionVariant>& p2,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
-        const double& lambda_f,
-        const double& mu_f,
-        const double& lambda_s,
-        const double& mu_s,
-        const double& gamma_f,
-        const double& p_inf_f,
-        const double& gamma_s,
-        const double& p_inf_s) const
+        const std::shared_ptr<const DiscreteFunctionVariant>& lambda1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& mu1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& gamma1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2) const
   {
     std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
     std::shared_ptr m1     = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
 
     auto [map1, map2] = _computeMapping(m1,m2,bc_descriptor_list1,bc_descriptor_list2);
 
-    //Préparation au calcul des lois d'états
-    size_t chi_solid = 0;
-    if((mu_s!=0) or (lambda_s!=0)){
-      chi_solid = 1;
-    }
-    const size_t chi_fluid = 1-chi_solid;
-
     auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr]                   = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,p1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,p2,bc_descriptor_list2,map1,map2);
     auto [ur1_o1, Fjr1_o1, ur2_o1, Fjr2_o1, CR_ur_o1, CR_Fjr_o1] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,bc_descriptor_list2,map1,map2);
-
+    
     auto [m1_app, rho1_app, u1_app, E1_app, CG1_app]   = apply_fluxes(dt1/2., rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);
-    const auto& [sigma1_app, aL1_app, aT1_app, p1_app] = _apply_state_law(law, rho1_app, u1_app, E1_app, CG1_app, chi_fluid, lambda_f, mu_f, gamma_f, p_inf_f);
 
-    auto [m2_app, rho2_app, u2_app, E2_app, CG2_app]   = compute_sub_iterations(solver_type,law,0.03,dt1/2.,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
-    const auto& [sigma2_app, aL2_app, aT2_app, p2_app] = _apply_state_law(law, rho2_app, u2_app, E2_app, CG2_app, chi_solid, lambda_s, mu_s, gamma_s, p_inf_s);
+    const std::shared_ptr<const DiscreteFunctionVariant>& lambda1_app = shallowCopy(m1_app,lambda1);
+    const std::shared_ptr<const DiscreteFunctionVariant>& mu1_app = shallowCopy(m1_app,mu1);
+    const std::shared_ptr<const DiscreteFunctionVariant>& gamma1_app = shallowCopy(m1_app,gamma1);
+    const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1_app = shallowCopy(m1_app,p_inf1);
+    const auto& [sigma1_app, p1_app, aL1_app, aT1_app] = apply_state_law<MeshType>(law, rho1_app, u1_app, E1_app, CG1_app, lambda1_app, mu1_app, gamma1_app, p_inf1_app);
+
+    auto [m2_app, rho2_app, u2_app, E2_app, CG2_app]   = compute_sub_iterations(solver_type,law,0.4,dt1/2.,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,lambda2,mu2,gamma2,p_inf2);
+
+    const std::shared_ptr<const DiscreteFunctionVariant>& lambda2_app = shallowCopy(m2_app,lambda2);
+    const std::shared_ptr<const DiscreteFunctionVariant>& mu2_app = shallowCopy(m2_app,mu2);
+    const std::shared_ptr<const DiscreteFunctionVariant>& gamma2_app = shallowCopy(m2_app,gamma2);
+    const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2_app = shallowCopy(m2_app,p_inf2);
+    const auto& [sigma2_app, p2_app, aL2_app, aT2_app] = apply_state_law<MeshType>(law, rho2_app, u2_app, E2_app, CG2_app, lambda2_app, mu2_app, gamma2_app, p_inf2_app);
 
     auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app]                   = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,p1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,p2_app,bc_descriptor_list2,map1,map2);
     auto [ur1_o1_app, Fjr1_o1_app, ur2_o1_app, Fjr2_o1_app, CR_ur_o1_app, CR_Fjr_o1_app] = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,bc_descriptor_list2,map1,map2);
 
     const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = order2_apply_fluxes(dt1, rho1, u1, E1, CG1, rho1_app, CG1_app, ur1_app, ur1_o1_app, Fjr1_app, Fjr1_o1_app);
-    auto [new_m2, new_rho2, new_u2, new_E2, new_CG2]       = compute_sub_iterations_midpoint_method(solver_type,law,0.03,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+    auto [new_m2, new_rho2, new_u2, new_E2, new_CG2]       = compute_sub_iterations_midpoint_method(solver_type,law,0.4,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur_app,CR_Fjr_app,map2,lambda2,mu2,gamma2,p_inf2);
 
-    std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) = mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
+    std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) = mesh_correction<MeshType>(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
 
     return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
 
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.hpp b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
index b0682e12a0875449ed3ac4dbf762a56a8ca9288f..8712feaf2d3e3bead44f2a5b83d61440ce5ef61c 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -89,14 +89,14 @@ class Order2LocalDtHyperelasticSolverHandler
         const std::shared_ptr<const DiscreteFunctionVariant>& p2,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
-        const double& lambda_f,
-        const double& mu_f,
-        const double& lambda_s,
-        const double& mu_s,
-        const double& gamma_f,
-        const double& p_inf_f,
-        const double& gamma_s,
-        const double& p_inf_s) const = 0;
+        const std::shared_ptr<const DiscreteFunctionVariant>& lambda1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& mu1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& gamma1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2) const = 0;
 
     IOrder2LocalDtHyperelasticSolver()                                        = default;
     IOrder2LocalDtHyperelasticSolver(IOrder2LocalDtHyperelasticSolver&&)            = default;