diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 40659530f8ddd6b8e26d753891080f63e4d5bd9a..7a6c5999463e66d46ad7b8282f542a9ffe47d936 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -190,4 +190,324 @@ EigenvalueSolver::computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
}
#endif // PUGS_HAS_SLEPC
+template <typename T>
+PUGS_INLINE TinyMatrix<3, 3, T>
+EigenvalueSolver::swap(TinyMatrix<3, 3, T>& matrix, size_t i, size_t j) const
+{
+ TinyMatrix<3, 3, T> swapMat = matrix;
+ for (size_t k = 0; k < 3; k++) {
+ double valuei = matrix(i, k);
+ swapMat(i, k) = matrix(j, k);
+ swapMat(j, k) = valuei;
+ }
+ return swapMat;
+}
+
+template <typename T>
+PUGS_INLINE constexpr TinyMatrix<3, 3, T>
+EigenvalueSolver::rowReduce(const TinyMatrix<3, 3, T>& matrix, const double epsilon) const
+{
+ TinyMatrix<3, 3> redmat = matrix;
+ for (size_t i = 0; i < 3; ++i) {
+ size_t pivotRow = i;
+ while (pivotRow < 3 && std::fabs(redmat(pivotRow, i)) < epsilon) {
+ ++pivotRow;
+ }
+ if (pivotRow == 3) {
+ continue;
+ }
+ std::cout << "before: "
+ << "i = " << i << " pivotRow = " << pivotRow << " " << redmat << "\n";
+ redmat = swap(redmat, i, pivotRow);
+ std::cout << "after: " << redmat << "\n";
+
+ double pivotValue = redmat(i, i);
+ for (size_t j = i; j < 3; ++j) {
+ redmat(i, j) /= pivotValue;
+ }
+
+ for (size_t k = i + 1; k < 3; ++k) {
+ double factor = redmat(k, i);
+ for (size_t j = i; j < 3; ++j) {
+ redmat(k, j) -= factor * redmat(i, j);
+ }
+ }
+ }
+ return redmat;
+}
+
+TinyMatrix<3, 3>
+EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
+ const double eigenvalue,
+ const size_t space_size,
+ const double epsilon) const
+{
+ TinyMatrix<3, 3> eigenvectors = zero;
+ if (space_size == 3) {
+ return eigenvectors0;
+ }
+ TinyMatrix<3, 3> B = A;
+
+ for (size_t i = 0; i < 3; ++i) {
+ B(i, i) -= eigenvalue;
+ }
+
+ std::cout << " Avant " << B << "\n"
+ << "\n";
+
+ B = rowReduce(B, epsilon);
+
+ std::cout << " Après " << B << "\n"
+ << "\n";
+ for (size_t i = 0; i < 3; ++i) {
+ bool isFreeVariableRow = true;
+ for (size_t j = 0; j < 3; ++j) {
+ if (std::fabs(B(i, j)) > epsilon) {
+ isFreeVariableRow = false;
+ std::cout << " B(" << i << "," << j << ") = " << B(i, j) << "\n";
+ break;
+ }
+ }
+
+ if (isFreeVariableRow) {
+ TinyVector<3> eigenvector = zero;
+ eigenvector[i] = 1.0;
+ std::cout << i << " is free "
+ << "\n";
+ std::cout << std::flush;
+ for (int k = i - 1; k >= 0; --k) {
+ double sum = 0.0;
+ for (size_t j = k + 1; j < 3; ++j) {
+ sum += B(k, j) * eigenvector[j];
+ }
+ eigenvector[k] = -sum;
+ }
+
+ eigenvector = 1. / l2Norm(eigenvector) * eigenvector;
+
+ for (size_t j = 0; j < 3; ++j) {
+ eigenvectors(i, j) = eigenvector[j];
+ }
+ }
+ }
+ std::cout << " Before orthorgonalization " << B << "\n";
+ std::cout << " Before orthorgonalization " << eigenvectors << "\n";
+
+ if (space_size == 2) {
+ TinyVector<3> first = zero;
+ TinyVector<3> second = zero;
+ // TinyVector<3> new_second =zero;
+ size_t rank1 = 0;
+ size_t rank2 = 1;
+ for (size_t j = 0; j < 3; ++j) {
+ first[j] = eigenvectors(0, j);
+ second[j] = eigenvectors(1, j);
+ }
+ if (l2Norm(first) < 1e-3) {
+ for (size_t j = 0; j < 3; ++j) {
+ first[j] = eigenvectors(2, j);
+ rank1 = 2;
+ }
+ }
+ if (l2Norm(second) < 1e-3) {
+ for (size_t j = 0; j < 3; ++j) {
+ second[j] = eigenvectors(2, j);
+ rank2 = 2;
+ }
+ }
+ double scalprod = dot(first, second);
+ double norm2 = dot(first, first);
+ TinyVector<3> normal_vector = second - scalprod / norm2 * first;
+ normal_vector = 1. / l2Norm(normal_vector) * normal_vector;
+ for (size_t j = 0; j < 3; ++j) {
+ eigenvectors(rank2, j) = normal_vector[j];
+ }
+ std::cout << " rank1 : " << rank1 << " rank2 " << rank2 << "\n";
+ }
+ std::cout << "In findEigenvectors for eigenvalue " << eigenvalue << " eigenvectors " << eigenvectors << "\n";
+ return eigenvectors;
+}
+
+bool
+EigenvalueSolver::isDiagonal(const TinyMatrix<3, 3>& A, const double epsilon)
+{
+ bool isDiag = false;
+ 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;
+}
+
+TinyMatrix<3, 3>
+EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eigenvalues, const double epsilon) const
+{
+ TinyMatrix<3> eigenMatrix = zero;
+ size_t count = 0;
+ TinyVector<3, int> space_size = zero;
+ for (size_t i = 0; i < 3; i++) {
+ space_size[i] = 1;
+ for (size_t j = i + 1; j < 3; j++) {
+ if (fabs(eigenvalues[i] - eigenvalues[j]) < epsilon) {
+ double temp = eigenvalues[i + 1];
+ eigenvalues[i + 1] = eigenvalues[j];
+ eigenvalues[j] = temp;
+ space_size[i] += 1;
+ i += 1;
+ }
+ }
+ }
+ std::cout << "space_size: " << space_size << "\n";
+ std::cout << "eigenvalues: " << eigenvalues << "\n";
+ size_t it = 0;
+ while (count < 3 && it < 3) {
+ TinyMatrix<3> Try = findEigenvector(A, eigenvalues[count], space_size[count], epsilon);
+ std::cout << eigenvalues[count] << " Frobenius norm try: " << frobeniusNorm(Try) << "\n";
+ if (frobeniusNorm(Try) < 1e-3) {
+ std::cout << " Error : no eigenvector corresponds to eigenvalue " << eigenvalues[count] << "\n";
+ }
+ TinyVector<3> eigenvector = zero;
+ int count2 = 0;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; j++) {
+ eigenvector[j] = Try(i, j);
+ }
+ std::cout << " count, i: " << count << ", " << i << " l2Norm(eigenvector) " << l2Norm(eigenvector) << "\n";
+ if (l2Norm(eigenvector) > 1e-3) {
+ for (size_t j = 0; j < 3; j++) {
+ eigenMatrix(count, j) = eigenvector[j];
+ }
+ count += 1;
+ count2 += 1;
+ }
+ }
+ if (count2 != space_size[count - count2]) {
+ std::cout << " error eigenvalue " << eigenvalues[count - 1] << " eigenvector space size " << count2
+ << " is not what was expected " << space_size[count - 1] << "\n";
+ }
+ it++;
+ }
+ std ::cout << " eigenMatrix " << eigenMatrix << "\n";
+ return eigenMatrix;
+}
+
+std::tuple<TinyVector<3>, TinyMatrix<3>>
+EigenvalueSolver::findEigen(const TinyMatrix<3> A)
+{
+ const double epsilon = 1.e-6; // * l2Norm(eigenvalues);
+ if (frobeniusNorm(A) < 1.e-15) {
+ return std::make_tuple(eigenvalues0, eigenvectors0);
+ }
+ TinyMatrix<3> C = 1. / frobeniusNorm(A) * A;
+ if (isDiagonal(C, epsilon)) {
+ 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) {
+ Diag(i, i) = eigenvalues[i];
+ }
+ TinyMatrix<3> B = eigenvectors0 * Diag * transpose(eigenvectors0);
+ std::cout << "\n"
+ << B << ", " << A << " Diff "
+ << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
+ return std::make_tuple(eigenvalues, eigenvectors0);
+ }
+ const TinyVector<3> eigenvalues = _findEigenValues(C, epsilon);
+ // std::cout << std::setprecision(15) << eigenvalues << "\n";
+ TinyMatrix<3> Eigenmatrix = transpose(findEigenvectors(C, eigenvalues, epsilon));
+ TinyMatrix<3> Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = frobeniusNorm(A) * eigenvalues[i];
+ }
+ TinyMatrix<3> B = Eigenmatrix * Diag * transpose(Eigenmatrix);
+ std::cout << "\n"
+ << B << ", " << A << " Diff "
+ << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
+ return std::make_tuple(eigenvalues, Eigenmatrix);
+}
+TinyVector<3>
+EigenvalueSolver::_findEigenValues(const TinyMatrix<3>& A, const double epsilon) const
+{
+ std::cout << " Matrix " << A << "\n";
+ 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) > 1e-11) {
+ eigenvalues[0] = _findFirstRoot(P);
+ Polynomial<1> Q1(-eigenvalues[0], 1);
+ Polynomial<3> S;
+ Polynomial<3> R;
+ divide(P, Q1, S, R);
+ std::cout << "Q1 : " << Q1 << "\n";
+ std::cout << "R : " << R << "\n";
+ std::cout << "S : " << 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, epsilon);
+ eigenvalues[1] = last_eigenvalues[0];
+ eigenvalues[2] = last_eigenvalues[1];
+ TinyVector<3, int> space_size = zero;
+ for (size_t i = 0; i < 3; i++) {
+ space_size[i] = 1;
+ for (size_t j = i + 1; j < 3; j++) {
+ if (eigenvalues[i] == eigenvalues[j]) {
+ double temp = eigenvalues[i + 1];
+ eigenvalues[i + 1] = eigenvalues[j];
+ eigenvalues[j] = temp;
+ space_size[i] += 1;
+ i += 1;
+ }
+ }
+ }
+ return eigenvalues;
+}
+
+double
+EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
+{
+ double guess = -P.coefficients()[2] / 3.;
+ double old_guess = -P.coefficients()[2] / 3.;
+ std::cout << " coefs P " << P.coefficients() << "\n";
+ // double delta = pow(2.*P.coefficients()[2],2) - 4*3.*P.coefficients()[3]*P.coefficients()[1];
+ Polynomial<2> Q = derivative(P);
+ Polynomial<1> R = derivative(Q);
+ size_t iteration = 0;
+ double residu = 0.;
+ do {
+ // guess -= P(guess) / Q(guess);
+ old_guess = guess;
+ guess -= 2 * P(guess) * Q(guess) / (2 * Q(guess) * Q(guess) - P(guess) * R(guess));
+ std::cout << "guess = " << guess << "\n";
+ std::cout << "g(guess) = " << Q(guess) << "\n";
+ residu = P(guess);
+ std::cout << "residu = " << residu << "\n";
+ iteration += 1;
+ } while ((fabs(residu) > 1.e-16) or (fabs(old_guess - guess) > 1.e-10));
+ std::cout << " nb Newton iterations " << iteration << "\n";
+ return guess;
+}
+
+TinyVector<2>
+EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const
+{
+ TinyVector<2> roots(0., 0.);
+ double delta = pow(P.coefficients()[1], 2) - 4 * P.coefficients()[2] * P.coefficients()[0];
+ Assert(delta > -epsilon, "Find roots is only for symmetric matrix");
+ if (fabs(delta) < epsilon) {
+ roots[0] = -P.coefficients()[1] / (2. * P.coefficients()[2]);
+ roots[1] = roots[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;
+}
+
EigenvalueSolver::EigenvalueSolver() {}
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index febea024a94b7a3b1e033005bf70003a55b9c293..6c01068e0c04d7fd89e97d2d45d8117367045934 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -5,6 +5,7 @@
#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
+#include <analysis/Polynomial.hpp>
#include <utils/Exceptions.hpp>
#include <utils/SmallArray.hpp>
@@ -12,6 +13,8 @@ class EigenvalueSolver
{
private:
struct Internals;
+ TinyMatrix<3> eigenvectors0{1, 0, 0, 0, 1, 0, 0, 0, 1};
+ TinyVector<3> eigenvalues0{0, 0, 0};
#ifdef PUGS_HAS_SLEPC
void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues);
@@ -75,6 +78,28 @@ class EigenvalueSolver
#endif // PUGS_HAS_SLEPC
}
+ template <typename T>
+ PUGS_INLINE TinyMatrix<3, 3, T> swap(TinyMatrix<3, 3, T>& matrix, size_t i, size_t j) const;
+
+ template <typename T>
+ PUGS_INLINE constexpr TinyMatrix<3, 3, T> rowReduce(const TinyMatrix<3, 3, T>& matrix, const double epsilon) const;
+
+ TinyMatrix<3, 3> findEigenvector(const TinyMatrix<3, 3>& A,
+ const double eigenvalue,
+ const size_t space_size,
+ const double epsilon) const;
+
+ bool isDiagonal(const TinyMatrix<3, 3>& A, const double epsilon);
+
+ TinyMatrix<3, 3> findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eigenvalues, const double epsilon) const;
+
+ std::tuple<TinyVector<3>, TinyMatrix<3>> findEigen(const TinyMatrix<3> A);
+ TinyVector<3> _findEigenValues(const TinyMatrix<3>& A, const double epsilon) const;
+
+ double _findFirstRoot(Polynomial<3> P) const;
+
+ TinyVector<2> _findLastRoots(Polynomial<2> P, const double epsilon) const;
+
EigenvalueSolver();
~EigenvalueSolver() = default;
};
diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index 0b422ed740b6ea34251a0c2966b0179467ee3e46..a69cd8ab9bf5f7e1eca3d507e4c18c10b1fd737a 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -200,306 +200,12 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
return expA;
}
- template <typename T>
- PUGS_INLINE TinyMatrix<3, 3, T>
- swap(TinyMatrix<3, 3, T>& matrix, size_t i, size_t j) const
- {
- TinyMatrix<3, 3, T> swapMat = matrix;
- for (size_t k = 0; k < 3; k++) {
- double valuei = matrix(i, k);
- swapMat(i, k) = matrix(j, k);
- swapMat(j, k) = valuei;
- }
- // std::cout << "In swap " << swapMat << "\n";
- return swapMat;
- }
- template <typename T>
- PUGS_INLINE constexpr TinyMatrix<3, 3, T>
- rowReduce(const TinyMatrix<3, 3, T>& matrix, const double epsilon) const
- {
- TinyMatrix<3, 3> redmat = matrix;
- for (size_t i = 0; i < 3; ++i) {
- size_t pivotRow = i;
- while (pivotRow < 3 && std::fabs(redmat(pivotRow, i)) < epsilon) {
- ++pivotRow;
- }
- if (pivotRow == 3) {
- continue;
- }
- std::cout << "before: "
- << "i = " << i << " pivotRow = " << pivotRow << " " << redmat << "\n";
- redmat = swap(redmat, i, pivotRow);
- std::cout << "after: " << redmat << "\n";
-
- double pivotValue = redmat(i, i);
- for (size_t j = i; j < 3; ++j) {
- redmat(i, j) /= pivotValue;
- }
-
- for (size_t k = i + 1; k < 3; ++k) {
- double factor = redmat(k, i);
- for (size_t j = i; j < 3; ++j) {
- redmat(k, j) -= factor * redmat(i, j);
- }
- }
- }
- return redmat;
- }
-
- TinyMatrix<3, 3>
- findEigenvector(const TinyMatrix<3, 3>& A,
- const double eigenvalue,
- const size_t space_size,
- const double epsilon) const
- {
- TinyMatrix<3, 3> eigenvectors = zero;
- if (space_size == 3) {
- return eigenvectors0;
- }
- TinyMatrix<3, 3> B = A;
-
- for (size_t i = 0; i < 3; ++i) {
- B(i, i) -= eigenvalue;
- }
-
- std::cout << " Avant " << B << "\n"
- << "\n";
-
- B = rowReduce(B, epsilon);
-
- std::cout << " Après " << B << "\n"
- << "\n";
- for (size_t i = 0; i < 3; ++i) {
- bool isFreeVariableRow = true;
- for (size_t j = 0; j < 3; ++j) {
- if (std::fabs(B(i, j)) > epsilon) {
- isFreeVariableRow = false;
- std::cout << " B(" << i << "," << j << ") = " << B(i, j) << "\n";
- break;
- }
- }
-
- if (isFreeVariableRow) {
- TinyVector<3> eigenvector = zero;
- eigenvector[i] = 1.0;
- std::cout << i << " is free "
- << "\n";
- std::cout << std::flush;
- for (int k = i - 1; k >= 0; --k) {
- double sum = 0.0;
- for (size_t j = k + 1; j < 3; ++j) {
- sum += B(k, j) * eigenvector[j];
- }
- eigenvector[k] = -sum;
- }
-
- eigenvector = 1. / l2Norm(eigenvector) * eigenvector;
-
- for (size_t j = 0; j < 3; ++j) {
- eigenvectors(i, j) = eigenvector[j];
- }
- }
- }
- std::cout << " Before orthorgonalization " << B << "\n";
- std::cout << " Before orthorgonalization " << eigenvectors << "\n";
-
- if (space_size == 2) {
- TinyVector<3> first = zero;
- TinyVector<3> second = zero;
- // TinyVector<3> new_second =zero;
- size_t rank1 = 0;
- size_t rank2 = 1;
- for (size_t j = 0; j < 3; ++j) {
- first[j] = eigenvectors(0, j);
- second[j] = eigenvectors(1, j);
- }
- if (l2Norm(first) < 1e-3) {
- for (size_t j = 0; j < 3; ++j) {
- first[j] = eigenvectors(2, j);
- rank1 = 2;
- }
- }
- if (l2Norm(second) < 1e-3) {
- for (size_t j = 0; j < 3; ++j) {
- second[j] = eigenvectors(2, j);
- rank2 = 2;
- }
- }
- double scalprod = dot(first, second);
- double norm2 = dot(first, first);
- TinyVector<3> normal_vector = second - scalprod / norm2 * first;
- normal_vector = 1. / l2Norm(normal_vector) * normal_vector;
- for (size_t j = 0; j < 3; ++j) {
- eigenvectors(rank2, j) = normal_vector[j];
- }
- std::cout << " rank1 : " << rank1 << " rank2 " << rank2 << "\n";
- }
- std::cout << "In findEigenvectors for eigenvalue " << eigenvalue << " eigenvectors " << eigenvectors << "\n";
- return eigenvectors;
- }
-
- bool
- isDiagonal(const TinyMatrix<3, 3>& A, const double epsilon) const
- {
- bool isDiag = false;
- 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;
- }
-
- TinyMatrix<3, 3>
- findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eigenvalues, const double epsilon) const
- {
- TinyMatrix<3> eigenMatrix = zero;
- size_t count = 0;
- TinyVector<3, int> space_size = zero;
- // eigenvalues[0] = -1;
- // eigenvalues[1] = -1;
- // eigenvalues[2] = 8;
- for (size_t i = 0; i < 3; i++) {
- space_size[i] = 1;
- for (size_t j = i + 1; j < 3; j++) {
- if (fabs(eigenvalues[i] - eigenvalues[j]) < epsilon) {
- double temp = eigenvalues[i + 1];
- eigenvalues[i + 1] = eigenvalues[j];
- eigenvalues[j] = temp;
- space_size[i] += 1;
- i += 1;
- }
- }
- }
- std::cout << "space_size: " << space_size << "\n";
- std::cout << "eigenvalues: " << eigenvalues << "\n";
- size_t it = 0;
- while (count < 3 && it < 3) {
- TinyMatrix<3> Try = findEigenvector(A, eigenvalues[count], space_size[count], epsilon);
- std::cout << eigenvalues[count] << " Frobenius norm try: " << frobeniusNorm(Try) << "\n";
- if (frobeniusNorm(Try) < 1e-3) {
- std::cout << " Error : no eigenvector corresponds to eigenvalue " << eigenvalues[count] << "\n";
- }
- TinyVector<3> eigenvector = zero;
- int count2 = 0;
- for (size_t i = 0; i < 3; ++i) {
- for (size_t j = 0; j < 3; j++) {
- eigenvector[j] = Try(i, j);
- }
- std::cout << " count, i: " << count << ", " << i << " l2Norm(eigenvector) " << l2Norm(eigenvector) << "\n";
- if (l2Norm(eigenvector) > 1e-3) {
- for (size_t j = 0; j < 3; j++) {
- eigenMatrix(count, j) = eigenvector[j];
- }
- count += 1;
- count2 += 1;
- }
- }
- if (count2 != space_size[count - count2]) {
- std::cout << " error eigenvalue " << eigenvalues[count - 1] << " eigenvector space size " << count2
- << " is not what was expected " << space_size[count - 1] << "\n";
- }
- it++;
- }
- std ::cout << " eigenMatrix " << eigenMatrix << "\n";
- return eigenMatrix;
- }
-
- std::tuple<TinyVector<3>, TinyMatrix<3>>
- _testMatrix(const TinyMatrix<3> A) const
- {
- std::cout << std::setprecision(15) << A << "\n";
- const double epsilon = 1.e-6; // * l2Norm(eigenvalues);
- if (frobeniusNorm(A) < 1.e-15) {
- return std::make_tuple(eigenvalues0, eigenvectors0);
- }
- TinyMatrix<3> C = 1. / frobeniusNorm(A) * A;
- if (isDiagonal(C, epsilon)) {
- 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) {
- Diag(i, i) = eigenvalues[i];
- }
- TinyMatrix<3> B = eigenvectors0 * Diag * transpose(eigenvectors0);
- std::cout << "\n"
- << B << ", " << A << " Diff "
- << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
- return std::make_tuple(eigenvalues, eigenvectors0);
- }
- const TinyVector<3> eigenvalues = _findEigenValues(C, epsilon);
- // std::cout << std::setprecision(15) << eigenvalues << "\n";
- TinyMatrix<3> Eigenmatrix = transpose(findEigenvectors(C, eigenvalues, epsilon));
- TinyMatrix<3> Diag = zero;
- for (size_t i = 0; i < 3; ++i) {
- Diag(i, i) = frobeniusNorm(A) * eigenvalues[i];
- }
- TinyMatrix<3> B = Eigenmatrix * Diag * transpose(Eigenmatrix);
- std::cout << "\n"
- << B << ", " << A << " Diff "
- << " A - B " << 1 / frobeniusNorm(A) * (A - B) << "\n";
- return std::make_tuple(eigenvalues, Eigenmatrix);
- }
-
double
_compute_chi(const R3x3 S, const double yield) const
{
const R3x3 Id = identity;
const R3x3 Sd = S - 1. / 3 * trace(S) * Id;
- // std::cout << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "Test A"
- // << "\n";
- // auto [eigenvalues, eigenmatrix] = _testMatrix(TestA);
- std::cout << "\n"
- << "\n"
- << "\n"
- << "\n"
- << "\n"
- << "\n"
- << "Test A2"
- << "\n";
- auto [eigenvalues, eigenmatrix] = _testMatrix(TestA2);
- // std::cout << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "Test B"
- // << "\n";
- // auto [eigenvalues2, eigenmatrix2] = _testMatrix(TestB);
- // std::cout << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "Test C"
- // << "\n";
- // auto [eigenvalues3, eigenmatrix3] = _testMatrix(TestC);
- // std::cout << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "Test D" << TestD << "\n";
- // auto [eigenvalues4, eigenmatrix4] = _testMatrix(TestD);
- // std::cout << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "\n"
- // << "Test E"
- // << "\n";
- // auto [eigenvalues5, eigenmatrix5] = _testMatrix(TestE);
- exit(0);
-
+ // auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestA);
double chi = 0;
const double normS = frobeniusNorm(Sd);
if ((normS - std::sqrt(2. / 3) * yield) > 0.) {
@@ -508,91 +214,6 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
return chi;
}
- TinyVector<3>
- _findEigenValues(const R3x3& A, const double epsilon) const
- {
- std::cout << " Matrix " << A << "\n";
- 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) > 1e-11) {
- eigenvalues[0] = _findFirstRoot(P);
- Polynomial<1> Q1(-eigenvalues[0], 1);
- Polynomial<3> S;
- Polynomial<3> R;
- divide(P, Q1, S, R);
- std::cout << "Q1 : " << Q1 << "\n";
- std::cout << "R : " << R << "\n";
- std::cout << "S : " << 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, epsilon);
- eigenvalues[1] = last_eigenvalues[0];
- eigenvalues[2] = last_eigenvalues[1];
- TinyVector<3, int> space_size = zero;
- for (size_t i = 0; i < 3; i++) {
- space_size[i] = 1;
- for (size_t j = i + 1; j < 3; j++) {
- if (eigenvalues[i] == eigenvalues[j]) {
- double temp = eigenvalues[i + 1];
- eigenvalues[i + 1] = eigenvalues[j];
- eigenvalues[j] = temp;
- space_size[i] += 1;
- i += 1;
- }
- }
- }
- return eigenvalues;
- }
-
- double
- _findFirstRoot(Polynomial<3> P) const
- {
- double guess = -P.coefficients()[2] / 3.;
- double old_guess = -P.coefficients()[2] / 3.;
- std::cout << " coefs P " << P.coefficients() << "\n";
- // double delta = pow(2.*P.coefficients()[2],2) - 4*3.*P.coefficients()[3]*P.coefficients()[1];
- Polynomial<2> Q = derivative(P);
- Polynomial<1> R = derivative(Q);
- size_t iteration = 0;
- double residu = 0.;
- do {
- // guess -= P(guess) / Q(guess);
- old_guess = guess;
- guess -= 2 * P(guess) * Q(guess) / (2 * Q(guess) * Q(guess) - P(guess) * R(guess));
- std::cout << "guess = " << guess << "\n";
- std::cout << "g(guess) = " << Q(guess) << "\n";
- residu = P(guess);
- std::cout << "residu = " << residu << "\n";
- iteration += 1;
- } while ((fabs(residu) > 1.e-16) or (fabs(old_guess - guess) > 1.e-10));
- std::cout << " nb Newton iterations " << iteration << "\n";
- return guess;
- }
-
- TinyVector<2>
- _findLastRoots(Polynomial<2> P, const double epsilon) const
- {
- TinyVector<2> roots(0., 0.);
- double delta = pow(P.coefficients()[1], 2) - 4 * P.coefficients()[2] * P.coefficients()[0];
- Assert(delta > -epsilon, "Find roots is only for symmetric matrix");
- if (fabs(delta) < epsilon) {
- roots[0] = -P.coefficients()[1] / (2. * P.coefficients()[2]);
- roots[1] = roots[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
{
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 323d2d6a54058359f1d004ef2f0cf60f7f57f2af..248af664b3ca26c294afed19fe45379fde99a209 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -6,6 +6,7 @@
#include <algebra/EigenvalueSolver.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyMatrix.hpp>
+#include <scheme/HyperelasticSolver.hpp>
#include <utils/pugs_config.hpp>
@@ -35,24 +36,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
CRSMatrix A{S.getCRSMatrix()};
- // Array<int> non_zeros2(3);
- // non_zeros2.fill(3);
- // CRSMatrixDescriptor<double> S2{3, 3, non_zeros};
-
- // S2(0, 0) = 1;
- // S2(0, 1) = 0;
- // S2(0, 2) = 0;
-
- // S2(1, 0) = 0;
- // S2(1, 1) = 0;
- // S2(1, 2) = 0;
-
- // S2(2, 0) = 0;
- // S2(2, 1) = 0;
- // S2(2, 2) = 1;
-
- // CRSMatrix B{S2.getCRSMatrix()};
-
SECTION("eigenvalues")
{
SmallArray<double> eigenvalues;
@@ -125,7 +108,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
D(i, i) = eigenvalue_list[i];
}
SmallMatrix PT = transpose(P);
- // double eigenval = -1.;
TinyMatrix<3, 3, double> TinyA;
for (size_t i = 0; i < 3; ++i) {
@@ -133,7 +115,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
TinyA(i, j) = S(i, j);
}
}
- // auto E = findEigenvectors(TinyA, eigenval);
SmallMatrix PDPT = P * D * PT;
for (size_t i = 0; i < 3; ++i) {
@@ -150,8 +131,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
{
SmallMatrix<double> expA{};
SmallMatrix<double> expS{3};
- // SmallMatrix<double> expB{};
- // SmallMatrix<double> expS2{3};
expS(0, 0) = 1325.074593930307812;
expS(0, 1) = 662.353357244568285;
@@ -172,19 +151,23 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
}
}
- // EigenvalueSolver{}.computeExpForSymmetricMatrix(B, expB);
-
- // for (size_t i = 0; i < 3; ++i) {
- // for (size_t j = 0; j < 3; ++j) {
- // REQUIRE(expB(i, j) - expS(i, j) == Catch::Approx(0).margin(1E-12));
- // }
- // }
-
#else // PUGS_HAS_SLEPC
REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA),
"not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif // PUGS_HAS_SLEPC
}
}
+ SECTION("symmetric tiny matrix")
+ {
+ TinyMatrix<3> TestA{3e10, 2e10, 4e10, 2e10, 0, 2e10, 4e10, 2e10, 3e10};
+ TinyMatrix<3> TestA2{4, 2, 3, 2, 0, 2, 3, 2, 4};
+ TinyMatrix<3> TestB{1, -1, 0, -1, 1, 0, 0, 0, 3};
+ 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);
+ }
}
}