code the eigenvalues calculation

src/scheme/HyperplasticSolver.cpp

-#include <scheme/HyperplasticSolver.hpp>
-
 #include <algebra/CRSMatrix.hpp>
 #include <algebra/CRSMatrixDescriptor.hpp>
 #include <algebra/EigenvalueSolver.hpp>
 #include <algebra/SmallMatrix.hpp>
+#include <analysis/Polynomial.hpp>
+#include <scheme/HyperplasticSolver.hpp>
 
 #include <language/utils/InterpolateItemValue.hpp>
 #include <mesh/ItemValueUtils.hpp>
@@ -292,22 +292,81 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
   double
   _compute_chi(const R3x3 S, const double yield) const
   {
+    const R3x3 Id = identity;
+    const R3x3 Sd = S - 1. / 3 * trace(S) * Id;
     double vp     = 2;
     TinyMatrix<3> A{1, -1, 0, -1, 1, 0, 0, 0, 3};
     std::cout << A << "\n";
+    const TinyVector<3> eigenvalues = _findEigenValues(Sd);
+    std::cout << eigenvalues << "\n";
     TinyMatrix<3> Try = findEigenvectors(A, vp);
     std::cout << Try << "\n";
     exit(0);
-    const R3x3 Id      = identity;
-    const R3x3 Sbar    = S - 1. / 3 * trace(S) * Id;
     double chi         = 0;
-    const double normS = frobeniusNorm(Sbar);
+    const double normS = frobeniusNorm(Sd);
     if ((normS - std::sqrt(2. / 3) * yield) > 0.) {
       chi = 1;
     }
     return chi;
   }
 
+  TinyVector<3>
+  _findEigenValues(const R3x3& A) const
+  {
+    TinyVector<3> eigenvalues(0., 0., 0.);
+    double b = -trace(A);
+    double c = 0.5 * (trace(A) * trace(A) - trace(A * A));
+    double d = -det(A);
+    Polynomial<3> P(d, c, b, 1.);
+    std::cout << P << "\n";
+    Polynomial<2> Q(c, b, 1.);
+    if (fabs(d) > 1.e-10) {
+      eigenvalues[0] = _findFirstRoot(P);
+      Polynomial<1> Q1(eigenvalues[0], 1);
+      Polynomial<3> S;
+      Polynomial<3> R;
+      divide(P, Q1, S, R);
+      std::cout << "S"
+                << "\n";
+      Q.coefficients()[0] = S.coefficients()[0];
+      Q.coefficients()[1] = S.coefficients()[1];
+      Q.coefficients()[2] = S.coefficients()[2];
+    }
+    TinyVector<2> last_eigenvalues = _findLastRoots(Q);
+    eigenvalues[1]                 = last_eigenvalues[0];
+    eigenvalues[2]                 = last_eigenvalues[1];
+    return eigenvalues;
+  }
+
+  double
+  _findFirstRoot(Polynomial<3> P) const
+  {
+    double guess = -P.coefficients()[2] / 3.;
+    //    double delta = pow(2.*P.coefficients()[2],2) - 4*3.*P.coefficients()[3]*P.coefficients()[1];
+    Polynomial<2> Q  = derivative(P);
+    size_t iteration = 0;
+    double residu    = 0;
+    while (residu > 1.e-14) {
+      guess -= P(guess) / Q(guess);
+      residu = P(guess);
+      iteration += 1;
+    }
+    std::cout << " nb Newton iterations " << iteration;
+    return guess;
+  }
+
+  TinyVector<2>
+  _findLastRoots(Polynomial<2> P) const
+  {
+    TinyVector<2> roots(0., 0.);
+    double delta = pow(P.coefficients()[1], 2) - 4 * P.coefficients()[2] * P.coefficients()[0];
+    if (delta >= 0.) {
+      roots[0] = (-P.coefficients()[1] - std::sqrt(delta)) / (2. * P.coefficients()[2]);
+      roots[1] = (-P.coefficients()[1] + std::sqrt(delta)) / (2. * P.coefficients()[2]);
+    }
+    return roots;
+  }
+
   NodeValuePerCell<const Rdxd>
   _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) const
   {