diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 7a6c5999463e66d46ad7b8282f542a9ffe47d936..41e4bf4aaf3d0e010ba967ec229db2e838af28b2 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -244,6 +244,8 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
{
TinyMatrix<3, 3> eigenvectors = zero;
if (space_size == 3) {
+ std::cout << "space_size = 3 "
+ << "\n";
return eigenvectors0;
}
TinyMatrix<3, 3> B = A;
@@ -331,10 +333,11 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
bool
EigenvalueSolver::isDiagonal(const TinyMatrix<3, 3>& A, const double epsilon)
{
- bool isDiag = false;
+ bool isDiag = true;
for (size_t i = 0; i < 2; i++) {
for (size_t j = i + 1; j < 3; j++) {
isDiag = !(fabs(A(i, j)) > epsilon);
+ return isDiag;
}
}
return isDiag;
@@ -397,10 +400,13 @@ EigenvalueSolver::findEigen(const TinyMatrix<3> A)
{
const double epsilon = 1.e-6; // * l2Norm(eigenvalues);
if (frobeniusNorm(A) < 1.e-15) {
+ std::cout << " Frobenius norm 0 " << frobeniusNorm(A) << "\n";
return std::make_tuple(eigenvalues0, eigenvectors0);
}
TinyMatrix<3> C = 1. / frobeniusNorm(A) * A;
if (isDiagonal(C, epsilon)) {
+ std::cout << "Matrix C " << C << " is diagonal "
+ << "\n";
const TinyVector<3> eigenvalues(A(0, 0), A(1, 1), A(2, 2));
TinyMatrix<3> Diag = zero;
for (size_t i = 0; i < 3; ++i) {
@@ -423,7 +429,7 @@ EigenvalueSolver::findEigen(const TinyMatrix<3> A)
std::cout << "\n"
<< B << ", " << A << " Diff "
<< " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
- return std::make_tuple(eigenvalues, Eigenmatrix);
+ return std::make_tuple(frobeniusNorm(A) * eigenvalues, Eigenmatrix);
}
TinyVector<3>
EigenvalueSolver::_findEigenValues(const TinyMatrix<3>& A, const double epsilon) const
diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index a69cd8ab9bf5f7e1eca3d507e4c18c10b1fd737a..e3f06ca9c279a8a4d91ed0ab6733dff29691ae41 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -203,9 +203,10 @@ 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;
- // auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestA);
+ const R3x3 Id = identity;
+ const R3x3 Sd = S - 1. / 3 * trace(S) * Id;
+ auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestB);
+ exit(0);
double chi = 0;
const double normS = frobeniusNorm(Sd);
if ((normS - std::sqrt(2. / 3) * yield) > 0.) {
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 248af664b3ca26c294afed19fe45379fde99a209..5efaa921b891b064679a3166d2d2221666661235 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -165,9 +165,78 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
TinyMatrix<3> TestC{0, 0, 0, 0, 0, 0, 0, 0, 0};
TinyMatrix<3> TestD{1e10, 0, 0, 0, 1, 0, 0, 0, 1};
TinyMatrix<3> TestE{3e-10, 2e-10, 4e-10, 2e-10, 0, 2e-10, 4e-10, 2e-10, 3e-10};
- TinyMatrix<3> eigenvectors0{1, 0, 0, 0, 1, 0, 0, 0, 1};
- TinyVector<3> eigenvalues0{0, 0, 0};
+
auto [eigenvalues, eigenmatrix] = EigenvalueSolver{}.findEigen(TestA);
+ TinyMatrix<3> Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = eigenvalues[i];
+ }
+ TinyMatrix<3> B = eigenmatrix * Diag * transpose(eigenmatrix);
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE((B(i, j) - TestA(i, j)) / frobeniusNorm(TestA) == Catch::Approx(0).margin(1E-8));
+ }
+ }
+
+ auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestA2);
+ Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = eigenvalues2[i];
+ }
+ B = eigenmatrix2 * Diag * transpose(eigenmatrix2);
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE((B(i, j) - TestA2(i, j)) / frobeniusNorm(TestA2) == Catch::Approx(0).margin(1E-8));
+ }
+ }
+
+ auto [eigenvalues3, eigenmatrix3] = EigenvalueSolver{}.findEigen(TestB);
+ Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = eigenvalues3[i];
+ }
+ B = eigenmatrix3 * Diag * transpose(eigenmatrix3);
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE((B(i, j) - TestB(i, j)) / frobeniusNorm(TestB) == Catch::Approx(0).margin(1E-8));
+ }
+ }
+
+ auto [eigenvalues4, eigenmatrix4] = EigenvalueSolver{}.findEigen(TestC);
+ Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = eigenvalues4[i];
+ }
+ B = eigenmatrix4 * Diag * transpose(eigenmatrix4);
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE((B(i, j) - TestC(i, j)) == Catch::Approx(0).margin(1E-8));
+ }
+ }
+
+ auto [eigenvalues5, eigenmatrix5] = EigenvalueSolver{}.findEigen(TestD);
+ Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = eigenvalues5[i];
+ }
+ B = eigenmatrix5 * Diag * transpose(eigenmatrix5);
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE((B(i, j) - TestD(i, j)) / frobeniusNorm(TestD) == Catch::Approx(0).margin(1E-8));
+ }
+ }
+
+ auto [eigenvalues6, eigenmatrix6] = EigenvalueSolver{}.findEigen(TestE);
+ Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = eigenvalues6[i];
+ }
+ B = eigenmatrix6 * Diag * transpose(eigenmatrix6);
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE((B(i, j) - TestE(i, j)) / frobeniusNorm(TestE) == Catch::Approx(0).margin(1E-8));
+ }
+ }
}
}
}