From 3638d7e4f931da3ec5761a54b6632348b5d61602 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Fri, 8 Nov 2024 19:06:58 +0100
Subject: [PATCH] partial commit, eigenvectors calculation
---
src/algebra/EigenvalueSolver.cpp | 1 -
src/algebra/TinyMatrix.hpp | 72 ++++++++++++++++
src/scheme/HyperplasticSolver.cpp | 131 +++++++++++++++++++++++++++---
tests/test_EigenvalueSolver.cpp | 41 ++++++++++
4 files changed, 234 insertions(+), 11 deletions(-)
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 842150bb5..40659530f 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -188,7 +188,6 @@ EigenvalueSolver::computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
expA = P * D * transpose(P);
EPSDestroy(&eps);
}
-
#endif // PUGS_HAS_SLEPC
EigenvalueSolver::EigenvalueSolver() {}
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 12122a874..29c70bc44 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -573,4 +573,76 @@ inverse(const TinyMatrix<3, 3, T>& A)
return A_cofactors_T *= 1. / determinent;
}
+// template <size_t N, typename T>
+// PUGS_INLINE constexpr TinyMatrix<N, N, T>
+// rowReduce(const TinyMatrix<N, N, T>& matrix)
+// {
+// for (size_t i = 0; i < N; ++i) {
+// size_t pivotRow = i;
+// while (pivotRow < N && std::fabs(matrix[pivotRow][i]) < 1e-10) {
+// ++pivotRow;
+// }
+// if (pivotRow == N) {
+// continue;
+// }
+
+// std::swap(matrix[i], matrix[pivotRow]);
+
+// double pivotValue = matrix[i][i];
+// for (size_t j = i; j < N; ++j) {
+// matrix[i][j] /= pivotValue;
+// }
+
+// for (size_t k = i + 1; k < N; ++k) {
+// double factor = matrix[k][i];
+// for (size_t j = i; j < N; ++j) {
+// matrix[k][j] -= factor * matrix[i][j];
+// }
+// }
+// }
+// }
+
+// template <size_t N, typename T>
+// PUGS_INLINE constexpr TinyMatrix<N, N, T>
+// findEigenvectors(const TinyMatrix<N, N, T>& A, double eigenvalue)
+// {
+// TinyMatrix<N, N, T> B = A;
+
+// for (size_t i = 0; i < N; ++i) {
+// B(i, i) -= eigenvalue;
+// }
+
+// rowReduce(B);
+// TinyMatrix<N, N, T> eigenvectors;
+
+// for (int i = 0; i < N; ++i) {
+// bool isFreeVariableRow = true;
+// for (int j = 0; j < N; ++j) {
+// if (std::abs(B[i][j]) > 1e-10) {
+// isFreeVariableRow = false;
+// break;
+// }
+// }
+
+// if (isFreeVariableRow) {
+// TinyVector<N, T> eigenvector(0.0);
+// eigenvector[i] = 1.0;
+
+// for (int k = i - 1; k >= 0; --k) {
+// double sum = 0.0;
+// for (int j = i; j < N; ++j) {
+// sum += B[k][j] * eigenvector[j];
+// }
+// eigenvector[k] = -sum;
+// }
+
+// eigenvector /= L2Norm(eigenvector);
+
+// eigenvectors.push_back(eigenvector);
+// }
+// }
+
+// return eigenvectors;
+// }
+
#endif // TINYMATRIX_HPP
diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index e56c12c39..d1d999b57 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -165,29 +165,139 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
R3x3
matrix_exp(const R3x3 A) const
{
- R3x3 expA = zero;
+ double normA = frobeniusNorm(A);
+ R3x3 expA = identity;
+ if (normA <= 1.e-12) {
+ return expA;
+ }
SmallMatrix<double> expB{3};
- Array<int> non_zeros(3);
- non_zeros.fill(3);
- CRSMatrixDescriptor<double> S{3, 3, non_zeros};
+ SmallMatrix<double> B(3, 3);
+ // B.fill(0.);
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
- S(i, j) = A(i, j);
+ B(i, i) = A(i, j);
}
}
- CRSMatrix B{S.getCRSMatrix()};
+
+ std::cout << B << "\n";
+ std::cout << std::flush;
+
EigenvalueSolver{}.computeExpForSymmetricMatrix(B, expB);
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
expA(i, j) = expB(i, j);
}
}
+ std::cout << expA << "\n";
return expA;
}
+ template <size_t N, typename T>
+ PUGS_INLINE TinyMatrix<N, N, T>
+ swap(TinyMatrix<N, N, T>& matrix, size_t i, size_t j) const
+ {
+ TinyMatrix<N, N, T> swapMat = matrix;
+ for (size_t k = 0; k < N; 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 <size_t N, typename T>
+ PUGS_INLINE constexpr TinyMatrix<N, N, T>
+ rowReduce(const TinyMatrix<N, N, T>& matrix) const
+ {
+ TinyMatrix<N, N> redmat = matrix;
+ for (size_t i = 0; i < N; ++i) {
+ size_t pivotRow = i;
+ while (pivotRow < N && std::fabs(redmat(pivotRow, i)) < 1e-10) {
+ ++pivotRow;
+ }
+ if (pivotRow == N) {
+ 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 < N; ++j) {
+ redmat(i, j) /= pivotValue;
+ }
+
+ for (size_t k = i + 1; k < N; ++k) {
+ double factor = redmat(k, i);
+ for (size_t j = i; j < N; ++j) {
+ redmat(k, j) -= factor * redmat(i, j);
+ }
+ }
+ }
+ return redmat;
+ }
+
+ template <size_t N>
+ PUGS_INLINE constexpr TinyMatrix<N, N>
+ findEigenvectors(const TinyMatrix<N, N>& A, const double eigenvalue) const
+ {
+ TinyMatrix<N, N> B = A;
+
+ for (size_t i = 0; i < N; ++i) {
+ B(i, i) -= eigenvalue;
+ }
+
+ std::cout << " Avant " << B << "\n"
+ << "\n";
+
+ B = rowReduce(B);
+
+ std::cout << " Après " << B << "\n"
+ << "\n";
+ TinyMatrix<N, N> eigenvectors = zero;
+
+ for (size_t i = 0; i < N; ++i) {
+ bool isFreeVariableRow = true;
+ for (size_t j = 0; j < N; ++j) {
+ if (std::fabs(B(i, j)) > 1e-10) {
+ isFreeVariableRow = false;
+ break;
+ }
+ }
+
+ if (isFreeVariableRow) {
+ TinyVector<N> eigenvector = zero;
+ eigenvector[i] = 1.0;
+
+ for (int k = i - 1; k >= 0; --k) {
+ double sum = 0.0;
+ for (size_t j = i; j < N; ++j) {
+ sum += B(k, j) * eigenvector[j];
+ }
+ eigenvector[k] = -sum;
+ }
+
+ eigenvector = 1. / l2Norm(eigenvector) * eigenvector;
+
+ for (size_t j = 0; j < N; ++j) {
+ eigenvectors(i, j) = eigenvector[j];
+ }
+ }
+ }
+
+ return eigenvectors;
+ }
+
double
_compute_chi(const R3x3 S, const double yield) const
{
+ double vp = 2;
+ TinyMatrix<3> A{1, -1, 0, -1, 1, 0, 0, 0, 3};
+ std::cout << A << "\n";
+ TinyMatrix<3> Try = findEigenvectors(A, vp);
+ std::cout << Try << "\n";
+ exit(0);
const R3x3 Id = identity;
const R3x3 Sbar = S - 1. / 3 * trace(S) * Id;
double chi = 0;
@@ -591,7 +701,8 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
CellValue<double> new_zeta = copy(zeta.cellValues());
CellValue<double> new_yield = copy(yield.cellValues());
CellValue<R3x3> sigma_s = copy(sigma.cellValues());
-
+ CellValue<double> chi{mesh.connectivity()};
+ chi.fill(0.);
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
@@ -612,11 +723,11 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
new_Fe[j] += dt_over_Vj * elastic_fluxes;
- double chi = _compute_chi(sigma[j], yield[j]);
- const R3x3 M = dt_over_Vj * chi * new_Fe[j];
+ chi[j] = _compute_chi(sigma[j], yield[j]);
+ const R3x3 M = dt * chi[j] * zeta[j] * new_Fe[j];
new_Fe[j] = new_Fe[j] * matrix_exp(M);
});
-
+ std::cout << "sum_chi " << sum(chi) << "\n";
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 11c630479..0cd3da4a0 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -5,6 +5,7 @@
#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/EigenvalueSolver.hpp>
#include <algebra/SmallMatrix.hpp>
+#include <algebra/TinyMatrix.hpp>
#include <utils/pugs_config.hpp>
@@ -34,6 +35,24 @@ 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;
@@ -73,6 +92,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
}
D(i, i) = eigenvalue_list[i];
}
+ SmallMatrix<double> eigen{3};
+ eigen = zero;
SmallMatrix PDPT = P * D * PT;
for (size_t i = 0; i < 3; ++i) {
@@ -104,6 +125,15 @@ 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) {
+ for (size_t j = 0; j < 3; ++j) {
+ TinyA(i, j) = S(i, j);
+ }
+ }
+ // auto E = findEigenvectors(TinyA, eigenval);
SmallMatrix PDPT = P * D * PT;
for (size_t i = 0; i < 3; ++i) {
@@ -120,6 +150,8 @@ 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;
@@ -139,6 +171,15 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
REQUIRE(expA(i, j) - expS(i, j) == Catch::Approx(0).margin(1E-12));
}
}
+
+ // 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");
--
GitLab