diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 33c1f7b501e7ddbc68e2a61acbc66a2ed54f5ce2..7512782764f0c31a201c8f637cca5401327c5f27 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -999,7 +999,7 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver_order2",
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
                             std::function(
 
                               [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -1020,7 +1020,9 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list2,
-                                 const double& mu, const double& lambda, const double& dt1)
+                                 const double& lambda1, const double& mu1,
+                                 const double& lambda2, const double& mu2, 
+                                 const double& gamma1,  const double& gamma2 , const double& dt1)
                                 -> std::tuple<
                                   std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
                                   std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1035,9 +1037,55 @@ SchemeModule::SchemeModule()
                                                                         getCommonMesh(
                                                                           {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
                                   .solver()
-                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1,
+                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
                                          u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
-                                         bc_descriptor_list2, mu, lambda);
+                                         bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, 0, gamma2, 0);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& lambda1, const double& mu1,
+                                 const double& lambda2, const double& mu2, 
+                                 const double& gamma1,  const double& p_inf1,
+                                 const double& gamma2,  const double& p_inf2, const double& dt1)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, 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 c249ab60f934d2ee241c620d526871e2e1080919..83b802f5185ad481eba9de4cb70553f0d0803873 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -185,7 +185,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
   _computeBr(const Mesh<Dimension>& mesh,
              const NodeValuePerCell<const Rdxd>& Ajr,
              DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
-             NodeValuePerCell<const Rdxd> DPk_sigma) const
+             DiscreteFunctionDPk<Dimension, const Rdxd> DPk_sigma) const
   {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
@@ -206,7 +206,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
         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) * DPk_u[J](xr[r]) - DPk_sigma(J, R) * Cjr(J, R);
+          br += Ajr(J, R) * DPk_u[J](xr[r]) - DPk_sigma[j](xr[r]) * Cjr(J, R);
         }
 
         b[r] = br;
@@ -338,7 +338,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
               const NodeValuePerCell<const Rdxd>& Ajr,
               const NodeValue<const Rd>& ur,
               DiscreteFunctionDPk<Dimension,const Rd> DPk_uh,
-              NodeValuePerCell<const Rdxd> DPk_sigma) const
+              DiscreteFunctionDPk<Dimension,const Rdxd> DPk_sigma) const
   {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
@@ -357,13 +357,97 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
         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];
-          F(J,R) = -Ajr(J,R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma(J,R) * Cjr(J,R);
+          F(J,R) = -Ajr(J,R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma[J](xr[r]) * Cjr(J,R);
         }
       });
 
     return F;
   }
 
+  void 
+  _tensor_limiter(const MeshType& mesh,
+                  const DiscreteFunctionP0<const Rdxd>& f,
+                  DiscreteFunctionDPk<Dimension, Rdxd>& DPk_fh) const
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
+
+    const auto xj = mesh_data.xj();
+
+    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+      const Rdxd fj = f[cell_id];
+
+      Rdxd f_min = fj;
+      Rdxd 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 row = 0; row < Dimension; ++row){
+          for(size_t col = 0; col < Dimension; ++col){
+            f_min(row,col) = std::min(f_min(row,col), f[cell_stencil[i_cell]](row,col));
+            f_max(row,col) = std::max(f_max(row,col), f[cell_stencil[i_cell]](row,col));
+          }
+        }
+      }
+
+      Rdxd f_bar_min = fj;
+      Rdxd 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 Rdxd f_xk        = DPk_fh[cell_id](xj[cell_k_id]);
+
+        for(size_t row = 0; row < Dimension; ++row){
+          for(size_t col = 0; col < Dimension; ++col){
+            f_bar_min(row,col) = std::min(f_bar_min(row,col), f_xk(row,col));
+            f_bar_max(row,col) = std::max(f_bar_max(row,col), f_xk(row,col));
+          }
+        }
+      }
+
+      const double eps = 1E-14;
+      Rdxd coef1;
+      for(size_t row = 0; row < Dimension; ++row){
+        for(size_t col = 0; col < Dimension; ++col){
+          coef1(row,col) = 1;
+          if(std::abs(f_bar_min(row,col) - fj(row,col)) > eps){
+            coef1(row,col) = (f_max(row,col) - fj(row,col)) / (f_bar_max(row,col) - fj(row,col));
+          }
+        }
+      }
+
+      Rdxd coef2;
+      for(size_t row = 0; row < Dimension; ++row){
+        for(size_t col = 0; col < Dimension; ++col){
+          coef2(row,col) = 1;
+          if (std::abs(f_bar_min(row,col) - fj(row,col)) > eps) {
+            coef2(row,col) = (f_min(row,col) - fj(row,col)) / ((f_bar_min(row,col) - fj(row,col)));
+          }
+        }
+      }
+
+      double min_coef1 = coef1(0,0);
+      double min_coef2 = coef2(0,0);
+      for(size_t row = 0; row < Dimension; ++row){
+        for(size_t col = 0; col < Dimension; ++col){
+          min_coef1 = std::min(min_coef1,coef1(row,col));
+          min_coef2 = std::min(min_coef2,coef2(row,col));
+        }
+      }
+
+      const double lambda = std::max(0., std::min(1.,std::min(min_coef1, min_coef2))); 
+
+      auto coefficients = DPk_fh.coefficients(cell_id);
+
+      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+      for (size_t i = 1; i < coefficients.size(); ++i) {
+        coefficients[i] *= lambda;
+      } 
+    });
+  }
+
    void
   _vector_limiter(const MeshType& mesh,
                   const DiscreteFunctionP0<const Rd>& f,
@@ -530,7 +614,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
     const TinyMatrix<Dimension> I = identity;
 
     const DiscreteScalarFunction& eps         = E - 0.5 * dot(u,u);
-    const DiscreteScalarFunction& p_d         = (1-chi_solid)*(gamma_f-1) * rho * eps - p_inf_f*gamma_f + chi_solid * (gamma_s - 1) * rho * eps - gamma_s * p_inf_s;
+    const DiscreteScalarFunction& p_d         = (1-chi_solid)*((gamma_f-1) * rho * eps)+ chi_solid * ((gamma_s - 1.) * rho * eps - gamma_s * p_inf_s);
     const DiscreteTensorFunction& sigma_d     = chi_solid * 2*mu/sqrt(det(CG))*(CG-1./3. * trace(CG)*I) - p_d * I;
     const DiscreteScalarFunction& aL_d        = sqrt(chi_solid * (2 * mu) / rho + ((1-chi_solid)*gamma_f * p_d + chi_solid * (gamma_s*(p_d+p_inf_s))) / rho);
     const DiscreteScalarFunction& aT_d        = sqrt(chi_solid * mu / rho);
@@ -567,7 +651,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
     if(law == 1){
       return _apply_NeoHook(rho_d,u_d,E_d,CG_d,chi_solid_d,lambda,mu,gamma_f,p_inf_f);
     } else {
-      return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid_d,lambda,gamma_f,p_inf_f,gamma_s,p_inf_s);
+      return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid_d,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
     }
   }
 
@@ -584,7 +668,8 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
     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");
-    }
+    }  //      - p* identity; // - chi_solid*P_zero_air_repo*identity;
+
 
     if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
       throw NormalError("hyperelastic solver expects P0 functions");
@@ -608,100 +693,109 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
     PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
                                                                 symmetry_boundary_descriptor_list);
 
-    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(mesh_v, 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);
-        }
-    }
-
-    NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-
-    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;
-            }
-          }
-        }
-      }
-    );
-
-    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];
-
-            for(size_t R = 0; R < cell_nodes.size(); ++R){
-                const NodeId r = cell_nodes[R];
-
-                for(size_t l = 0; l < sigma_coef.size(); ++l){
-                    const size_t& i = row[l];
-                    const size_t& k = col[l];
-
-                    auto coefficients = sigma_coef[l].coefficients(j);
-
-                    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];
-                    }
-                    
-                    sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
-                    if(i != k){
-                      sigma_lim(j,R)(i,k) *= 2;
-                    }
-                }
-
-                sigma_lim(j,R) += transpose(sigma_lim(j,R));
-                sigma_lim(j,R) *= 0.5;
-            }
-
-        }
-    );
-
-    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
-    auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
-
-    DiscreteFunctionDPk<Dimension, Rd> u_lim  = copy(DPk_uh);
+    //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(mesh_v, 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);
+    //    }
+    //}
+    //
+    //NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
+    //const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    //
+    //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;
+    //        }
+    //      }
+    //    }
+    //  }
+    //);
+    //
+    //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];
+    //
+    //        for(size_t R = 0; R < cell_nodes.size(); ++R){
+    //            const NodeId r = cell_nodes[R];
+    //
+    //            for(size_t l = 0; l < sigma_coef.size(); ++l){
+    //                const size_t& i = row[l];
+    //                const size_t& k = col[l];
+    //
+    //                auto coefficients = sigma_coef[l].coefficients(j);
+    //
+    //                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];
+    //                }
+    //                
+    //                sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
+    //                if(i != k){
+    //                  sigma_lim(j,R)(i,k) *= 2;
+    //                }
+    //            }
+    //
+    //            sigma_lim(j,R) += transpose(sigma_lim(j,R));
+    //            sigma_lim(j,R) *= 0.5;
+    //        }
+    //
+    //    }
+    //);
+
+    std::cout << "Reconstruction" << "\n";
+    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v,sigma_v);
+    auto DPk_uh     = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto DPk_sigmah = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>(); 
+    std::cout << "Reconstruction Ok ! " << "\n";
+
+    std::cout << "Limitation" << "\n";
+    DiscreteFunctionDPk<Dimension, Rd> u_lim        = copy(DPk_uh);
+    DiscreteFunctionDPk<Dimension, Rdxd> sigma_lim  = copy(DPk_sigmah);
     this->_vector_limiter(mesh,u,u_lim);
+    this->_tensor_limiter(mesh,sigma,sigma_lim);
+    std::cout << "Limitation OK !" << "\n";
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+    std::cout << "br" << "\n";
     NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
+    std::cout << "br Ok !" << "\n";
 
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -710,7 +804,9 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
     synchronize(br);
 
     NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
+    std::cout << "Fjr" << "\n";
     NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
+    std::cout << "Fjr Ok ! " << "\n";
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -956,14 +1052,14 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
     std::shared_ptr<const DiscreteFunctionVariant> CG_app  = CG;
 
     auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
-    std::tie(m_app,rho_app,u_app,E_app,CG_app) = apply_fluxes(dt/2, rho, u, E, CG, ur, Fjr);
+    //auto [m_app, rho_app, u_app, E_app, CG_app] = apply_fluxes(dt/2., rho, u, E, CG, ur, Fjr);
 
-    std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
+    //std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
 
-    auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
+    //auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
 
-    auto [ur_app, Fjr_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, Fjr_app);
+    //auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
+    return apply_fluxes(dt, rho, u, E, CG, ur, Fjr);
   }
 
   Order2HyperelasticSolver()                     = default;
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index d74835e0b08d7b8d2e6d290057797b5c0f76ffb2..a80ae83fed4234da7fd65f7ef9248406d496fbf5 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -672,6 +672,77 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
         return sigma_lim;
     }
 
+  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+  _apply_NeoHook(const DiscreteScalarFunction rho, 
+                  const DiscreteVectorFunction u, 
+                  const DiscreteScalarFunction E, 
+                  const DiscreteTensorFunction CG, 
+                  const size_t& chi_solid,
+                  const double& lambda,
+                  const double& mu,
+                  const double& gamma,
+                  const double& p_inf) const
+  {
+    const TinyMatrix<Dimension> I = identity;
+
+    const DiscreteScalarFunction& eps         = E - 0.5 * dot(u,u);
+    const DiscreteScalarFunction& p_d         = (1-chi_solid)*(gamma-1) * rho * eps - p_inf*gamma;
+    const DiscreteTensorFunction& sigma_d     = chi_solid * (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I) - p_d * I;
+    const DiscreteScalarFunction& aL_d        = sqrt(chi_solid * (lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
+    const DiscreteScalarFunction& aT_d        = sqrt(chi_solid * mu / rho);
+
+    return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+            std::make_shared<DiscreteFunctionVariant>(aL_d),
+            std::make_shared<DiscreteFunctionVariant>(aT_d)};
+  }
+
+  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+  _apply_NeoHook2(const DiscreteScalarFunction rho, 
+                  const DiscreteVectorFunction u, 
+                  const DiscreteScalarFunction E, 
+                  const DiscreteTensorFunction CG, 
+                  const size_t& chi_solid,
+                  const double& mu,
+                  const double& gamma,
+                  const double& p_inf) const
+  {
+    const TinyMatrix<Dimension> I = identity;
+
+    const DiscreteScalarFunction& eps         = E - 0.5 * dot(u,u);
+    const DiscreteScalarFunction& p_d         = (gamma-1) * rho * eps - p_inf*gamma;
+    const DiscreteTensorFunction& sigma_d     = chi_solid * 2*mu/sqrt(det(CG))*(CG-1./3. * trace(CG)*I) - p_d * I;
+    const DiscreteScalarFunction& aL_d        = sqrt(chi_solid * (2 * mu) / rho + ((1-chi_solid)*gamma * p_d + chi_solid * (gamma*(p_d+p_inf))) / rho);
+    const DiscreteScalarFunction& aT_d        = sqrt(chi_solid * mu / rho);
+
+    return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+            std::make_shared<DiscreteFunctionVariant>(aL_d),
+            std::make_shared<DiscreteFunctionVariant>(aT_d)};
+  }
+
+  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+  _apply_state_law(const size_t law,
+                   const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+                   const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                   const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+                   const std::shared_ptr<const DiscreteFunctionVariant>& CG_v,
+                   const size_t& chi_solid,
+                   const double& lambda,
+                   const double& mu,
+                   const double& gamma,
+                   const double& p_inf) const 
+  {
+    const DiscreteScalarFunction& rho_d       = rho_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u_d         = u_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& E_d         = E_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction& CG_d        = CG_v->get<DiscreteTensorFunction>();
+
+    if(law == 1){
+      return _apply_NeoHook(rho_d,u_d,E_d,CG_d,chi_solid,lambda,mu,gamma,p_inf);
+    } else {
+      return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid,lambda,gamma,p_inf);
+    }
+  }
+
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
@@ -1307,6 +1378,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
              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,
@@ -1319,9 +1391,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                                    NodeValue<Rd> CR_ur,
                                    NodeValuePerCell<Rd> CR_Fjr,
                                    std::vector<NodeId>& map2,
-                                   size_t fluid,
+                                   size_t chi_solid,
                                    const double& lambda,
-                                   const double& mu) const
+                                   const double& mu,
+                                   const double& gamma,
+                                   const double& p_inf) const
   {
     
     std::shared_ptr m_app = getCommonMesh({rho, aL, aT, u, sigma});
@@ -1339,15 +1413,12 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     std::shared_ptr<const DiscreteFunctionVariant> new_aL   = aL;
     std::shared_ptr<const DiscreteFunctionVariant> new_aT   = aT;
 
-    DiscreteScalarFunction new_aL_d = aL->get<DiscreteScalarFunction>();
-    DiscreteScalarFunction new_aT_d = aT->get<DiscreteScalarFunction>();
-
-    const TinyMatrix<Dimension> I = identity;
-    const double& gamma = 1.4;
     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);
@@ -1358,39 +1429,12 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
       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 DiscreteScalarFunction& rho_d    = rho_app->get<DiscreteScalarFunction>();
-      const DiscreteVectorFunction& u_d      = u_app->get<DiscreteVectorFunction>();
-      const DiscreteScalarFunction& E_d      = E_app->get<DiscreteScalarFunction>();
-      const DiscreteScalarFunction& eps      = E_d - 0.5 * dot(u_d,u_d);
-      const DiscreteTensorFunction& CG_d     = CG_app->get<DiscreteTensorFunction>();
-      const DiscreteScalarFunction& p_d      = 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_d * I;
-      const DiscreteScalarFunction& aL_d_app = sqrt((lambda + 2 * mu) / rho_d + gamma * p_d / rho_d);
-      const DiscreteScalarFunction& aT_d_app = sqrt(mu / rho_d);
-
-      const std::shared_ptr<const DiscreteFunctionVariant>& sigma_app =
-        std::make_shared<const DiscreteFunctionVariant>(sigma_d);
-      const std::shared_ptr<const DiscreteFunctionVariant>& aL_app =
-        std::make_shared<const DiscreteFunctionVariant>(aL_d_app);
-      const std::shared_ptr<const DiscreteFunctionVariant>& aT_app =
-        std::make_shared<const DiscreteFunctionVariant>(aT_d_app);
+      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); 
 
-      const DiscreteScalarFunction& new_rho_d    = new_rho->get<DiscreteScalarFunction>();
-      const DiscreteVectorFunction& new_u_d      = new_u->get<DiscreteVectorFunction>();
-      const DiscreteScalarFunction& new_E_d      = new_E->get<DiscreteScalarFunction>();
-      const DiscreteScalarFunction& new_eps      = new_E_d - 0.5 * dot(new_u_d,new_u_d);
-      const DiscreteTensorFunction& new_CG_d     = CG_app->get<DiscreteTensorFunction>();
-      const DiscreteScalarFunction& new_p_d      = fluid*(gamma-1) * new_rho_d * new_eps;
-      const DiscreteTensorFunction& new_sigma_d  = (mu / sqrt(det(new_CG_d)) * (new_CG_d - I) + lambda / sqrt(det(new_CG_d)) * log(sqrt(det(new_CG_d))) * I) - new_p_d * I;
-      new_aL_d     = sqrt((lambda + 2 * mu) / new_rho_d + gamma * new_p_d / new_rho_d);
-      new_aT_d     = sqrt(mu / new_rho_d);
-
-      new_sigma = std::make_shared<const DiscreteFunctionVariant>(new_sigma_d); 
-      new_aL    = std::make_shared<const DiscreteFunctionVariant>(new_aL_d);
-      new_aT    = std::make_shared<const DiscreteFunctionVariant>(new_aT_d);
+      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;
 
@@ -1507,6 +1551,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>>
   apply(const SolverType& solver_type,
+        const size_t& law, 
         const double& dt1,
         const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
         const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
@@ -1524,8 +1569,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
         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
+        const double& lambda_f,
+        const double& mu_f,
+        const double& lambda_s,
+        const double& mu_s,
+        const double& gamma_f,
+        const double& p_inf_f,
+        const double& gamma_s,
+        const double& p_inf_s) const
   {
     std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
     std::shared_ptr m1     = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
@@ -1564,35 +1615,25 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
 
     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
 
-    const double gamma            = 1.4; //à modif
-    const TinyMatrix<Dimension> I = identity;
-    double sum_dt                 = dt2;
+    double sum_dt = dt2;
 
-    size_t fluid = 0;
-    if((mu==0) and (lambda==0)){
-      fluid = 1;
+    size_t chi_solid = 0;
+    if((mu_s!=0) and (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
-      const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
-      const DiscreteVectorFunction& u2_d   = u2_app->get<DiscreteVectorFunction>();
-      const DiscreteScalarFunction& E2_d   = E2_app->get<DiscreteScalarFunction>();
-      const DiscreteScalarFunction& eps2   = E2_d - 0.5 * dot(u2_d, u2_d);
-      const DiscreteTensorFunction& CG2_d  = CG2_app->get<DiscreteTensorFunction>();
-      const DiscreteScalarFunction& p2 = fluid * (gamma - 1) * rho2_d * eps2;
-      const DiscreteTensorFunction& sigma2_d = mu / sqrt(det(CG2_d)) * (CG2_d - I) + lambda / sqrt(det(CG2_d)) * log(sqrt(det(CG2_d))) * I - p2 * I;
-      const DiscreteScalarFunction& aL2_d_app = sqrt((lambda + 2 * mu) / rho2_d + gamma * p2 / rho2_d);
-      const DiscreteScalarFunction& aT2_d_app = sqrt(mu / rho2_d);
-
-      sigma2_app = std::make_shared<const DiscreteFunctionVariant>(sigma2_d);
-      aL2_app    = std::make_shared<const DiscreteFunctionVariant>(aL2_d_app);
-      aT2_app    = std::make_shared<const DiscreteFunctionVariant>(aT2_d_app);
+      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;
+	      dt2 = dt1/2 - sum_dt;
       }
 
       auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, rho2_app, aL2_app, aT2_app, u2_app, sigma2_app,
@@ -1604,40 +1645,15 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
       sum_dt += dt2;
     }
 
-    const DiscreteScalarFunction& rho1_d       = rho1_app->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u1_d         = u1_app->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& E1_d         = E1_app->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& eps1         = E1_d - 0.5 * dot(u1_d,u1_d);
-    const DiscreteTensorFunction& CG1_d        = CG1_app->get<DiscreteTensorFunction>();
-    const DiscreteScalarFunction& p1_d         = (1-fluid)*(gamma-1) * rho1_d * eps1;
-    const DiscreteTensorFunction& sigma1_d     = fluid * (mu / sqrt(det(CG1_d)) * (CG1_d - I) + lambda / sqrt(det(CG1_d)) * log(sqrt(det(CG1_d))) * I) - p1_d * I;
-    const DiscreteScalarFunction& aL1_d        = sqrt(fluid * (lambda + 2 * mu) / rho1_d + gamma * p1_d / rho1_d);
-    const DiscreteScalarFunction& aT1_d        = sqrt(fluid * mu / rho1_d);
-
-    const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_app = std::make_shared<const DiscreteFunctionVariant>(sigma1_d);
-    const std::shared_ptr<const DiscreteFunctionVariant>& aL1_app    = std::make_shared<const DiscreteFunctionVariant>(aL1_d);
-    const std::shared_ptr<const DiscreteFunctionVariant>& aT1_app    = std::make_shared<const DiscreteFunctionVariant>(aT1_d);
-
-    const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u2_d   = u2_app->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& E2_d   = E2_app->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& eps2   = E2_d - 0.5 * dot(u2_d, u2_d);
-    const DiscreteTensorFunction& CG2_d  = CG2_app->get<DiscreteTensorFunction>();
-    const DiscreteScalarFunction& p2 = fluid * (gamma - 1) * rho2_d * eps2;
-    const DiscreteTensorFunction& sigma2_d = mu / sqrt(det(CG2_d)) * (CG2_d - I) + lambda / sqrt(det(CG2_d)) * log(sqrt(det(CG2_d))) * I - p2 * I;
-    const DiscreteScalarFunction& aL2_d_app = sqrt((lambda + 2 * mu) / rho2_d + gamma * p2 / rho2_d);
-    const DiscreteScalarFunction& aT2_d_app = sqrt(mu / rho2_d);
-
-    sigma2_app = std::make_shared<const DiscreteFunctionVariant>(sigma2_d);
-    aL2_app    = std::make_shared<const DiscreteFunctionVariant>(aL2_d_app);
-    aT2_app    = std::make_shared<const DiscreteFunctionVariant>(aT2_d_app);
+    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_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, Fjr1_app);
-    auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,bc_descriptor_list2,
-                                                                                        CR_ur_app,CR_Fjr_app,map2,fluid,lambda,mu);
+    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);
 
     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 58fe0b85242dcf0dc8f6012222d960ed5b59162a..da09ae6138414ed0a8adcf6999b4789d9549be36 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -93,6 +93,7 @@ class Order2LocalDtHyperelasticSolverHandler
                        std::shared_ptr<const DiscreteFunctionVariant>,
                        std::shared_ptr<const DiscreteFunctionVariant>>
   apply(const SolverType& solver_type,
+        const size_t& law,
         const double& dt1,
         const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
         const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
@@ -110,8 +111,14 @@ class Order2LocalDtHyperelasticSolverHandler
         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;
+        const double& lambda_f,
+        const double& mu_f,
+        const double& lambda_s,
+        const double& mu_s,
+        const double& gamma_f,
+        const double& p_inf_f,
+        const double& gamma_s,
+        const double& p_inf_s) const = 0;
 
     IOrder2LocalDtHyperelasticSolver()                                        = default;
     IOrder2LocalDtHyperelasticSolver(IOrder2LocalDtHyperelasticSolver&&)            = default;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2af32e157475b0fe99b9f2cc9cc6c445e05681b3..9c8a7bf1857070e93b3f5665e1a7a24c56314b12 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -41,6 +41,16 @@ symmetrize_coordinates(const TinyVector<Dimension>& origin,
   return u - 2 * dot(u - origin, normal) * normal;
 }
 
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_tensor(const TinyVector<Dimension>& normal,
+                  const TinyMatrix<Dimension>& A)
+{
+  const TinyMatrix<Dimension>& I = identity;
+  const TinyMatrix<Dimension>& S = (I - 2*tensorProduct(normal,normal));
+  return S*A*S;
+}
+
 class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
 {
  public:
@@ -927,12 +937,21 @@ PolynomialReconstruction::_build(
                     } else if constexpr (is_tiny_matrix_v<DataType>) {
                       if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
                                     (DataType::NumberOfColumns == MeshType::Dimension)) {
-                        throw NotImplementedError("TinyMatrix symmetries for reconstruction");
-                      }
-                      const DataType& qi_qj = discrete_function[cell_i_id] - qj;
-                      for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
-                        for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
-                          B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+                        const Rd& normal = symmetry_normal_list[i_symmetry];
+
+                        const DataType& qi = discrete_function[cell_i_id];
+                        const DataType& qi_qj = symmetrize_tensor(normal,qi) - qj;
+                        for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                          for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                            B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+                          }
+                        }
+                      } else {
+                        const DataType& qi_qj = discrete_function[cell_i_id] - qj;
+                        for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                          for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                            B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+                          }
                         }
                       }
                     }