diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index 500dce7d7c2411b0859e7397b5b3c2bfb2a747fb..96f1db1a0872e31e481c651f0a90e7533d1dfdf4 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -1,3 +1,4 @@
+#include <algebra/EigenvalueSolver.hpp>
 #include <scheme/Order2HyperelasticSolver.hpp>
 
 #include <language/utils/InterpolateItemValue.hpp>
@@ -667,7 +668,14 @@ void
       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<SmallArray<double>> eigunvalues{mesh.connectivity()};
+    CellValue<std::vector<SmallVector<double>>> eigunvectors{mesh.connectivity()};
+    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
+      EigenvalueSolver{}.computeForSymmetricMatrix(sym_grad_u[j],eigunvalues[j],eigunvectors[j]);
+    }
+    
+    //CellValue<Rdxd> V = _computeEigenVectorsMatrix(mesh,sym_grad_u,1E-4);
 
     CellValue<Rd> n{mesh.connectivity()};
     CellValue<Rd> t{mesh.connectivity()};
@@ -677,13 +685,16 @@ void
         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;
-          Rdxd Vj = V[cell_id];
+          //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] = Vj(i,0);
+            // tj[i] = Vj(i,1);
+            nj[i] = eigunvectors[cell_id][0][i];
+            tj[i] = eigunvectors[cell_id][1][i];
           }
           n[cell_id] = nj;
           t[cell_id] = tj;
@@ -699,11 +710,14 @@ void
             Rd nj;
             Rd tj;
             Rd lj;
-            Rdxd Vj = V[cell_id];
+        //  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] = 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];
             }
             n[cell_id] = nj;
             t[cell_id] = tj;
@@ -719,8 +733,9 @@ void
             }
           }
         }
-      });
-
+        }
+    });
+    
     parallel_for(
       mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
         const double fn = dot(f[cell_id],n[cell_id]);