From df485b076da2f6ff4ba724b4c88411d7318422af Mon Sep 17 00:00:00 2001
From: Alexandre Gangloff <alexandre.gangloff@cea.fr>
Date: Thu, 20 Feb 2025 16:36:42 +0100
Subject: [PATCH] Add space second order extension for couplig method

---
 src/language/modules/LocalFSIModule.cpp       |   12 +-
 src/scheme/Order2HyperelasticSolver.cpp       |   83 -
 .../Order2LocalDtHyperelasticSolver.cpp       | 1645 +++++++++++------
 .../Order2LocalDtHyperelasticSolver.hpp       |   88 +-
 4 files changed, 1084 insertions(+), 744 deletions(-)

diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 3ed298e6e..76d2bc03f 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -402,6 +402,7 @@ LocalFSIModule::LocalFSIModule()
                                  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,
@@ -411,6 +412,7 @@ LocalFSIModule::LocalFSIModule()
                                  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 double& lambda1, const double& mu1,
@@ -431,13 +433,13 @@ LocalFSIModule::LocalFSIModule()
                                                                           {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
                                   .solver()
                                   .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         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);
                               }
 
                               ));
 
-  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook2_solver_order2",
                             std::function(
 
                               [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -447,6 +449,7 @@ LocalFSIModule::LocalFSIModule()
                                  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,
@@ -456,6 +459,7 @@ LocalFSIModule::LocalFSIModule()
                                  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 double& lambda1, const double& mu1,
@@ -476,8 +480,8 @@ LocalFSIModule::LocalFSIModule()
                                                                         getCommonMesh(
                                                                           {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
                                   .solver()
-                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 2, 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, p_inf1, gamma2, p_inf2);
                               }
 
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index 96f1db1a0..6d0dd07db 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -326,72 +326,6 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
     this->_applyVelocityBC(bc_list, Ar, br);
   }
 
-  CellValue<Rdxd>
-  _computeEigenVectorsMatrix(const MeshType& mesh,
-                             const CellValue<Rdxd>& grad_v,
-                             const double& eps) const
-  {
-    CellValue<Rdxd> A = copy(grad_v);
-
-    CellValue<Rdxd> V{mesh.connectivity()};
-    const TinyMatrix<Dimension> I = identity;
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { 
-        V[j] = I; 
-
-        while(true){
-          double max_val = 0;
-          size_t p = 0;
-          size_t q = 1;
-          for(size_t i = 0; i < Dimension; ++i){
-            for(size_t k = i + 1; k < Dimension; ++k){
-              if(abs(A[j](i,k)) > max_val){
-                max_val = abs(A[j](i,k));
-                p = i;
-                q = k;
-              }
-            }
-          }
-          if(max_val < eps){
-            break;
-          }
-
-          const double& theta = 0.5*(A[j](q,q) - A[j](p,p)) / A[j](p, q);
-          double t;
-          if(theta >= 0){
-            t = 1/(theta  + sqrt(theta*theta + 1));
-          }else{
-            t = -1/(-theta + sqrt(theta*theta + 1));
-          }
-          const double& c = 1/sqrt(t*t + 1);
-          const double& s = t * c;
-
-          for(size_t i = 0; i < Dimension; ++i){
-            if((i != p) and (i!=q)){
-              const double& temp_p = c * A[j](i,p) - s * A[j](i,q);
-              const double& temp_q = q * A[j](i,p) + c * A[j](i,q);
-              A[j](i,p) = temp_p;
-              A[j](i,q) = temp_q;
-              A[j](p,i) = temp_p;
-              A[j](q,i) = temp_q;
-            }
-          }
-          A[j](p,p) = A[j](p,p) - t * A[j](p,q);
-          A[j](q,q) = A[j](q,q) + t * A[j](p,q);
-          A[j](p,q) = 0.;
-          A[j](q,p) = 0.;
-
-          for(size_t i = 0; i < Dimension; ++i){
-            V[j](i,p) = c * V[j](i,p) - s * V[j](i,q);
-            V[j](i,q) = s * V[j](i,p) + c * V[j](i,q);
-          }
-        }
-         
-    });
-
-    return V;
-  }
-
   CellValue<Rdxd>
   _computeGradU(const SolverType& solver_type,
                const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
@@ -674,8 +608,6 @@ void
     for(CellId j = 0; j < mesh.numberOfCells(); ++j){
       EigenvalueSolver{}.computeForSymmetricMatrix(sym_grad_u[j],eigunvalues[j],eigunvectors[j]);
     }
-    
-    //CellValue<Rdxd> V = _computeEigenVectorsMatrix(mesh,sym_grad_u,1E-4);
 
     CellValue<Rd> n{mesh.connectivity()};
     CellValue<Rd> t{mesh.connectivity()};
@@ -689,10 +621,7 @@ void
         if constexpr(Dimension == 2){
           Rd nj;
           Rd tj;
-          //Rdxd Vj = V[cell_id];
           for(size_t i = 0; i < Dimension; ++i){
-            // nj[i] = Vj(i,0);
-            // tj[i] = Vj(i,1);
             nj[i] = eigunvectors[cell_id][0][i];
             tj[i] = eigunvectors[cell_id][1][i];
           }
@@ -710,11 +639,7 @@ void
             Rd nj;
             Rd tj;
             Rd lj;
-        //  Rdxd Vj = V[cell_id];
             for(size_t i = 0; i < Dimension; ++i){
-        //      nj[i] = Vj(i,0);
-        //      tj[i] = Vj(i,1);
-        //      lj[i] = Vj(i,2);
               nj[i] = eigunvectors[cell_id][0][i];
               tj[i] = eigunvectors[cell_id][1][i];
               lj[i] = eigunvectors[cell_id][2][i];
@@ -1268,7 +1193,6 @@ void
 
     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);
-    //this->_vector_limiter(mesh, u, u_lim, DPk_ph);
     CellValue<double> lambda = this->_matrix_limiter(mesh, S, S_lim);
     this->_scalar_limiter(mesh, p, p_lim);
 
@@ -1706,12 +1630,6 @@ void
         const double& gamma_s,
         const double& p_inf_s) const
   {
-    //std::shared_ptr m_app                                  = getCommonMesh({rho, aL, aT, u, sigma});
-    //std::shared_ptr<const DiscreteFunctionVariant> rho_app = rho;
-    //std::shared_ptr<const DiscreteFunctionVariant> u_app   = u;
-    //std::shared_ptr<const DiscreteFunctionVariant> E_app   = E;
-    //std::shared_ptr<const DiscreteFunctionVariant> CG_app  = CG;
-
     auto [ur, Fjr]       = compute_fluxes(solver_type, rho, aL, aT, u, sigma, p, bc_descriptor_list);
     auto [ur_o1, Fjr_o1] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
 
@@ -1725,7 +1643,6 @@ void
     auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app, p_app, bc_descriptor_list);
     auto [ur_o1_app, Fjr_o1_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
     return order2_apply_fluxes(dt, rho, u, E, CG, rho_app, CG_app, ur_app, ur_o1_app, Fjr_app, Fjr_o1_app);
-    //return apply_fluxes(dt, rho, u, E, CG, ur, ur_o1, Fjr, Fjr_o1);
   }
 
   Order2HyperelasticSolver()                           = default;
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index a80ae83fe..a98d525a1 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -1,3 +1,4 @@
+#include <algebra/EigenvalueSolver.hpp>
 #include <scheme/Order2LocalDtHyperelasticSolver.hpp>
 
 #include <language/utils/InterpolateItemValue.hpp>
@@ -221,6 +222,78 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     return b;
   }
 
+  NodeValue<Rd>
+  _computeBr(const Mesh<Dimension>& mesh,
+             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
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+    const auto& xr                        = mesh.xr();
+    const TinyMatrix<Dimension> I = identity;
+
+    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+    NodeValue<Rd> b{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+        Rd br = zero;
+        for (size_t j = 0; j < node_to_cell.size(); ++j) {
+          const CellId J       = node_to_cell[j];
+          const unsigned int R = node_local_number_in_its_cells[j];
+          const Rdxd sigma_r = sqrt(lambda[J]+(1-lambda[J])*coef(J,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);
+        }
+
+        b[r] = br;
+      });
+
+    return b;
+  }
+
+    NodeValue<Rd>
+  _computeBr(const Mesh<Dimension>& mesh,
+             const NodeValuePerCell<const Rdxd>& Ajr,
+             const DiscreteVectorFunction& u,
+             const DiscreteTensorFunction& sigma) const
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+
+    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+    NodeValue<Rd> b{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+        Rd br = zero;
+        for (size_t j = 0; j < node_to_cell.size(); ++j) {
+          const CellId J       = node_to_cell[j];
+          const unsigned int R = node_local_number_in_its_cells[j];
+          br += Ajr(J, R) * u[J] - sigma[J] * Cjr(J, R);
+        }
+
+        b[r] = br;
+      });
+
+    return b;
+  }
+
   BoundaryConditionList
   _getBCList(const MeshType& mesh,
              const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
@@ -354,6 +427,63 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     this->_applyVelocityBC(bc_list, Ar, br);
   }
 
+  CellValue<Rdxd>
+  _computeGradU(const SolverType& solver_type,
+               const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+  {
+    std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+    if (not mesh_v) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperelastic solver expects P0 functions");
+    }
+
+    const MeshType& mesh                = *mesh_v->get<MeshType>();
+    const DiscreteScalarFunction& rho   = rho_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u     = u_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+
+    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, sigma);
+
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+
+    synchronize(Ar);
+    synchronize(br);
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+    const NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
+
+    CellValue<Rdxd> grad_u{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        auto cell_nodes = cell_to_node_matrix[j];
+
+        Rdxd gradv           = zero;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+        }
+        grad_u[j] = gradv;
+      });
+
+    return grad_u;
+  }
+
   NodeValue<Rd>
   _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
   {
@@ -395,6 +525,64 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     return F;
   }
 
+  NodeValuePerCell<Rd>
+  _computeFjr(const MeshType& mesh,
+              const NodeValuePerCell<const Rdxd>& Ajr,
+              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
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const NodeValue<const Rd>& xr        = mesh.xr();
+    const TinyMatrix<Dimension> I = identity;
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    NodeValuePerCell<Rd> F{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        for(size_t r = 0; r < cell_nodes.size(); ++r) {
+          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;
+          F(j,r) = -Ajr(j,r) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + sigma_r*Cjr(j,r);
+        }
+      });
+
+    return F;
+  }
+
+  NodeValuePerCell<Rd>
+  _computeFjr(const MeshType& mesh,
+              const NodeValuePerCell<const Rdxd>& Ajr,
+              const NodeValue<const Rd>& ur,
+              const DiscreteVectorFunction& u,
+              const DiscreteTensorFunction& sigma) const
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    NodeValuePerCell<Rd> F{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        for (size_t r = 0; r < cell_nodes.size(); ++r) {
+          F(j, r) = -Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + sigma[j] * Cjr(j, r);
+        }
+      });
+
+    return F;
+  }
+
   std::tuple<std::vector<NodeId>, std::vector<NodeId>>
   _computeMapping(std::shared_ptr<const MeshVariant>& i_mesh1,
                   std::shared_ptr<const MeshVariant>& i_mesh2,
@@ -450,10 +638,10 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     return std::make_tuple(map1, map2);
   }
 
-     void
-  _vector_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const Rd>& f,
-                  DiscreteFunctionDPk<Dimension, Rd>& DPk_fh) const
+  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;
@@ -462,58 +650,44 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
     const auto xj = mesh_data.xj();
 
+    CellValue<double> alph_J2{mesh.connectivity()};
+
     parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-      const Rd fj = f[cell_id];
+      const double fj = f[cell_id];
 
-      Rd f_min = fj;
-      Rd f_max = fj;
+      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){
-        for(size_t dim = 0; dim < Dimension; ++dim){
-          f_min[dim] = std::min(f_min[dim], f[cell_stencil[i_cell]][dim]);
-          f_max[dim] = std::max(f_max[dim], f[cell_stencil[i_cell]][dim]);
-        }
+        f_min = std::min(f_min, f[cell_stencil[i_cell]]);
+        f_max = std::max(f_max, f[cell_stencil[i_cell]]);
       }
 
-      Rd f_bar_min = fj;
-      Rd f_bar_max = fj;
+      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 Rd f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
+        const double f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
 
-        for(size_t dim = 0; dim < Dimension; ++dim){
-          f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]);
-          f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]);
-        }
+        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;
-      Rd coef1;
-      for(size_t dim = 0; dim < Dimension; ++dim){
-        coef1[dim] = 1;
-        if (std::abs(f_bar_max[dim] - fj[dim]) > eps) {
-          coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim]));
-        }
-      }
-
-      Rd coef2;
-      for(size_t dim = 0; dim < Dimension; ++dim){
-        coef2[dim] = 1;
-        if (std::abs(f_bar_min[dim] - fj[dim]) > eps) {
-          coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim]));
-        }
+      double coef1     = 1;
+      if (std::abs(f_bar_max - fj) > eps) {
+        coef1 = (f_max - fj) / ((f_bar_max - fj));
       }
 
-      double min_coef1 = coef1[0];
-      double min_coef2 = coef2[0];
-      for(size_t dim = 1; dim < Dimension; ++dim){
-        min_coef1 = std::min(min_coef1,coef1[dim]);
-        min_coef2 = std::min(min_coef2,coef2[dim]);
+      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(min_coef1, min_coef2)));
+      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);
 
@@ -524,155 +698,244 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     });
   }
 
-  CellValue<double> 
-  _scalar_limiter_coefficient(const MeshType& mesh,
-                              const DiscreteFunctionP0<const double>& f,
-                              const DiscreteFunctionDPk<Dimension, const double>& DPk_fh) const
+  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);
-
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+    
     const auto xj = mesh_data.xj();
 
-    CellValue<double> lambda{mesh.connectivity()};
-
-    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-      const double fj = f[cell_id];
+    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]);
+    });
 
-      double f_min = fj;
-      double f_max = fj;
+    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]);
+    }
 
-      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]]);
-      }
+    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];
+            }
+          }
+        }
+        }
+    });
 
-      double f_bar_min = fj;
-      double f_bar_max = fj;
+    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);
+        }
 
-      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]);
+        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;
 
-        f_bar_min = std::min(f_bar_min, f_xk);
-        f_bar_max = std::max(f_bar_max, f_xk);
-      }
+        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 eps = 1E-14;
-      double coef1     = 1;
-      if (std::abs(f_bar_max - fj) > eps) {
-        coef1 = (f_max - fj) / ((f_bar_max - fj));
-      }
+          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);
 
-      double coef2 = 1.;
-      if (std::abs(f_bar_min - fj) > eps) {
-        coef2 = (f_min - fj) / ((f_bar_min - fj));
-      }
+          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);
 
-      lambda[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
-    });
+          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);
+        }
 
-    return lambda;
-  }
+        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);
+        }
 
-  NodeValuePerCell<Rdxd>
-  _compute_tensor_reconstruction(const std::shared_ptr<const MeshVariant>& v_mesh,
-                                 DiscreteTensorFunction sigma,
-                                 PolynomialReconstructionDescriptor reconstruction_descriptor) const
-    {
-        const MeshType& mesh                = *v_mesh->get<MeshType>();
-
-        CellValue<double> limiter_j{mesh.connectivity()};
-        parallel_for(
-          mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
-            limiter_j[j] = 1;
-          }
-        );
-
-        const CellValue<const Rdxd>& sigma_j = copy(sigma.cellValues());
-        std::vector<DiscreteFunctionDPk<Dimension, double>> sigma_coef;
-        std::vector<size_t> row;    
-        std::vector<size_t> col;
-        for(size_t i = 0; i < Dimension; ++i){
-            for(size_t k = i; k < Dimension; ++k){
-
-                CellValue<double> coef{mesh.connectivity()};
-                for(CellId j = 0; j <mesh.numberOfCells(); ++j){
-                    coef[j] = sigma_j[j](i,k);
-                }
-
-                const auto& coef_scalar_function = DiscreteScalarFunction(v_mesh, coef);
-                auto coef_v = std::make_shared<DiscreteFunctionVariant>(coef_scalar_function);
-                auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(coef_v);
-                auto DPk_coef = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
-
-                const CellValue<const double>& limiter_ik = _scalar_limiter_coefficient(mesh,coef_scalar_function,DPk_coef);
-                parallel_for(
-                  mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-                    limiter_j[j] = std::min(limiter_j[j],limiter_ik[j]);
-                });
-
-                sigma_coef.push_back(copy(DPk_coef));
-                row.push_back(i);
-                col.push_back(k);
-            }
+        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);
         }
 
-        NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
-        const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+        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];
+        }
+    });
+  }
 
-        parallel_for(
-          mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
-            const auto& cell_nodes = cell_to_node_matrix[j];
-            for(size_t i = 0; i < Dimension; ++i){
-              for(size_t k = 0; k < Dimension; ++k){
-                for(size_t R = 0; R < cell_nodes.size(); ++R){
-                  sigma_lim(j,R)(i,k) = 0;
-                }
-              }
-            }
-          }
-        );
+    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 NodeValue<const Rd>& xr = copy(mesh.xr());
-        parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-                const auto& cell_nodes = cell_to_node_matrix[j];
+    const auto xj = mesh_data.xj();
 
-                for(size_t R = 0; R < cell_nodes.size(); ++R){
-                    const NodeId r = cell_nodes[R];
+    //A retirer + tard
+    const auto xr = mesh.xr();
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    //
 
-                    for(size_t l = 0; l < sigma_coef.size(); ++l){
-                        const size_t& i = row[l];
-                        const size_t& k = col[l];
+    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;
 
-                        auto coefficients = sigma_coef[l].coefficients(j);
+        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]]);
+        }
 
-                        coefficients[0] = (1-limiter_j[j])*sigma_j[j](i,k) + limiter_j[j] * coefficients[0];
-                        for (size_t indice = 1; indice < coefficients.size(); ++indice) {
-                          coefficients[indice] *= limiter_j[j];
-                        }
+        double inv2_bar_min = inv2j;
+        double inv2_bar_max = inv2j;
 
-                        sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
-                        if(i != k){
-                          sigma_lim(j,R)(i,k) *= 2;
-                        }
-                    }
+        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]]));
 
-                    sigma_lim(j,R) += transpose(sigma_lim(j,R));
-                    sigma_lim(j,R) *= 0.5;
-                }
+          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);
+        }
 
-        return sigma_lim;
-    }
+        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>>
+  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, 
@@ -693,10 +956,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
     return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
             std::make_shared<DiscreteFunctionVariant>(aL_d),
-            std::make_shared<DiscreteFunctionVariant>(aT_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>>
+  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, 
@@ -716,10 +983,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
     return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
             std::make_shared<DiscreteFunctionVariant>(aL_d),
-            std::make_shared<DiscreteFunctionVariant>(aT_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>>
+  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,
@@ -739,70 +1010,187 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     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,lambda,gamma,p_inf);
+      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>>
+  std::tuple<const std::shared_ptr<const ItemValueVariant>,
+             const std::shared_ptr<const SubItemValuePerItemVariant>,
+             const std::shared_ptr<const ItemValueVariant>,
+             const std::shared_ptr<const SubItemValuePerItemVariant>,
+             NodeValue<Rd>,
+             NodeValuePerCell<Rd>>
   compute_fluxes(const SolverType& solver_type,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
-                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+                 const std::vector<NodeId>& map1,
+                 const std::vector<NodeId>& map2) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
-    if (not i_mesh) {
+    std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
+    std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
+    if (not i_mesh1) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
-    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+    if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
       throw NormalError("hyperelastic solver expects P0 functions");
     }
+    if (not i_mesh2) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
 
-    const MeshType& mesh                = *i_mesh->get<MeshType>();
-    const DiscreteScalarFunction& rho   = rho_v->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u     = u_v->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
-    const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+    if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
 
-    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+    const MeshType& mesh1                = *i_mesh1->get<MeshType>();
+    const DiscreteScalarFunction& rho1   = rho1_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u1     = u1_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& aL1    = aL1_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& aT1    = aT1_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
+    const DiscreteScalarFunction& p1     = p1_v->get<DiscreteScalarFunction>();
 
-    for(auto&& bc_descriptor : bc_descriptor_list){
+    const MeshType& mesh2                = *i_mesh2->get<MeshType>();
+    const DiscreteScalarFunction& rho2   = rho2_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u2     = u2_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& aL2    = aL2_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& aT2    = aT2_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
+    const DiscreteScalarFunction& p2     = p2_v->get<DiscreteScalarFunction>();
+
+    const TinyMatrix<Dimension> I    = identity;
+    const DiscreteTensorFunction& S1 = sigma1 + p1*I;
+    const DiscreteTensorFunction& S2 = sigma2 + p2*I;
+
+    const std::shared_ptr<const DiscreteFunctionVariant>& S1_v = std::make_shared<DiscreteFunctionVariant>(S1);
+    const std::shared_ptr<const DiscreteFunctionVariant>& S2_v = std::make_shared<DiscreteFunctionVariant>(S2);
+
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;
+
+    for(auto&& bc_descriptor : bc_descriptor_list1){
       if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
-        symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+        symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+    for(auto&& bc_descriptor : bc_descriptor_list2){
+      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+        symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
       }
     }
 
-    PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
-                                                                symmetry_boundary_descriptor_list);
+    PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list1);
+    PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list2);
 
-    NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor);
-    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
-    auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v, S1_v, p1_v);
+    auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto DPk_Sh1 = reconstruction1[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+    auto DPk_ph1 = reconstruction1[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
 
-    DiscreteFunctionDPk<Dimension, Rd> u_lim  = copy(DPk_uh);
-    this->_vector_limiter(mesh,u,u_lim);
+    DiscreteFunctionDPk<Dimension, Rd> u1_lim       = copy(DPk_uh1);
+    DiscreteFunctionDPk<Dimension, Rdxd> S1_lim     = copy(DPk_Sh1);
+    DiscreteFunctionDPk<Dimension, double> p1_lim   = copy(DPk_ph1);
 
-    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
+    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);
 
-    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
-    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
+    auto copy_DPk_Sh1 = copy(DPk_Sh1);
+    const NodeValuePerCell<double> coef1 = _compute_inv2_coefficient(mesh1, S1, copy_DPk_Sh1);
 
-    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
-    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+    auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v, S2_v, p2_v);
+    auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto DPk_Sh2 = reconstruction2[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+    auto DPk_ph2 = reconstruction2[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
 
-    synchronize(Ar);
-    synchronize(br);
+    DiscreteFunctionDPk<Dimension, Rd> u2_lim       = copy(DPk_uh2);
+    DiscreteFunctionDPk<Dimension, Rdxd> S2_lim     = copy(DPk_Sh2);
+    DiscreteFunctionDPk<Dimension, double> p2_lim   = copy(DPk_ph2);
 
-    NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
+    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);
 
-    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
-                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+    auto copy_DPk_Sh2 = copy(DPk_Sh2);
+    const NodeValuePerCell<double> coef2 = _compute_inv2_coefficient(mesh2, S2, copy_DPk_Sh2);
+
+    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<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, S2_lim, p2_lim, lambda2, coef2);
+
+    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+    this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+    this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+
+    synchronize(Ar1);
+    synchronize(br1);
+    synchronize(Ar2);
+    synchronize(br2);
+
+    NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
+    NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);
+
+    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);
+
+    NodeValue<Rd> CR_ur         = copy(ur2);
+    NodeValuePerCell<Rd> CR_Fjr = copy(Fjr2);
+
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+                           std::make_shared<const ItemValueVariant>(ur2),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr2), CR_ur, CR_Fjr);
   }
 
   std::tuple<const std::shared_ptr<const ItemValueVariant>,
@@ -872,32 +1260,13 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
       }
     }
 
-    PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
-                                                                 symmetry_boundary_descriptor_list1);
-    PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
-                                                                 symmetry_boundary_descriptor_list2);
-
-    NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1);
-    NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2);
-
-    auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v);
-    auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
-    auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v);
-    auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
-
-    DiscreteFunctionDPk<Dimension, Rd> u1_lim  = copy(DPk_uh1);
-    this->_vector_limiter(mesh1,u1,u1_lim);
-    DiscreteFunctionDPk<Dimension, Rd> u2_lim  = copy(DPk_uh2);
-    this->_vector_limiter(mesh2,u2,u2_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, sigma1_lim);
+    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1, sigma1);
     NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
-    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim);
-    this->_applyCouplingBC(Ar1, Ar2, br1, br2, map1, map2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2, sigma2);
 
     const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
     const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
@@ -910,11 +1279,15 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     synchronize(br2);
 
     NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
-    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim);
     NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);
-    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
-    NodeValue<Rd> CR_ur         = this->_computeUr(mesh2, Ar2, br2);
-    NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
+
+    this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+
+    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
+    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+
+    NodeValue<Rd> CR_ur         = copy(ur2);
+    NodeValuePerCell<Rd> CR_Fjr = copy(Fjr2);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
@@ -929,6 +1302,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                  const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
                  const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p_v, 
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                  NodeValue<Rd> CR_ur,
                  NodeValuePerCell<Rd> CR_Fjr,
@@ -949,6 +1323,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
     const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
     const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+    const DiscreteScalarFunction& p     = p_v->get<DiscreteScalarFunction>();
+
+    const TinyMatrix<Dimension> I = identity;
+    const DiscreteTensorFunction& S     = sigma + p*I;
+    const std::shared_ptr<const DiscreteFunctionVariant>& S_v = std::make_shared<DiscreteFunctionVariant>(S);
 
     std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
 
@@ -961,148 +1340,89 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
                                                                 symmetry_boundary_descriptor_list);
 
-    NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor);
-    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
-    auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v, S_v, p_v);
+    auto DPk_uh     = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto DPk_Sh     = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+    auto DPk_ph     = reconstruction[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
+
+    DiscreteFunctionDPk<Dimension, Rd> u_lim       = copy(DPk_uh);
+    DiscreteFunctionDPk<Dimension, Rdxd> S_lim     = copy(DPk_Sh);
+    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);
 
-    DiscreteFunctionDPk<Dimension, Rd> u_lim  = copy(DPk_uh);
-    this->_vector_limiter(mesh,u,u_lim);
+    auto copy_DPk_Sh = copy(DPk_Sh);
+    const NodeValuePerCell<double> coef = _compute_inv2_coefficient(mesh, S, copy_DPk_Sh);
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
-    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, S_lim, p_lim, lambda, coef);
 
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
 
-    synchronize(Ar);
-    synchronize(br);
-
-    NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
-
-    this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
-
-    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
-                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
-  }
-
-    std::tuple<const std::shared_ptr<const ItemValueVariant>,
-             const std::shared_ptr<const SubItemValuePerItemVariant>,
-             const std::shared_ptr<const ItemValueVariant>,
-             const std::shared_ptr<const SubItemValuePerItemVariant>,
-             NodeValue<Rd>,
-             NodeValuePerCell<Rd>>
-  compute_fluxes2(const SolverType& solver_type,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
-                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
-                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
-                 const std::vector<NodeId>& map1,
-                 const std::vector<NodeId>& map2) const
-  {
-    std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
-    std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
-    if (not i_mesh1) {
-      throw NormalError("discrete functions are not defined on the same mesh ici2");
-    }
-
-    if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
-      throw NormalError("hyperelastic solver expects P0 functions");
-    }
-    if (not i_mesh2) {
-      throw NormalError("discrete functions are not defined on the same mesh ici2 mesh2");
-    }
-
-    if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
-      throw NormalError("acoustic solver expects P0 functions");
-    }
+    synchronize(Ar);
+    synchronize(br);
 
-    const MeshType& mesh1                = *i_mesh1->get<MeshType>();
-    const DiscreteScalarFunction& rho1   = rho1_v->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u1     = u1_v->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& aL1    = aL1_v->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& aT1    = aT1_v->get<DiscreteScalarFunction>();
-    const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
+    NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
+    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim, lambda, coef);
 
-    const MeshType& mesh2                = *i_mesh2->get<MeshType>();
-    const DiscreteScalarFunction& rho2   = rho2_v->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u2     = u2_v->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& aL2    = aL2_v->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& aT2    = aT2_v->get<DiscreteScalarFunction>();
-    const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
+    this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
 
-    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
-    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+  }
 
-    for(auto&& bc_descriptor : bc_descriptor_list1){
-      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
-        symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
-      }
-    }
-    for(auto&& bc_descriptor : bc_descriptor_list2){
-      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
-        symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
-      }
+  std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+  compute_fluxes(const SolverType& solver_type,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                 NodeValue<Rd> CR_ur,
+                 NodeValuePerCell<Rd> CR_Fjr,
+                 const std::vector<NodeId> map2) const
+  {
+    std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+    if (not i_mesh) {
+      throw NormalError("discrete functions are not defined on the same mesh");
     }
 
-    PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
-                                                                 symmetry_boundary_descriptor_list1);
-    PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
-                                                                 symmetry_boundary_descriptor_list2);
-
-    NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1);
-    NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2);
+    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperelastic solver expects P0 functions");
+    }
 
-    auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v);
-    auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
-    auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v);
-    auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    const MeshType& mesh                = *i_mesh->get<MeshType>();
+    const DiscreteScalarFunction& rho   = rho_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u     = u_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
 
-    DiscreteFunctionDPk<Dimension, Rd> u1_lim  = copy(DPk_uh1);
-    this->_vector_limiter(mesh1,u1,u1_lim);
-    DiscreteFunctionDPk<Dimension, Rd> u2_lim  = copy(DPk_uh2);
-    this->_vector_limiter(mesh2,u2,u2_lim);
+    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
-    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> Ar = this->_computeAr(mesh, Ajr);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, sigma);
 
-    NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
-    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1_lim, sigma1_lim);
-    NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
-    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim);
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
 
-    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
-    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
-    this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
-    this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+    synchronize(Ar);
+    synchronize(br);
 
-    synchronize(Ar1);
-    synchronize(br1);
-    synchronize(Ar2);
-    synchronize(br2);
+    NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
+    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
 
-    NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
-    NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);
-    this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
-    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim);
-    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
+    this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
 
-    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
-                           std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
-                           std::make_shared<const ItemValueVariant>(ur2),
-                           std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
-                           ur2,
-                           Fjr2);
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
   }
 
   std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1116,21 +1436,24 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                const DiscreteVectorFunction& u,
                const DiscreteScalarFunction& E,
                const DiscreteTensorFunction& CG,
-               const NodeValue<const Rd>& ur,
-               const NodeValuePerCell<const Rd>& Fjr) const
+               const NodeValue<const Rd>& ur_o2,
+               const NodeValue<const Rd>& ur_o1,
+               const NodeValuePerCell<const Rd>& Fjr_o2,
+               const NodeValuePerCell<const Rd>& Fjr_o1) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
+    NodeValue<Rd> ur         = copy(ur_o2);
+    NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);
+
     if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
         (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
       throw NormalError("fluxes are not defined on the same connectivity than the mesh");
     }
 
-    NodeValue<Rd> new_xr = copy(mesh.xr());
-    parallel_for(
-      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
-
-    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+    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);
 
     CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
     const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
@@ -1140,29 +1463,105 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     CellValue<double> new_E   = copy(E.cellValues());
     CellValue<Rdxd> new_CG    = copy(CG.cellValues());
 
+    CellValue<double> new_rho_stencil = copy(rho.cellValues());
+    CellValue<Rd> new_u_stencil       = copy(u.cellValues());
+    CellValue<double> new_E_stencil   = copy(E.cellValues());
+    CellValue<Rdxd> new_CG_stencil    = copy(CG.cellValues());
+
+    CellValue<int> cell_flag{mesh.connectivity()};
+
     parallel_for(
       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-        const auto& cell_nodes = cell_to_node_matrix[j];
+        auto cell_nodes = cell_to_node_matrix[j];
 
         Rd momentum_fluxes   = zero;
         double energy_fluxes = 0;
         Rdxd gradv           = zero;
         for (size_t R = 0; R < cell_nodes.size(); ++R) {
-          const NodeId r = cell_nodes[R];
+          NodeId r = cell_nodes[R];
           gradv += tensorProduct(ur[r], Cjr(j, R));
           momentum_fluxes += Fjr(j, R);
           energy_fluxes += dot(Fjr(j, R), ur[r]);
         }
-        const Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
-        const double dt_over_Mj        = dt / (rho[j] * Vj[j]);
-        const double dt_over_Vj        = dt / Vj[j];
+        Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
+        double dt_over_Mj        = dt / (rho[j] * Vj[j]);
+        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;
+
+        double new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
+        if(new_eps < 0.){
+          cell_flag[j] = 1;
+        }else{
+          cell_flag[j] = 0;
+        }
       });
 
+    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
+      if(cell_flag[j] == 1){
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        Rdxd gradv           = zero;
+        std::vector<NodeId> NodeList;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          NodeId r = cell_nodes[R];
+          NodeList.push_back(r);
+          Fjr(j,R) = Fjr_o1(j,R);
+          ur[r] = ur_o1[r];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += dot(Fjr(j, R), ur[r]);
+        }
+        Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
+        double dt_over_Mj        = dt / (rho[j] * Vj[j]);
+        double dt_over_Vj        = dt / Vj[j];
+        new_u[j] = new_u_stencil[j] +  dt_over_Mj * momentum_fluxes;
+        new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
+        new_CG[j] = new_CG_stencil[j] + dt_over_Vj * cauchy_green_fluxes;
+        new_CG[j] += transpose(new_CG[j]);
+        new_CG[j] *= 0.5;
+
+        auto cell_stencil = stencil[j];
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          CellId J = cell_stencil[i_cell];
+          auto i_cell_nodes = cell_to_node_matrix[J];
+          momentum_fluxes = zero;
+          energy_fluxes   = 0;
+          gradv           = zero;
+          for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
+            NodeId i_node = i_cell_nodes[R];
+            for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
+              if(NodeList[node_test] == i_node){
+                Fjr(J,R) = Fjr_o1(J,R);
+              }
+            }
+            gradv += tensorProduct(ur[i_node], Cjr(J, R));
+            momentum_fluxes += Fjr(J, R);
+            energy_fluxes += dot(Fjr(J, R), ur[i_node]);
+          }
+          cauchy_green_fluxes = gradv * CG[J] + CG[J] * transpose(gradv);
+          dt_over_Mj          = dt / (rho[J] * Vj[J]);
+          dt_over_Vj          = dt / Vj[J];
+          new_u[J] = new_u_stencil[J] +  dt_over_Mj * momentum_fluxes;
+          new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
+          new_CG[J] = new_CG_stencil[J] + dt_over_Vj * cauchy_green_fluxes;
+          new_CG[J] += transpose(new_CG[J]);
+          new_CG[J] *= 0.5;
+        } 
+      }
+    }
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
     CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
 
     parallel_for(
@@ -1186,11 +1585,13 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                const std::shared_ptr<const DiscreteFunctionVariant>& E,
                const std::shared_ptr<const DiscreteFunctionVariant>& CG,
                const std::shared_ptr<const ItemValueVariant>& ur,
-               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
+               const std::shared_ptr<const ItemValueVariant>& ur_o1,
+               const std::shared_ptr<const SubItemValuePerItemVariant> Fjr,
+               const std::shared_ptr<const SubItemValuePerItemVariant> Fjr_o1) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
-    if (not i_mesh) {
-      throw NormalError("discrete functions are not defined on the same mesh ici3");
+    std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
+    if (not mesh_v) {
+      throw NormalError("discrete functions are not defined on the same mesh");
     }
 
     if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
@@ -1198,13 +1599,16 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     }
 
     return this->apply_fluxes(dt,                                   //
-                              *i_mesh->get<MeshType>(),             //
+                              *mesh_v->get<MeshType>(),             //
                               rho->get<DiscreteScalarFunction>(),   //
                               u->get<DiscreteVectorFunction>(),     //
                               E->get<DiscreteScalarFunction>(),     //
                               CG->get<DiscreteTensorFunction>(),    //
-                              ur->get<NodeValue<const Rd>>(),       //
-                              Fjr->get<NodeValuePerCell<const Rd>>());
+                              ur->get<NodeValue<const Rd>>(),
+                              ur_o1->get<NodeValue<const Rd>>(),       //
+                              Fjr->get<NodeValuePerCell<const Rd>>(),
+                              Fjr_o1->get<NodeValuePerCell<const Rd>>()
+                              );
   }
 
   std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1213,28 +1617,31 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>>
   order2_apply_fluxes(const double& dt,
-               const MeshType& mesh,
-               const DiscreteScalarFunction& rho,
-               const DiscreteVectorFunction& u,
-               const DiscreteScalarFunction& E,
-               const DiscreteTensorFunction& CG,
-               const DiscreteScalarFunction& rho_app,
-               const DiscreteTensorFunction& CG_app,
-               const NodeValue<const Rd>& ur,
-               const NodeValuePerCell<const Rd>& Fjr) const
+                      const MeshType& mesh,
+                      const DiscreteScalarFunction& rho,
+                      const DiscreteVectorFunction& u,
+                      const DiscreteScalarFunction& E,
+                      const DiscreteTensorFunction& CG,
+                      const DiscreteScalarFunction& rho_app,
+                      const DiscreteTensorFunction& CG_app,
+                      const NodeValue<const Rd>& ur_o2,
+                      const NodeValue<const Rd>& ur_o1,
+                      const NodeValuePerCell<const Rd>& Fjr_o2,
+                      const NodeValuePerCell<const Rd>& Fjr_o1) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
+    NodeValue<Rd> ur         = copy(ur_o2);
+    NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);
+
     if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
         (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
       throw NormalError("fluxes are not defined on the same connectivity than the mesh");
     }
 
-    NodeValue<Rd> new_xr = copy(mesh.xr());
-    parallel_for(
-      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
-
-    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+    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);
 
     CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
     const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
@@ -1244,6 +1651,13 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     CellValue<double> new_E   = copy(E.cellValues());
     CellValue<Rdxd> new_CG    = copy(CG.cellValues());
 
+    CellValue<double> new_rho_stencil = copy(rho.cellValues());
+    CellValue<Rd> new_u_stencil       = copy(u.cellValues());
+    CellValue<double> new_E_stencil   = copy(E.cellValues());
+    CellValue<Rdxd> new_CG_stencil    = copy(CG.cellValues());
+
+    CellValue<int> cell_flag{mesh.connectivity()};
+
     parallel_for(
       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
@@ -1264,8 +1678,75 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
         new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
         new_CG[j] += transpose(new_CG[j]);
         new_CG[j] *= 0.5;
+
+        const double& new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
+        if(new_eps < 0.){
+          cell_flag[j] = 1;
+        }else{
+          cell_flag[j] = 0;
+        }
       });
 
+    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
+      if(cell_flag[j] == 1){
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        Rdxd gradv           = zero;
+        std::vector<NodeId> NodeList;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          NodeId r = cell_nodes[R];
+          NodeList.push_back(r);
+          Fjr(j,R) = Fjr_o1(j,R);
+          ur[r] = ur_o1[r];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += dot(Fjr(j, R), ur[r]);
+        }
+        Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
+        double dt_over_Mj        = dt / (rho[j] * Vj[j]);
+        new_u[j] = new_u_stencil[j] + dt_over_Mj * momentum_fluxes;
+        new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
+        new_CG[j] = new_CG_stencil[j] + dt_over_Mj * cauchy_green_fluxes;
+        new_CG[j] += transpose(new_CG[j]);
+        new_CG[j] *= 0.5;
+
+        auto cell_stencil = stencil[j];
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          CellId J = cell_stencil[i_cell];
+          auto i_cell_nodes = cell_to_node_matrix[J];
+          momentum_fluxes = zero;
+          energy_fluxes   = 0;
+          gradv           = zero;
+          for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
+            NodeId i_node = i_cell_nodes[R];
+            for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
+              if(NodeList[node_test] == i_node){
+                Fjr(J,R) = Fjr_o1(J,R);
+              }
+            }
+            gradv += tensorProduct(ur[i_node], Cjr(J, R));
+            momentum_fluxes += Fjr(J, R);
+            energy_fluxes += dot(Fjr(J, R), ur[i_node]);
+          }
+          cauchy_green_fluxes = rho_app[J] * (gradv * CG_app[J] + CG_app[J] * transpose(gradv));
+          dt_over_Mj        = dt / (rho[J] * Vj[J]);
+          new_u[J] = new_u_stencil[J] + dt_over_Mj * momentum_fluxes;
+          new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
+          new_CG[J] = new_CG_stencil[J] + dt_over_Mj * cauchy_green_fluxes;
+          new_CG[J] += transpose(new_CG[J]);
+          new_CG[J] *= 0.5;
+        } 
+      }
+    }
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
     CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
 
     parallel_for(
@@ -1284,34 +1765,38 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>>
   order2_apply_fluxes(const double& dt,
-               const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-               const std::shared_ptr<const DiscreteFunctionVariant>& u,
-               const std::shared_ptr<const DiscreteFunctionVariant>& E,
-               const std::shared_ptr<const DiscreteFunctionVariant>& CG,
-               const std::shared_ptr<const DiscreteFunctionVariant>& rho_app,
-               const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
-               const std::shared_ptr<const ItemValueVariant>& ur,
-               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
+                      const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                      const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                      const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                      const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+                      const std::shared_ptr<const DiscreteFunctionVariant>& rho_app,
+                      const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
+                      const std::shared_ptr<const ItemValueVariant>& ur,
+                      const std::shared_ptr<const ItemValueVariant>& ur_o1,
+                      const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+                      const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_o1) const
   {
     std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
     if (not mesh_v) {
-      throw NormalError("discrete functions are not defined on the same mesh ici4");
+      throw NormalError("discrete functions are not defined on the same mesh");
     }
 
     if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
       throw NormalError("hyperelastic solver expects P0 functions");
     }
 
-    return this->order2_apply_fluxes(dt,                                   //
-                              *mesh_v->get<MeshType>(),             //
-                              rho->get<DiscreteScalarFunction>(),   //
-                              u->get<DiscreteVectorFunction>(),     //
-                              E->get<DiscreteScalarFunction>(),     //
-                              CG->get<DiscreteTensorFunction>(),    //
-                              rho_app->get<DiscreteScalarFunction>(),   //
-                              CG_app->get<DiscreteTensorFunction>(),    //
-                              ur->get<NodeValue<const Rd>>(),       //
-                              Fjr->get<NodeValuePerCell<const Rd>>());
+    return this->order2_apply_fluxes(dt,                                       //
+                                     *mesh_v->get<MeshType>(),                 //
+                                     rho->get<DiscreteScalarFunction>(),       //
+                                     u->get<DiscreteVectorFunction>(),         //
+                                     E->get<DiscreteScalarFunction>(),         //
+                                     CG->get<DiscreteTensorFunction>(),        //
+                                     rho_app->get<DiscreteScalarFunction>(),   //
+                                     CG_app->get<DiscreteTensorFunction>(),    //
+                                     ur->get<NodeValue<const Rd>>(),           //
+                                     ur_o1->get<NodeValue<const Rd>>(),
+                                     Fjr->get<NodeValuePerCell<const Rd>>(),
+                                     Fjr_o1->get<NodeValuePerCell<const Rd>>());
   }
 
   std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1373,173 +1858,155 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
   }
 
   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_classic_midpoint_method(const SolverType& solver_type,
-                                   const size_t& law,
-                                   const double& dt1,
-                                   const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                   const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                   const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                   const std::shared_ptr<const DiscreteFunctionVariant>& CG,
-                                   const std::shared_ptr<const DiscreteFunctionVariant>& aL,
-                                   const std::shared_ptr<const DiscreteFunctionVariant>& aT,
-                                   const std::shared_ptr<const DiscreteFunctionVariant>& sigma,                               
-                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                                   NodeValue<Rd> CR_ur,
-                                   NodeValuePerCell<Rd> CR_Fjr,
-                                   std::vector<NodeId>& map2,
-                                   size_t chi_solid,
-                                   const double& lambda,
-                                   const double& mu,
-                                   const double& gamma,
-                                   const double& p_inf) const
-  {
-    
-    std::shared_ptr m_app = getCommonMesh({rho, aL, aT, u, sigma});
-    std::shared_ptr<const DiscreteFunctionVariant> rho_app = rho;
-    std::shared_ptr<const DiscreteFunctionVariant> u_app   = u;
-    std::shared_ptr<const DiscreteFunctionVariant> E_app   = E;
-    std::shared_ptr<const DiscreteFunctionVariant> CG_app  = CG;
-
-    std::shared_ptr new_m = getCommonMesh({rho, aL, aT, u, sigma});
-    std::shared_ptr<const DiscreteFunctionVariant> new_rho = rho;
-    std::shared_ptr<const DiscreteFunctionVariant> new_u   = u;
-    std::shared_ptr<const DiscreteFunctionVariant> new_E   = E;
-    std::shared_ptr<const DiscreteFunctionVariant> new_CG  = CG;
-    std::shared_ptr<const DiscreteFunctionVariant> new_sigma = sigma;
-    std::shared_ptr<const DiscreteFunctionVariant> new_aL   = aL;
-    std::shared_ptr<const DiscreteFunctionVariant> new_aT   = aT;
-
-    double sum_dt = 0;
-
-    while(sum_dt < dt1){
-      
-      DiscreteScalarFunction new_aL_d = new_aL->get<DiscreteScalarFunction>();
-      DiscreteScalarFunction new_aT_d = new_aT->get<DiscreteScalarFunction>();
-      const std::shared_ptr<const DiscreteFunctionVariant>& new_c = std::make_shared<const DiscreteFunctionVariant>(new_aL_d + new_aT_d);
-
-      double dt2 = 0.4 * hyperelastic_dt(new_c);
-      if(sum_dt + dt2 > dt1){
-        dt2 = dt1 - sum_dt;
-      }
-
-      auto [ur, Fjr] = compute_fluxes(solver_type,new_rho,new_aL,new_aT,new_u,new_sigma,bc_descriptor_list,CR_ur,CR_Fjr,map2);
-      std::tie(m_app, rho_app, u_app, E_app, CG_app) = apply_fluxes(dt2/2,new_rho,new_u,new_E,new_CG,ur,Fjr);
-
-      const auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid,lambda,mu,gamma,p_inf);
-
-      auto [ur_app, Fjr_app] = compute_fluxes(solver_type,rho_app,aL_app,aT_app,u_app,sigma_app,bc_descriptor_list,CR_ur,CR_Fjr,map2);
-      std::tie(new_m,new_rho,new_u,new_E,new_CG) = order2_apply_fluxes(dt2,new_rho,new_u,new_E,new_CG,rho_app,CG_app,ur_app,Fjr_app); 
-
-      std::tie(new_sigma,new_aL,new_aT) = _apply_state_law(law,new_rho,new_u,new_E,new_CG,chi_solid,lambda,mu,gamma,p_inf);
-        
-      sum_dt += dt2;
-
-    }
-
-    return {new_m, new_rho, new_u, new_E, 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>,
-             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
+  compute_sub_iterations(const SolverType& solver_type,
+                         const size_t& law,
+                         const double& CFL,
+                         const double& dt,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& 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 size_t& chi_solid,
+                         const double& lambda,
+                         const double& mu,
+                         const double& gamma,
+                         const double& p_inf) 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);
+    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;
 
-    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;
+    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);
 
-    size_t fluid = 0;
-    if ((mu == 0) and (lambda == 0)) {
-      fluid = 1;
-    }
+      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);
 
-    for (size_t i = 0; i < q; i++) {
-      if (i == q - 1) {
-        dt2 = dt1 - sum_dt;
+      double sub_dt = CFL * hyperelastic_dt(c);
+      if(sum_dt + sub_dt > dt){
+	      sub_dt = dt - 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);
+      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);
 
-      const DiscreteScalarFunction& aT_d = sqrt(mu / rho_d);
+      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);
 
-      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;
+      sum_dt += sub_dt;
     }
 
-    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};
+    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>,
@@ -1567,6 +2034,8 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
         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 double& lambda_f,
@@ -1581,79 +2050,51 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
     std::shared_ptr m1     = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
 
-    std::shared_ptr m1_app = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
-    std::shared_ptr<const DiscreteFunctionVariant> rho1_app = rho1;
-    std::shared_ptr<const DiscreteFunctionVariant> u1_app   = u1;
-    std::shared_ptr<const DiscreteFunctionVariant> E1_app   = E1;
-    std::shared_ptr<const DiscreteFunctionVariant> CG1_app  = CG1;
-
-    std::shared_ptr m2_app = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
-    std::shared_ptr<const DiscreteFunctionVariant> rho2_app   = rho2;
-    std::shared_ptr<const DiscreteFunctionVariant> u2_app     = u2;
-    std::shared_ptr<const DiscreteFunctionVariant> E2_app     = E2;
-    std::shared_ptr<const DiscreteFunctionVariant> CG2_app    = CG2;
+    //std::shared_ptr m1_app = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
+    //std::shared_ptr<const DiscreteFunctionVariant> rho1_app = rho1;
+    //std::shared_ptr<const DiscreteFunctionVariant> u1_app   = u1;
+    //std::shared_ptr<const DiscreteFunctionVariant> E1_app   = E1;
+    //std::shared_ptr<const DiscreteFunctionVariant> CG1_app  = CG1;
+
+    //std::shared_ptr m2_app = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+    //std::shared_ptr<const DiscreteFunctionVariant> rho2_app   = rho2;
+    //std::shared_ptr<const DiscreteFunctionVariant> u2_app     = u2;
+    //std::shared_ptr<const DiscreteFunctionVariant> E2_app     = E2;
+    //std::shared_ptr<const DiscreteFunctionVariant> CG2_app    = CG2;
     std::shared_ptr<const DiscreteFunctionVariant> aL2_app    = aL2;
     std::shared_ptr<const DiscreteFunctionVariant> aT2_app    = aT2;
     std::shared_ptr<const DiscreteFunctionVariant> sigma2_app = sigma2;
+    std::shared_ptr<const DiscreteFunctionVariant> p2_app     = p2;
 
 
     auto [map1, map2] = _computeMapping(m1,m2,bc_descriptor_list1,bc_descriptor_list2);
 
-    auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] = 
-      compute_fluxes2(solver_type, rho1, aL1, aT1, u1, sigma1, bc_descriptor_list1, rho2, aL2, aT2, u2, sigma2,
-                     bc_descriptor_list2, map1, map2);                                                              //Fluxes t^n
-    std::tie(m1_app,rho1_app,u1_app,E1_app,CG1_app) = apply_fluxes(dt1/2, rho1, u1, E1, CG1, ur1, Fjr1);            //Fluid t^n+(dt*0.5)
-
-    DiscreteScalarFunction aL_d = aL2->get<DiscreteScalarFunction>();
-    DiscreteScalarFunction aT_d = aT2->get<DiscreteScalarFunction>();
-    const std::shared_ptr<const DiscreteFunctionVariant>& c = std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
-
-    double dt2 = 0.4 * hyperelastic_dt(c);
-    if(dt2 > dt1/2){
-	    dt2 = dt1/2;
-    }
-
-    std::tie(m2_app,rho2_app,u2_app,E2_app,CG2_app) = apply_fluxes(dt2, rho2, u2, E2, CG2, ur2, Fjr2);              //Solid t^n + dt2
+    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);
 
-    double sum_dt = dt2;
+    //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);
 
     size_t chi_solid = 0;
-    if((mu_s!=0) and (lambda_s!=0)){
+    if((mu_s!=0) or (lambda_s!=0)){
       chi_solid = 1;
     }
     const size_t chi_fluid = 1-chi_solid;
 
-    while(sum_dt < dt1/2){  //Partie a revoir avec loi de comportement diff
-      std::tie(sigma2_app, aL2_app, aT2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
-
-      DiscreteScalarFunction aL2_d_app = aL2_app->get<DiscreteScalarFunction>();
-      DiscreteScalarFunction aT2_d_app = aT2_app->get<DiscreteScalarFunction>();
-      const std::shared_ptr<const DiscreteFunctionVariant>& c_app =
-        std::make_shared<const DiscreteFunctionVariant>(aL2_d_app + aT2_d_app);
-
-      dt2 = 0.4 * hyperelastic_dt(c_app);
-      if(sum_dt + dt2 > dt1/2){
-	      dt2 = dt1/2 - sum_dt;
-      }
-
-      auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, rho2_app, aL2_app, aT2_app, u2_app, sigma2_app,
-                                                bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+    auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = compute_sub_iterations(solver_type,law,0.04,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
 
-      std::tie(m2_app, rho2_app, u2_app, E2_app, CG2_app) = 
-        apply_fluxes(dt2, rho2_app, u2_app, E2_app, CG2_app, sub_ur2, sub_Fjr2);
-
-      sum_dt += dt2;
-    }
+    //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);
+    //std::tie(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 auto& [sigma1_app, aL1_app, aT1_app] = _apply_state_law(law,rho1_app,u1_app,E1_app,CG1_app,chi_fluid,lambda_f,mu_f,gamma_f,p_inf_f);
-    std::tie(sigma2_app, aL2_app, aT2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+    //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_app_o1, Fjr1_app_o1, ur2_app_o1, Fjr2_app_o1, CR_ur_app_o1, CR_Fjr_app_o1] = 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);
 
-    auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app] = compute_fluxes2(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); //Fluxes t^n + dt*0.5
+    //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_app_o1, Fjr1_app, Fjr1_app_o1);
+    //auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,law,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,p2,bc_descriptor_list2,
+    //                                                                                    CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
 
-    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, Fjr1_app);
-    auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,law,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,bc_descriptor_list2,
-                                                                                        CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+    const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);
 
     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);
 
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.hpp b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
index da09ae613..b0682e12a 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -26,61 +26,37 @@ class Order2LocalDtHyperelasticSolverHandler
  private:
   struct IOrder2LocalDtHyperelasticSolver
   {
-    virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
-                       const std::shared_ptr<const SubItemValuePerItemVariant>>
-    compute_fluxes(
-      const SolverType& solver_type,
-      const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-      const std::shared_ptr<const DiscreteFunctionVariant>& aL,
-      const std::shared_ptr<const DiscreteFunctionVariant>& aT,
-      const std::shared_ptr<const DiscreteFunctionVariant>& u,
-      const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
-      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
-
-    virtual std::tuple<std::shared_ptr<const MeshVariant>,
-                       std::shared_ptr<const DiscreteFunctionVariant>,
-                       std::shared_ptr<const DiscreteFunctionVariant>,
-                       std::shared_ptr<const DiscreteFunctionVariant>,
-                       std::shared_ptr<const DiscreteFunctionVariant>>
-    apply_fluxes(const double& dt,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                 const std::shared_ptr<const DiscreteFunctionVariant>& CG,
-                 const std::shared_ptr<const ItemValueVariant>& ur,
-                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) 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 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 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>,
@@ -109,6 +85,8 @@ class Order2LocalDtHyperelasticSolverHandler
         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 double& lambda_f,
-- 
GitLab