diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index 4b2c337542bec368e42914b0cdd8599c223f60f5..5611f9dbb9f4f7cd82ae63a9febf7868c98ae48d 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -192,45 +192,45 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
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
+ template <typename T>
+ PUGS_INLINE TinyMatrix<3, 3, T>
+ swap(TinyMatrix<3, 3, T>& matrix, size_t i, size_t j) const
{
- TinyMatrix<N, N, T> swapMat = matrix;
- for (size_t k = 0; k < N; k++) {
+ 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";
+ // 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
+ template <typename T>
+ PUGS_INLINE constexpr TinyMatrix<3, 3, T>
+ rowReduce(const TinyMatrix<3, 3, T>& matrix) const
{
- TinyMatrix<N, N> redmat = matrix;
- for (size_t i = 0; i < N; ++i) {
+ TinyMatrix<3, 3> redmat = matrix;
+ for (size_t i = 0; i < 3; ++i) {
size_t pivotRow = i;
- while (pivotRow < N && std::fabs(redmat(pivotRow, i)) < 1e-10) {
+ while (pivotRow < 3 && std::fabs(redmat(pivotRow, i)) < 1e-10) {
++pivotRow;
}
- if (pivotRow == N) {
+ if (pivotRow == 3) {
continue;
}
- std::cout << "before: "
- << "i = " << i << " pivotRow = " << pivotRow << " " << redmat << "\n";
+ // std::cout << "before: "
+ // << "i = " << i << " pivotRow = " << pivotRow << " " << redmat << "\n";
redmat = swap(redmat, i, pivotRow);
- std::cout << "after: " << redmat << "\n";
+ // std::cout << "after: " << redmat << "\n";
double pivotValue = redmat(i, i);
- for (size_t j = i; j < N; ++j) {
+ for (size_t j = i; j < 3; ++j) {
redmat(i, j) /= pivotValue;
}
- for (size_t k = i + 1; k < N; ++k) {
+ for (size_t k = i + 1; k < 3; ++k) {
double factor = redmat(k, i);
- for (size_t j = i; j < N; ++j) {
+ for (size_t j = i; j < 3; ++j) {
redmat(k, j) -= factor * redmat(i, j);
}
}
@@ -238,13 +238,12 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
return redmat;
}
- template <size_t N>
- PUGS_INLINE constexpr TinyMatrix<N, N>
- findEigenvectors(const TinyMatrix<N, N>& A, const double eigenvalue) const
+ TinyMatrix<3, 3>
+ findEigenvector(const TinyMatrix<3, 3>& A, const double eigenvalue, const double space_size) const
{
- TinyMatrix<N, N> B = A;
+ TinyMatrix<3, 3> B = A;
- for (size_t i = 0; i < N; ++i) {
+ for (size_t i = 0; i < 3; ++i) {
B(i, i) -= eigenvalue;
}
@@ -255,11 +254,21 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
std::cout << " Après " << B << "\n"
<< "\n";
- TinyMatrix<N, N> eigenvectors = zero;
-
- for (size_t i = 0; i < N; ++i) {
+ TinyMatrix<3, 3> eigenvectors = zero;
+ if (space_size == 3) {
+ eigenvectors(0, 0) = 1;
+ eigenvectors(1, 0) = 0;
+ eigenvectors(2, 0) = 0;
+ eigenvectors(0, 1) = 0;
+ eigenvectors(1, 1) = 1;
+ eigenvectors(2, 1) = 0;
+ eigenvectors(0, 2) = 0;
+ eigenvectors(1, 2) = 0;
+ eigenvectors(2, 2) = 1;
+ }
+ for (size_t i = 0; i < 3; ++i) {
bool isFreeVariableRow = true;
- for (size_t j = 0; j < N; ++j) {
+ for (size_t j = 0; j < 3; ++j) {
if (std::fabs(B(i, j)) > 1e-10) {
isFreeVariableRow = false;
break;
@@ -267,12 +276,13 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
}
if (isFreeVariableRow) {
- TinyVector<N> eigenvector = zero;
+ TinyVector<3> eigenvector = zero;
eigenvector[i] = 1.0;
-
+ std::cout << i << " is free "
+ << "\n";
for (int k = i - 1; k >= 0; --k) {
double sum = 0.0;
- for (size_t j = i; j < N; ++j) {
+ for (size_t j = k + 1; j < 3; ++j) {
sum += B(k, j) * eigenvector[j];
}
eigenvector[k] = -sum;
@@ -280,28 +290,126 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
eigenvector = 1. / l2Norm(eigenvector) * eigenvector;
- for (size_t j = 0; j < N; ++j) {
+ for (size_t j = 0; j < 3; ++j) {
eigenvectors(i, j) = eigenvector[j];
}
}
}
-
+ if (space_size == 2) {
+ TinyVector<3> first = zero;
+ TinyVector<3> second = zero;
+ // TinyVector<3> new_second =zero;
+ 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-10) {
+ for (size_t j = 0; j < 3; ++j) {
+ first[j] = eigenvectors(2, j);
+ }
+ }
+ if (l2Norm(second) < 1e-10) {
+ for (size_t j = 0; j < 3; ++j) {
+ second[j] = eigenvectors(2, j);
+ }
+ rank2 = 1;
+ }
+ 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 << "In findEigenvectors for eigenvalue " << eigenvalue << " eigenvectors " << eigenvectors << "\n";
return eigenvectors;
}
+ TinyMatrix<3, 3>
+ findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eigenvalues) 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]) < 1e-10) {
+ 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]);
+ std::cout << "Frobenius norm try: " << frobeniusNorm(Try) << "\n";
+ if (frobeniusNorm(Try) < 1e-10) {
+ 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-10) {
+ for (size_t j = 0; j < 3; j++) {
+ eigenMatrix(count, j) = eigenvector[j];
+ }
+ count += 1;
+ count2 += 1;
+ }
+ }
+ if (count2 != space_size[count - 1]) {
+ 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 TinyVector<3> eigenvalues = _findEigenValues(A);
+ // std::cout << std::setprecision(15) << eigenvalues << "\n";
+ TinyMatrix<3> Eigenmatrix = transpose(findEigenvectors(A, eigenvalues));
+ TinyMatrix<3> Diag = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ Diag(i, i) = eigenvalues[i];
+ }
+ TinyMatrix<3> B = Eigenmatrix * Diag * transpose(Eigenmatrix);
+ std::cout << "\n"
+ << B << ", " << A << " Diff "
+ << " A - B " << A - B << "\n";
+ exit(0);
+ 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;
- double vp = 2;
- TinyMatrix<3> A{1, -1, 0, -1, 1, 0, 0, 0, 3};
- std::cout << A << "\n";
- const TinyVector<3> eigenvalues = _findEigenValues(Sd);
- std::cout << eigenvalues << "\n";
- TinyMatrix<3> Try = findEigenvectors(A, vp);
- std::cout << Try << "\n";
- exit(0);
+ TinyMatrix<3> A{3, 2, 4, 2, 0, 2, 4, 2, 3};
+ // TinyMatrix<3> A{-1, 1, 0, 1, -1, 0, 0, 0, 3};
+ auto [eigenvalues, eigenmatrix] = _testMatrix(A);
+
double chi = 0;
const double normS = frobeniusNorm(Sd);
if ((normS - std::sqrt(2. / 3) * yield) > 0.) {
@@ -320,14 +428,15 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
Polynomial<3> P(d, c, b, 1.);
std::cout << P << "\n";
Polynomial<2> Q(c, b, 1.);
- if (fabs(d) > 1.e-10) {
+ if (fabs(d) > 1e-10) {
eigenvalues[0] = _findFirstRoot(P);
- Polynomial<1> Q1(eigenvalues[0], 1);
+ Polynomial<1> Q1(-eigenvalues[0], 1);
Polynomial<3> S;
Polynomial<3> R;
divide(P, Q1, S, R);
- std::cout << "S"
- << "\n";
+ 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];
@@ -335,6 +444,19 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
TinyVector<2> last_eigenvalues = _findLastRoots(Q);
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;
}
@@ -342,16 +464,22 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
_findFirstRoot(Polynomial<3> P) const
{
double 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;
- while (residu > 1.e-14) {
- guess -= P(guess) / Q(guess);
+ double residu = 0.;
+ do {
+ // guess -= P(guess) / Q(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;
- }
- std::cout << " nb Newton iterations " << iteration;
+ } while (fabs(residu) > 1.e-16);
+ std::cout << " nb Newton iterations " << iteration << "\n";
return guess;
}
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 0cd3da4a07272caff13228eca562e826e1612768..323d2d6a54058359f1d004ef2f0cf60f7f57f2af 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -35,23 +35,23 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
CRSMatrix A{S.getCRSMatrix()};
- Array<int> non_zeros2(3);
- non_zeros2.fill(3);
- CRSMatrixDescriptor<double> S2{3, 3, non_zeros};
+ // 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(0, 0) = 1;
+ // S2(0, 1) = 0;
+ // S2(0, 2) = 0;
- S2(1, 0) = 0;
- S2(1, 1) = 0;
- S2(1, 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;
+ // S2(2, 0) = 0;
+ // S2(2, 1) = 0;
+ // S2(2, 2) = 1;
- CRSMatrix B{S2.getCRSMatrix()};
+ // CRSMatrix B{S2.getCRSMatrix()};
SECTION("eigenvalues")
{
@@ -150,7 +150,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
{
SmallMatrix<double> expA{};
SmallMatrix<double> expS{3};
- SmallMatrix<double> expB{};
+ // SmallMatrix<double> expB{};
// SmallMatrix<double> expS2{3};
expS(0, 0) = 1325.074593930307812;
@@ -174,11 +174,11 @@ 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));
- }
- }
+ // 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),