diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index 6875c5005faaf096a2b61a0625a29a9ed389057f..500dce7d7c2411b0859e7397b5b3c2bfb2a747fb 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -325,6 +325,129 @@ 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,
+               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<const Rd>
   _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
   {
@@ -525,6 +648,175 @@ void
     });
   }
 
+  void
+  _vector_limiter(const MeshType& mesh,
+                  const DiscreteFunctionP0<const Rd>& f,
+                  DiscreteFunctionDPk<Dimension, Rd> DPk_fh,
+                  const CellValue<const Rdxd>& grad_u) const 
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+    
+    const auto xj = mesh_data.xj();
+
+    CellValue<Rdxd> sym_grad_u{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId j){
+        sym_grad_u[j] = grad_u[j] + transpose(grad_u[j]);
+    });
+    CellValue<Rdxd> V = _computeEigenVectorsMatrix(mesh,sym_grad_u,1E-4);
+
+    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 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);
+          }
+          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;
+            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);
+            }
+            n[cell_id] = nj;
+            t[cell_id] = tj;
+            l[cell_id] = lj;
+            if(l2Norm(nj) != 0){
+            n[cell_id] = (1/l2Norm(nj))*n[cell_id];
+            }
+            if(l2Norm(tj) != 0){
+              t[cell_id] = (1/l2Norm(tj))*t[cell_id];
+            }
+            if(l2Norm(lj) != 0){
+              l[cell_id] = (1/l2Norm(lj))*l[cell_id];
+            }
+          }
+        }
+      });
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+        const double fn = dot(f[cell_id],n[cell_id]);
+        const double ft = dot(f[cell_id],t[cell_id]);
+        const double fl = dot(f[cell_id],l[cell_id]);
+
+        double fn_min = fn;
+        double fn_max = fn;
+        double ft_min = ft;
+        double ft_max = ft;
+        double fl_min = fl;
+        double fl_max = fl;
+
+        const auto cell_stencil = stencil[cell_id];
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_id]);
+          fn_min            = std::min(fn_min,fn_k);
+          fn_max            = std::max(fn_max,fn_k);
+
+          const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_id]);
+          ft_min            = std::min(ft_min,ft_k);
+          ft_max            = std::max(ft_max,ft_k);
+
+          const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_id]);
+          fl_min            = std::min(fl_min,fl_k);
+          fl_max            = std::max(fl_max,fl_k);
+        }
+
+        double fn_bar_min = fn;
+        double fn_bar_max = fn;
+        double ft_bar_min = ft;
+        double ft_bar_max = ft;
+        double fl_bar_min = fl;
+        double fl_bar_max = fl;
+
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          const CellId cell_k_id = cell_stencil[i_cell];
+          const Rd f_xk          = DPk_fh[cell_id](xj[cell_k_id]);
+
+          const double fn_xk = dot(f_xk,n[cell_id]);
+          fn_bar_min         = std::min(fn_bar_min,fn_xk);
+          fn_bar_max         = std::max(fn_bar_max,fn_xk);
+
+          const double ft_xk = dot(f_xk,t[cell_id]);
+          ft_bar_min         = std::min(ft_bar_min,ft_xk);
+          ft_bar_max         = std::max(ft_bar_max,ft_xk);
+
+          const double fl_xk = dot(f_xk,l[cell_id]);
+          fl_bar_min         = std::min(fl_bar_min,fl_xk);
+          fl_bar_max         = std::max(fl_bar_max,fl_xk);
+        }
+
+        const double eps = 1E-14;
+        double coef1_n   = 1;
+        if(std::abs(fn_bar_max - fn) > eps){
+          coef1_n = (fn_max - fn) / (fn_bar_max - fn);
+        } 
+        double coef2_n = 1;
+        if(std::abs(fn_bar_min - fn) > eps){
+          coef2_n = (fn_min - fn) / (fn_bar_min - fn);
+        }
+        
+        double coef1_t = 1;
+        if(std::abs(ft_bar_max - ft) > eps){
+          coef1_t = (ft_max - ft) / (ft_bar_max - ft);
+        }
+        double coef2_t = 1;
+        if(std::abs(ft_bar_min - ft) > eps){
+          coef2_t = (ft_min - ft) / (ft_bar_min - ft);
+        }
+
+        double coef1_l = 1;
+        if(std::abs(fl_bar_max - fl) > eps){
+          coef1_l = (fl_max - fl) / (fl_bar_max - fl);
+        }
+        double coef2_l = 1;
+        if(std::abs(fl_bar_min - fl) > eps){
+          coef2_l = (fl_min - fl) / (fl_bar_min - fl);
+        }
+
+        const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
+        const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
+        const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
+
+        auto coefficients = DPk_fh.coefficients(cell_id);
+        coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id] 
+                        + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id] 
+                        + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
+        for(size_t i = 1; i < coefficients.size(); ++i){
+          coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id] 
+                          + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id] 
+                          + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
+        }
+    });
+  }
+
   void 
   _vector_limiter(const MeshType& mesh,
                   const DiscreteFunctionP0<const Rd>& f,
@@ -959,7 +1251,9 @@ void
     DiscreteFunctionDPk<Dimension, Rdxd> S_lim     = copy(DPk_Sh);
     DiscreteFunctionDPk<Dimension, double> p_lim   = copy(DPk_ph);
 
-    this->_vector_limiter(mesh, u, u_lim, 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);
+    //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);