diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 0df2b95344a231490914a62adea82a55eb7dba32..eeffc83f5486b0dcfe09480e3c300ff91f268532 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -440,6 +440,63 @@ LocalFSIModule::LocalFSIModule()
 
                               ));
 
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 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,
+                                 const size_t& q)
+                                -> 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>, 
+                                  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 Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, q, 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, p_inf1, gamma2, p_inf2);
+                              }
+
+                              ));
+
+
   this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook2_solver_order2",
                             std::function(
 
@@ -493,6 +550,62 @@ LocalFSIModule::LocalFSIModule()
 
                               ));
 
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook2_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 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,
+                                 const size_t& q)
+                                -> 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>, 
+                                  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 Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt1, q, 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, p_inf1, gamma2, p_inf2);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
                             std::function(
 
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index 2def47e16c7ba64386df3566558da343cfef6353..dfcf6d6c8d200b23e5947f97629ab9ec2fdafbba 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -229,8 +229,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
              const NodeValuePerCell<const Rdxd>& Ajr,
              DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
              DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
-             DiscreteFunctionDPk<Dimension, const double> DPk_p,
-             CellValue<double> lambda, NodeValuePerCell<double> coef) const
+             DiscreteFunctionDPk<Dimension, const double> DPk_p) const
   {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
@@ -252,7 +251,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
         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];
-          const Rdxd sigma_r = sqrt(lambda[J]+(1-lambda[J])*coef(J,R))*DPk_S[J](xr[r]) - DPk_p[J](xr[r])*I;
+          const Rdxd sigma_r = DPk_S[J](xr[r]) - DPk_p[J](xr[r])*I;
 
           br += Ajr(J, R) * DPk_u[J](xr[r]) - sigma_r * Cjr(J, R);
         }
@@ -533,8 +532,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
               const NodeValue<const Rd>& ur,
               DiscreteFunctionDPk<Dimension, const Rd> DPk_uh,
               DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
-              DiscreteFunctionDPk<Dimension, const double> DPk_p,
-              CellValue<double> lambda, NodeValuePerCell<double> coef) const
+              DiscreteFunctionDPk<Dimension, const double> DPk_p) const
   {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
@@ -551,7 +549,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         for(size_t r = 0; r < cell_nodes.size(); ++r) {
-          Rdxd sigma_r = sqrt(lambda[j]+(1-lambda[j])*coef(j,r))*DPk_S[j](xr[cell_nodes[r]]) - DPk_p[j](xr[cell_nodes[r]])*I;
+          Rdxd sigma_r = DPk_S[j](xr[cell_nodes[r]]) - DPk_p[j](xr[cell_nodes[r]])*I;
           F(j,r) = -Ajr(j,r) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + sigma_r*Cjr(j,r);
         }
       });
@@ -640,408 +638,6 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     return std::make_tuple(map1, map2);
   }
 
-  void
-  _scalar_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const double>& f,
-                  DiscreteFunctionDPk<Dimension, double>& DPk_fh) 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);
-
-    const auto xj = mesh_data.xj();
-
-    CellValue<double> alph_J2{mesh.connectivity()};
-
-    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-      const double fj = f[cell_id];
-
-      double f_min = fj;
-      double f_max = fj;
-
-      const auto cell_stencil = stencil[cell_id];
-      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-        f_min = std::min(f_min, f[cell_stencil[i_cell]]);
-        f_max = std::max(f_max, f[cell_stencil[i_cell]]);
-      }
-
-      double f_bar_min = fj;
-      double f_bar_max = fj;
-
-      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]);
-
-        f_bar_min = std::min(f_bar_min, f_xk);
-        f_bar_max = std::max(f_bar_max, f_xk);
-      }
-
-      const double eps = 1E-14;
-      double coef1     = 1;
-      if (std::abs(f_bar_max - fj) > eps) {
-        coef1 = (f_max - fj) / ((f_bar_max - fj));
-      }
-
-      double coef2 = 1.;
-      if (std::abs(f_bar_min - fj) > eps) {
-        coef2 = (f_min - fj) / ((f_bar_min - fj));
-      }
-
-      const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
-      alph_J2[cell_id] = lambda;
-
-      auto coefficients = DPk_fh.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;
-      }
-    });
-  }
-
-  void
-  _vector_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const Rd>& f,
-                  DiscreteFunctionDPk<Dimension, Rd> DPk_fh,
-                  const CellValue<const Rdxd>& grad_u) 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);
-    
-    const auto xj = mesh_data.xj();
-
-    CellValue<Rdxd> sym_grad_u{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId j){
-        sym_grad_u[j] = grad_u[j] + transpose(grad_u[j]);
-    });
-
-    CellValue<SmallArray<double>> eigunvalues{mesh.connectivity()};
-    CellValue<std::vector<SmallVector<double>>> eigunvectors{mesh.connectivity()};
-    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
-      EigenvalueSolver{}.computeForSymmetricMatrix(sym_grad_u[j],eigunvalues[j],eigunvectors[j]);
-    }
-
-    CellValue<Rd> n{mesh.connectivity()};
-    CellValue<Rd> t{mesh.connectivity()};
-    CellValue<Rd> l{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-        n[cell_id] = zero;
-        t[cell_id] = zero;
-        l[cell_id] = zero;
-        if(frobeniusNorm(sym_grad_u[cell_id]) > 1E-5){
-        if constexpr(Dimension == 2){
-          Rd nj;
-          Rd tj;
-          for(size_t i = 0; i < Dimension; ++i){
-            nj[i] = eigunvectors[cell_id][0][i];
-            tj[i] = eigunvectors[cell_id][1][i];
-          }
-          n[cell_id] = nj;
-          t[cell_id] = tj;
-          if(l2Norm(nj) != 0){
-            n[cell_id] = (1/l2Norm(nj))*n[cell_id];
-          }
-          if(l2Norm(tj) != 0){
-            t[cell_id] = (1/l2Norm(tj))*t[cell_id];
-          }
-        }
-        else{
-          if constexpr(Dimension == 3){
-            Rd nj;
-            Rd tj;
-            Rd lj;
-            for(size_t i = 0; i < Dimension; ++i){
-              nj[i] = eigunvectors[cell_id][0][i];
-              tj[i] = eigunvectors[cell_id][1][i];
-              lj[i] = eigunvectors[cell_id][2][i];
-            }
-            n[cell_id] = nj;
-            t[cell_id] = tj;
-            l[cell_id] = lj;
-            if(l2Norm(nj) != 0){
-            n[cell_id] = (1/l2Norm(nj))*n[cell_id];
-            }
-            if(l2Norm(tj) != 0){
-              t[cell_id] = (1/l2Norm(tj))*t[cell_id];
-            }
-            if(l2Norm(lj) != 0){
-              l[cell_id] = (1/l2Norm(lj))*l[cell_id];
-            }
-          }
-        }
-        }
-    });
-
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-        const double fn = dot(f[cell_id],n[cell_id]);
-        const double ft = dot(f[cell_id],t[cell_id]);
-        const double fl = dot(f[cell_id],l[cell_id]);
-
-        double fn_min = fn;
-        double fn_max = fn;
-        double ft_min = ft;
-        double ft_max = ft;
-        double fl_min = fl;
-        double fl_max = fl;
-
-        const auto cell_stencil = stencil[cell_id];
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_id]);
-          fn_min            = std::min(fn_min,fn_k);
-          fn_max            = std::max(fn_max,fn_k);
-
-          const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_id]);
-          ft_min            = std::min(ft_min,ft_k);
-          ft_max            = std::max(ft_max,ft_k);
-
-          const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_id]);
-          fl_min            = std::min(fl_min,fl_k);
-          fl_max            = std::max(fl_max,fl_k);
-        }
-
-        double fn_bar_min = fn;
-        double fn_bar_max = fn;
-        double ft_bar_min = ft;
-        double ft_bar_max = ft;
-        double fl_bar_min = fl;
-        double fl_bar_max = fl;
-
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          const CellId cell_k_id = cell_stencil[i_cell];
-          const Rd f_xk          = DPk_fh[cell_id](xj[cell_k_id]);
-
-          const double fn_xk = dot(f_xk,n[cell_id]);
-          fn_bar_min         = std::min(fn_bar_min,fn_xk);
-          fn_bar_max         = std::max(fn_bar_max,fn_xk);
-
-          const double ft_xk = dot(f_xk,t[cell_id]);
-          ft_bar_min         = std::min(ft_bar_min,ft_xk);
-          ft_bar_max         = std::max(ft_bar_max,ft_xk);
-
-          const double fl_xk = dot(f_xk,l[cell_id]);
-          fl_bar_min         = std::min(fl_bar_min,fl_xk);
-          fl_bar_max         = std::max(fl_bar_max,fl_xk);
-        }
-
-        const double eps = 1E-14;
-        double coef1_n   = 1;
-        if(std::abs(fn_bar_max - fn) > eps){
-          coef1_n = (fn_max - fn) / (fn_bar_max - fn);
-        } 
-        double coef2_n = 1;
-        if(std::abs(fn_bar_min - fn) > eps){
-          coef2_n = (fn_min - fn) / (fn_bar_min - fn);
-        }
-        
-        double coef1_t = 1;
-        if(std::abs(ft_bar_max - ft) > eps){
-          coef1_t = (ft_max - ft) / (ft_bar_max - ft);
-        }
-        double coef2_t = 1;
-        if(std::abs(ft_bar_min - ft) > eps){
-          coef2_t = (ft_min - ft) / (ft_bar_min - ft);
-        }
-
-        double coef1_l = 1;
-        if(std::abs(fl_bar_max - fl) > eps){
-          coef1_l = (fl_max - fl) / (fl_bar_max - fl);
-        }
-        double coef2_l = 1;
-        if(std::abs(fl_bar_min - fl) > eps){
-          coef2_l = (fl_min - fl) / (fl_bar_min - fl);
-        }
-
-        const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
-        const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
-        const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
-
-        auto coefficients = DPk_fh.coefficients(cell_id);
-        coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id] 
-                        + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id] 
-                        + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
-        for(size_t i = 1; i < coefficients.size(); ++i){
-          coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id] 
-                          + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id] 
-                          + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
-        }
-    });
-  }
-
-    CellValue<double> 
-  _matrix_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const Rdxd>& S,
-                  DiscreteFunctionDPk<Dimension, Rdxd>& DPk_Sh) 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);
-
-    const auto xj = mesh_data.xj();
-
-    //A retirer + tard
-    const auto xr = mesh.xr();
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-    //
-
-    CellValue<double> lambda{mesh.connectivity()};
-    const DiscreteScalarFunction& inv2 = 0.5*trace(transpose(S)*S);
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-        const double inv2j = inv2[cell_id];
-        
-        double inv2_min = inv2j;
-        double inv2_max = inv2j;
-
-        const auto cell_stencil = stencil[cell_id];
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          inv2_min = std::min(inv2_min, inv2[cell_stencil[i_cell]]);
-          inv2_max = std::max(inv2_max, inv2[cell_stencil[i_cell]]);
-        }
-
-        double inv2_bar_min = inv2j;
-        double inv2_bar_max = inv2j;
-
-        const auto& cell_nodes = cell_to_node_matrix[cell_id];
-        for(size_t r = 0; r < cell_nodes.size(); ++r){
-          const double inv2_xr = 0.5*trace(transpose(DPk_Sh[cell_id](xr[cell_nodes[r]]))*DPk_Sh[cell_id](xr[cell_nodes[r]]));
-
-          inv2_bar_min = std::min(inv2_bar_min, inv2_xr);
-          inv2_bar_max = std::max(inv2_bar_max, inv2_xr);
-        }
-
-        const double eps =  1E-14;
-        double coef1     = 1.;
-        double coef2     = 1.;
-        if(std::abs(inv2_bar_max - inv2j) > eps){
-          coef1 = (inv2_max - inv2j) / (inv2_bar_max - inv2j);
-        } 
-        if(std::abs(inv2_bar_min - inv2j) > eps){
-          coef2 = (inv2_min - inv2j) / (inv2_bar_min - inv2j);
-        }
-
-        const double lambda_inv2 = std::max(0., std::min(1., std::min(coef1, coef2)));
-        lambda[cell_id] = lambda_inv2;
-    });
-    return lambda;
-  } 
-
-  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, 
-             const std::shared_ptr<const DiscreteFunctionVariant>, 
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>>
-  _apply_NeoHook(const DiscreteScalarFunction rho, 
-                  const DiscreteVectorFunction u, 
-                  const DiscreteScalarFunction E, 
-                  const DiscreteTensorFunction CG, 
-                  const size_t& chi_solid,
-                  const double& lambda,
-                  const double& mu,
-                  const double& gamma,
-                  const double& p_inf) const
-  {
-    const TinyMatrix<Dimension> I = identity;
-
-    const DiscreteScalarFunction& eps         = E - 0.5 * dot(u,u);
-    const DiscreteScalarFunction& p_d         = (1-chi_solid)*(gamma-1) * rho * eps - p_inf*gamma;
-    const DiscreteTensorFunction& sigma_d     = chi_solid * (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I) - p_d * I;
-    const DiscreteScalarFunction& aL_d        = sqrt(chi_solid * (lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
-    const DiscreteScalarFunction& aT_d        = sqrt(chi_solid * mu / rho);
-
-    return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
-            std::make_shared<DiscreteFunctionVariant>(aL_d),
-            std::make_shared<DiscreteFunctionVariant>(aT_d),
-            std::make_shared<DiscreteFunctionVariant>(p_d)};
-  }
-
-  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, 
-             const std::shared_ptr<const DiscreteFunctionVariant>, 
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>>
-  _apply_NeoHook2(const DiscreteScalarFunction rho, 
-                  const DiscreteVectorFunction u, 
-                  const DiscreteScalarFunction E, 
-                  const DiscreteTensorFunction CG, 
-                  const size_t& chi_solid,
-                  const double& mu,
-                  const double& gamma,
-                  const double& p_inf) const
-  {
-    const TinyMatrix<Dimension> I = identity;
-
-    const DiscreteScalarFunction& eps         = E - 0.5 * dot(u,u);
-    const DiscreteScalarFunction& p_d         = (gamma-1) * rho * eps - p_inf*gamma;
-    const DiscreteTensorFunction& sigma_d     = chi_solid * 2*mu/sqrt(det(CG))*(CG-1./3. * trace(CG)*I) - p_d * I;
-    const DiscreteScalarFunction& aL_d        = sqrt(chi_solid * (2 * mu) / rho + ((1-chi_solid)*gamma * p_d + chi_solid * (gamma*(p_d+p_inf))) / rho);
-    const DiscreteScalarFunction& aT_d        = sqrt(chi_solid * mu / rho);
-
-    return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
-            std::make_shared<DiscreteFunctionVariant>(aL_d),
-            std::make_shared<DiscreteFunctionVariant>(aT_d),
-            std::make_shared<DiscreteFunctionVariant>(p_d)};
-  }
-
-  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, 
-             const std::shared_ptr<const DiscreteFunctionVariant>, 
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>>
-  _apply_state_law(const size_t law,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& CG_v,
-                   const size_t& chi_solid,
-                   const double& lambda,
-                   const double& mu,
-                   const double& gamma,
-                   const double& p_inf) const 
-  {
-    const DiscreteScalarFunction& rho_d       = rho_v->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u_d         = u_v->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& E_d         = E_v->get<DiscreteScalarFunction>();
-    const DiscreteTensorFunction& CG_d        = CG_v->get<DiscreteTensorFunction>();
-
-    if(law == 1){
-      return _apply_NeoHook(rho_d,u_d,E_d,CG_d,chi_solid,lambda,mu,gamma,p_inf);
-    } else {
-      return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid,mu,gamma,p_inf);
-    }
-  }
-
-  NodeValuePerCell<double>
-  _compute_inv2_coefficient(const MeshType& mesh,
-                            const DiscreteTensorFunction& S,
-                            const DiscreteFunctionDPk<Dimension, Rdxd>& DPk_Sh) const
-  {
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-    const NodeValue<const Rd>& xr   = mesh.xr();
-
-    const DiscreteScalarFunction& J2 = 0.5 * trace(transpose(S)*S);
-    NodeValuePerCell<double> coef{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId j){
-        const auto& cell_nodes = cell_to_node_matrix[j];
-        for(size_t r = 0; r < cell_nodes.size(); ++r){
-          const double J2_p = 0.5*trace(transpose(DPk_Sh[j](xr[cell_nodes[r]]))*DPk_Sh[j](xr[cell_nodes[r]]));
-          if(J2_p != 0.){
-            coef(j,r) = J2[j]/J2_p;
-          }else{
-          coef(j,r) = 0;
-          }
-        }
-    });
-
-    return coef;
-  }
-
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>,
              const std::shared_ptr<const SubItemValuePerItemVariant>,
@@ -1136,12 +732,9 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     DiscreteFunctionDPk<Dimension, double> p1_lim   = copy(DPk_ph1);
 
     const CellValue<const Rdxd> grad_u1 = this->_computeGradU(solver_type,rho1_v,aL1_v,aT1_v,u1_v,sigma1_v,bc_descriptor_list1);
-    this->_vector_limiter(mesh1,u1,u1_lim,grad_u1);
-    CellValue<double> lambda1 = this->_matrix_limiter(mesh1, S1, S1_lim);
-    this->_scalar_limiter(mesh1, p1, p1_lim);
-
-    auto copy_DPk_Sh1 = copy(DPk_Sh1);
-    const NodeValuePerCell<double> coef1 = _compute_inv2_coefficient(mesh1, S1, copy_DPk_Sh1);
+    vector_limiter(mesh1,u1,u1_lim,grad_u1);
+    matrix_limiter(mesh1, S1, S1_lim);
+    scalar_limiter(mesh1, p1, p1_lim);
 
     auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v, S2_v, p2_v);
     auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
@@ -1153,20 +746,17 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     DiscreteFunctionDPk<Dimension, double> p2_lim   = copy(DPk_ph2);
 
     const CellValue<const Rdxd> grad_u2 = this->_computeGradU(solver_type,rho2_v,aL2_v,aT2_v,u2_v,sigma2_v,bc_descriptor_list2);
-    this->_vector_limiter(mesh2,u2,u2_lim,grad_u2);
-    CellValue<double> lambda2 = this->_matrix_limiter(mesh2, S2, S2_lim);
-    this->_scalar_limiter(mesh2, p2, p2_lim);
-
-    auto copy_DPk_Sh2 = copy(DPk_Sh2);
-    const NodeValuePerCell<double> coef2 = _compute_inv2_coefficient(mesh2, S2, copy_DPk_Sh2);
+    vector_limiter(mesh2,u2,u2_lim,grad_u2);
+    matrix_limiter(mesh2, S2, S2_lim);
+    scalar_limiter(mesh2, p2, p2_lim);
 
     NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
     NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
 
     NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
-    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1_lim, S1_lim, p1_lim, lambda1, coef1);
+    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1_lim, S1_lim, p1_lim);
     NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
-    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, S2_lim, p2_lim, lambda2, coef2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, S2_lim, p2_lim);
 
     const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
     const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
@@ -1183,8 +773,8 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
     this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
 
-    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, S1_lim, p1_lim, lambda1, coef1);
-    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, S2_lim, p2_lim, lambda2, coef2);
+    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, S1_lim, p1_lim);
+    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, S2_lim, p2_lim);
 
     NodeValue<Rd> CR_ur         = copy(ur2);
     NodeValuePerCell<Rd> CR_Fjr = copy(Fjr2);
@@ -1352,17 +942,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     DiscreteFunctionDPk<Dimension, double> p_lim   = copy(DPk_ph);
 
     const CellValue<const Rdxd> grad_u = this->_computeGradU(solver_type,rho_v,aL_v,aT_v,u_v,sigma_v,bc_descriptor_list);
-    this->_vector_limiter(mesh,u,u_lim,grad_u);
-    CellValue<double> lambda = this->_matrix_limiter(mesh, S, S_lim);
-    this->_scalar_limiter(mesh, p, p_lim);
-
-    auto copy_DPk_Sh = copy(DPk_Sh);
-    const NodeValuePerCell<double> coef = _compute_inv2_coefficient(mesh, S, copy_DPk_Sh);
+    vector_limiter(mesh,u,u_lim,grad_u);
+    matrix_limiter(mesh, S, S_lim);
+    scalar_limiter(mesh, p, p_lim);
 
     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, S_lim, p_lim, lambda, coef);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, S_lim, p_lim);
 
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -1371,7 +958,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     synchronize(br);
 
     NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim, lambda, coef);
+    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim);
 
     this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
 
@@ -1864,6 +1451,71 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     return {sub_m, sub_rho, sub_u, sub_E, sub_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>>
+  compute_sub_iterations(const SolverType& solver_type,
+                         const size_t& law,
+                         const double& dt,
+                         const size_t& q,
+                         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::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                         NodeValue<Rd> CR_ur,
+                         NodeValuePerCell<Rd> CR_Fjr,
+                         std::vector<NodeId> map2,
+                         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, 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 sub_dt = dt / q;
+    double sum_dt = 0.;
+    for(size_t i = 0; i < q; ++i){
+
+      if( i == q - 1){
+        sub_dt = dt - sum_dt;
+      }
+      
+      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>();
+      const std::shared_ptr<const DiscreteFunctionVariant>& c =
+        std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+
+      auto [sub_ur2, sub_Fjr2]       = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+      auto [sub_ur2_o1, sub_Fjr2_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+
+      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;
+
+    }
+
+    return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
+  }
+
   std::tuple<std::shared_ptr<const MeshVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1932,101 +1584,149 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     return {sub_m, sub_rho, sub_u, sub_E, sub_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>,
-  //           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 double& dt1,
-  //      const size_t& q,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-  //      const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-  //      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
-  //      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
-  //      const double& mu,
-  //      const double& lambda) const
-  //{
-  //  std::shared_ptr new_m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
-  //  std::shared_ptr m1     = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
-//
-  //  std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
-  //  std::shared_ptr<const DiscreteFunctionVariant> new_u2   = u2;
-  //  std::shared_ptr<const DiscreteFunctionVariant> new_E2   = E2;
-  //  std::shared_ptr<const DiscreteFunctionVariant> new_CG2  = CG2;
-//
-  //  auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
-//
-  //  auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] =
-  //    compute_fluxes(solver_type, rho1, aL1, aT1, u1, sigma1, bc_descriptor_list1, rho2, aL2, aT2, u2, sigma2,
-  //                   bc_descriptor_list2, map1, map2);
-  //  const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, Fjr1);
-//
-  //  double dt2                    = dt1 / q;
-  //  const double gamma            = 1.4;
-  //  const TinyMatrix<Dimension> I = identity;
-  //  double sum_dt                 = 0;
-//
-  //  size_t fluid = 0;
-  //  if ((mu == 0) and (lambda == 0)) {
-  //    fluid = 1;
-  //  }
-//
-  //  for (size_t i = 0; i < q; i++) {
-  //    if (i == q - 1) {
-  //      dt2 = dt1 - sum_dt;
-  //    }
-//
-  //    const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
-  //    const DiscreteVectorFunction& u_d   = new_u2->get<DiscreteVectorFunction>();
-  //    const DiscreteScalarFunction& E_d   = new_E2->get<DiscreteScalarFunction>();
-  //    const DiscreteScalarFunction& eps   = E_d - 0.5 * dot(u_d, u_d);
-  //    const DiscreteTensorFunction& CG_d  = new_CG2->get<DiscreteTensorFunction>();
-  //    const DiscreteScalarFunction& p     = fluid * (gamma - 1) * rho_d * eps;
-  //    const DiscreteTensorFunction& sigma_d =
-  //      mu / sqrt(det(CG_d)) * (CG_d - I) + lambda / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I - p * I;
-  //    const DiscreteScalarFunction& aL_d = sqrt((lambda + 2 * mu) / rho_d + gamma * p / rho_d);
-//
-  //    const DiscreteScalarFunction& aT_d = sqrt(mu / rho_d);
-//
-  //    const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 =
-  //      std::make_shared<const DiscreteFunctionVariant>(sigma_d);
-  //    const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 =
-  //      std::make_shared<const DiscreteFunctionVariant>(aL_d);
-  //    const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 =
-  //      std::make_shared<const DiscreteFunctionVariant>(aT_d);
-//
-  //    auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, new_rho2, new_aL2, new_aT2, new_u2, new_sigma2,
-  //                                      bc_descriptor_list2, CR_ur, CR_Fjr, map2);
-//
-  //    std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
-  //      apply_fluxes(dt2, new_rho2, new_u2, new_E2, new_CG2, sub_ur2, sub_Fjr2);
-//
-  //    sum_dt += dt2;
-  //  }
-//
-  //  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);
-//
-  //  return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
-  //}
+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>>
+compute_sub_iterations_midpoint_method(const SolverType& solver_type,
+                              const size_t& law,
+                              const double& dt,
+                              const size_t q,
+                              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::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                              NodeValue<Rd> CR_ur,
+                              NodeValuePerCell<Rd> CR_Fjr,
+                              std::vector<NodeId> map2,
+                              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;
+  std::shared_ptr<const DiscreteFunctionVariant> sub_u   = u;
+  std::shared_ptr<const DiscreteFunctionVariant> sub_E   = E;
+  std::shared_ptr<const DiscreteFunctionVariant> sub_CG  = CG;
+
+  double sub_dt = dt / q;
+  double sum_dt = 0;
+  for(size_t i = 0; i < q; ++i){
+
+    if( i == q - 1){
+      sub_dt = dt - sum_dt;
+    }
+
+    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>();
+    const std::shared_ptr<const DiscreteFunctionVariant>& c = std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+
+    auto [sub_ur, sub_Fjr]       = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+    auto [sub_ur_o1, sub_Fjr_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+
+    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 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);
+
+    std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = order2_apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_rho_app, sub_CG_app, sub_ur_app, sub_ur_o1_app, sub_Fjr_app, sub_Fjr_o1_app); 
+
+    sum_dt += sub_dt;
+  }
+
+  return {sub_m, sub_rho, sub_u, sub_E, sub_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>,
+             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 size_t& law,
+        const double& dt1,
+        const size_t& q,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+        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 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);
+
+    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 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,dt1/2.,q,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,dt1,q,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<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};
+  }
 
   std::tuple<std::shared_ptr<const MeshVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
@@ -2101,7 +1801,6 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     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};
-
   }
 
   Order2LocalDtHyperelasticSolver()                            = default;
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.hpp b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
index 8712feaf2d3e3bead44f2a5b83d61440ce5ef61c..633a87ccc2e93b074d4abee64b6d2b323d10e44d 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -26,37 +26,46 @@ class Order2LocalDtHyperelasticSolverHandler
  private:
   struct IOrder2LocalDtHyperelasticSolver
   {
-    //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>,
-    //                   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 double& dt1,
-    //      const size_t& q,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-    //      const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-    //      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
-    //      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
-    //      const double& mu,
-    //      const double& lambda) 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>,
+                       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 size_t& law,
+        const double& dt1,
+        const size_t& q,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+        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 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;
 
     virtual std::tuple<std::shared_ptr<const MeshVariant>,
                        std::shared_ptr<const DiscreteFunctionVariant>,