Merge branch 'master' of https://github.com/sofa-framework/sofa

Author: SperryBot

Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/material/Ogden.h

@@ -66,9 +66,7 @@ class Ogden: public HyperelasticMaterial<DataTypes>
                                 const MaterialParameters<DataTypes>& param)
     {
         const MatrixSym& C = sinfo->deformationTensor;
-        const Real mu1 = param.parameterArray[0];
         const Real alpha1 = param.parameterArray[1];
-        const Real k0 = param.parameterArray[2];
 
         m_FJ = pow(sinfo->J, -alpha1/static_cast<Real>(3));
         m_logJ = log(sinfo->J);
@@ -79,17 +77,16 @@ class Ogden: public HyperelasticMaterial<DataTypes>
             for (sofa::Index j = 0; j < 3; ++j) 
                 CEigen(i, j) = C[MatrixSym::voigtID(i, j)];
 
-        // 17/11/2025: Disable /*Eigen::SelfAdjointEigenSolver<EigenMatrix>*/
-        // due to incorrect eigenvector computation for 3x3 matrices.
-        Eigen::EigenSolver<Eigen::Matrix<Real, 3, 3> > EigenProblemSolver(CEigen, true);
+        Eigen::SelfAdjointEigenSolver<Eigen::Matrix<Real,3,3> > 
+            EigenProblemSolver(CEigen, Eigen::ComputeEigenvectors);
         if (EigenProblemSolver.info() != Eigen::Success)
         {
             dmsg_warning("Ogden") << "EigenSolver iterations failed to converge";
             return;
         }
 
-        sinfo->Evalue = EigenProblemSolver.eigenvalues().real().eval();
-        sinfo->Evect = EigenProblemSolver.eigenvectors().real().eval();
+        sinfo->Evalue = EigenProblemSolver.eigenvalues();
+        sinfo->Evect = EigenProblemSolver.eigenvectors();
 
         const Real aBy2{alpha1/static_cast<Real>(2)};
         m_trCaBy2 = static_cast<Real>(0);
@@ -194,8 +191,9 @@ class Ogden: public HyperelasticMaterial<DataTypes>
                     {
                         if (eJ == eI) continue;
 
-                        const bool isDegenerate = std::fabs(Evalue[eI] - Evalue[eJ]) <
-                                                   std::numeric_limits<Real>::epsilon();
+                        constexpr Real tol = static_cast<Real>(1e-8);
+                        const Real maxEval = std::max(std::fabs(Evalue[eI]), std::fabs(Evalue[eJ]));
+                        const bool isDegenerate = std::fabs(Evalue[eI] - Evalue[eJ]) < maxEval * tol;
                         const Real coefRot = isDegenerate 
                             ? aBy2Minus1 * evalPowI2
                             : (evalPowI - pow(Evalue[eJ], aBy2Minus1))/(Evalue[eI] - Evalue[eJ]);