From 6a37e06a8c2a40c25e9805d3ba79e133d05bd6f5 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Mon, 18 Nov 2024 17:59:36 +0100
Subject: [PATCH] add exponential function and test
---
src/algebra/EigenvalueSolver.cpp | 11 +++++++++++
src/algebra/EigenvalueSolver.hpp | 2 ++
src/scheme/HyperplasticSolver.cpp | 9 ++++++---
tests/test_EigenvalueSolver.cpp | 14 +++++++++++++-
4 files changed, 32 insertions(+), 4 deletions(-)
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 41e4bf4aa..efaaa19a2 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -516,4 +516,15 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const
return roots;
}
+void
+EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA)
+{
+ auto [eigenvalues, eigenvectors] = findEigen(A);
+ TinyMatrix<3> D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = std::exp(eigenvalues[i]);
+ }
+ expA = eigenvectors * D * transpose(eigenvectors);
+}
+
EigenvalueSolver::EigenvalueSolver() {}
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index 6c01068e0..728cacdfa 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -100,6 +100,8 @@ class EigenvalueSolver
TinyVector<2> _findLastRoots(Polynomial<2> P, const double epsilon) const;
+ void computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA);
+
EigenvalueSolver();
~EigenvalueSolver() = default;
};
diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index e3f06ca9c..9e8a27278 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -203,9 +203,12 @@ 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(TestB);
+ const R3x3 Id = identity;
+ const R3x3 Sd = S - 1. / 3 * trace(S) * Id;
+ // auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestB);
+ R3x3 B;
+ EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, B);
+ std::cout << B << "\n";
exit(0);
double chi = 0;
const double normS = frobeniusNorm(Sd);
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 5efaa921b..d2316d884 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -166,6 +166,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
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> expA2;
+
auto [eigenvalues, eigenmatrix] = EigenvalueSolver{}.findEigen(TestA);
TinyMatrix<3> Diag = zero;
for (size_t i = 0; i < 3; ++i) {
@@ -177,7 +179,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
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) {
@@ -189,6 +190,17 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
REQUIRE((B(i, j) - TestA2(i, j)) / frobeniusNorm(TestA2) == Catch::Approx(0).margin(1E-8));
}
}
+ EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, expA2);
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = std::exp(eigenvalues2[i]);
+ }
+ B = eigenmatrix2 * Diag * transpose(eigenmatrix2);
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(expA2(i, j) - B(i, j) == Catch::Approx(0).margin(1E-12));
+ }
+ }
auto [eigenvalues3, eigenmatrix3] = EigenvalueSolver{}.findEigen(TestB);
Diag = zero;
--
GitLab