Incorrect eigen decomposition in the Ogden material model

[th-skam]

Problem

Lately, I've been working on the Ogden material model which fails the unit tests in Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/Material_test.cpp.

This model uses the SelfAdjointEigenSolver from the Eigen lib to compute eigenpairs of a matrix. The eigenvalues are computed correctly but the eigenvectors it produces are wrong. To progress, I switched to the generic EigenSolver which does the job although it is not the optimal choice for symmetric matrices. I wanted to document this in case someone else comes across it. Maybe it gets fixed with Eigen 5.0 (I am currently using Eigen-3.4.0).

An example

In MaterialTest_cpp, a random tensor is generated, then adjusted to become positive-definite. On my machine, the random matrix is:

$$C = \begin{bmatrix} 6.03838452029702 & -0.500804028823228 & 1.0692913506807 \\ -0.500804028823228 & 7.01386290516793 & 2.22215388836855 \\\ 1.0692913506807 & 2.22215388836855 & 1.53182086321567 \end{bmatrix}$$

For this matrix, the eigenpairs ($\lambda, \vec{v}$) should satisfy the equation $\lvert C \vec{v} - \lambda\vec{v} \rvert = 0$. The SelfAdjointEigenSolver fails the check.

SelfAdjointEigenSolver	EigenSolver
$\lvert C \vec{v_1} - \lambda_1\vec{v_1} \rvert \approx 4.1$	$\lvert C \vec{v_1} - \lambda_1\vec{v_1} \rvert < 1e-15$
$\lvert C \vec{v_2} - \lambda_2\vec{v_2} \rvert \approx 0.7$	$\lvert C \vec{v_2} - \lambda_2\vec{v_2} \rvert < 2e-15$
$\lvert C \vec{v_3} - \lambda_3\vec{v_3} \rvert \approx 7.4$	$\lvert C \vec{v_3} - \lambda_3\vec{v_3} \rvert < 1e-15$

Code to reproduce this

Eigen::Matrix<double,3,3> CEigen;
CEigen(0,0) = 6.03838452029702; 
CEigen(0,1) = -0.500804028823228; 
CEigen(1,0) = -0.500804028823228; 
CEigen(1,1) = 7.01386290516793;
CEigen(1,2) = 2.22215388836855; 
CEigen(2,1) = 2.22215388836855; 
CEigen(2,0) = 1.0692913506807; 
CEigen(0,2) = 1.0692913506807; 
CEigen(2,2) = 1.53182086321567;

std::cout << "C = " << std::endl;
for (int i = 0; i < 3; i++) {for (int j = 0; j < 3; j++) std::cout << std::setprecision(15) << CEigen(i, j) << " "; std::cout << std::endl;}

Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> SelfAdjEigenSolver(CEigen, true);
Eigen::EigenSolver<Eigen::Matrix<double, 3, 3> > EigenSolver(CEigen, true);

// Evaluate ||Av - λv|| should be = 0
for (int i = 0; i < 3; i++)
    std::cout << "||Cv"<< i << " - EigVal" << i << "v" << i << "|| = "
    << "SelfAdjointEigenSolver: "
    << (CEigen * SelfAdjEigenSolver.eigenvectors().col(i) - SelfAdjEigenSolver.eigenvalues()(i) * SelfAdjEigenSolver.eigenvectors().col(i)).norm()
    << " EigenSolver: "
    << (CEigen * EigenSolver.eigenvectors().real().col(i) - EigenSolver.eigenvalues().real()(i) * EigenSolver.eigenvectors().real().col(i)).norm()
    << std::endl;

[alxbilger]

Hi @th-skam,

does it mean the Ogden materiel works and the error comes from the test?

[th-skam]

Hi, @alxbilger

No, currently, this error affects only the Ogden model. There seem to be other inconsistencies in the model which I'm working on fixing.

To me, it seems the tests are set up correctly. In the test file, the SelfAdjointEigenSolver is used only to compute the eigenvalues and the solver does so correctly.