Skip to content
Snippets Groups Projects
Commit d15023ee authored by Emmanuel Labourasse's avatar Emmanuel Labourasse
Browse files

fix a bug in the plasticity management

parent e380b4c0
No related branches found
No related tags found
No related merge requests found
...@@ -368,6 +368,10 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige ...@@ -368,6 +368,10 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige
TinyMatrix<3> Try = findEigenvector(A, eigenvalues[count], space_size[count], epsilon); TinyMatrix<3> Try = findEigenvector(A, eigenvalues[count], space_size[count], epsilon);
// std::cout << eigenvalues[count] << " Frobenius norm try: " << frobeniusNorm(Try) << "\n"; // std::cout << eigenvalues[count] << " Frobenius norm try: " << frobeniusNorm(Try) << "\n";
Assert(frobeniusNorm(Try) > 1e-3, " Error : no eigenvector corresponds to eigenvalue "); Assert(frobeniusNorm(Try) > 1e-3, " Error : no eigenvector corresponds to eigenvalue ");
// if (frobeniusNorm(Try) < 1e-3) {
// std::cout << " Error : no eigenvector corresponds to eigenvalue \n";
// exit(0);
// }
TinyVector<3> eigenvector = zero; TinyVector<3> eigenvector = zero;
int count2 = 0; int count2 = 0;
...@@ -385,6 +389,10 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige ...@@ -385,6 +389,10 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige
} }
} }
Assert(count2 == space_size[count - count2], " eigenvector space size is not what was expected "); Assert(count2 == space_size[count - count2], " eigenvector space size is not what was expected ");
// if (count2 != space_size[count - count2]) {
// std::cout << " eigenvector space size is not what was expected \n";
// exit(0);
// }
it++; it++;
} }
// std ::cout << " eigenMatrix " << eigenMatrix << "\n"; // std ::cout << " eigenMatrix " << eigenMatrix << "\n";
...@@ -479,7 +487,7 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const ...@@ -479,7 +487,7 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
// double delta = pow(2.*P.coefficients()[2],2) - 4*3.*P.coefficients()[3]*P.coefficients()[1]; // double delta = pow(2.*P.coefficients()[2],2) - 4*3.*P.coefficients()[3]*P.coefficients()[1];
Polynomial<2> Q = derivative(P); Polynomial<2> Q = derivative(P);
Polynomial<1> R = derivative(Q); Polynomial<1> R = derivative(Q);
size_t iteration = 0; // size_t iteration = 0;
double residu = 0.; double residu = 0.;
do { do {
// guess -= P(guess) / Q(guess); // guess -= P(guess) / Q(guess);
...@@ -489,7 +497,7 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const ...@@ -489,7 +497,7 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
// std::cout << "g(guess) = " << Q(guess) << "\n"; // std::cout << "g(guess) = " << Q(guess) << "\n";
residu = P(guess); residu = P(guess);
// std::cout << "residu = " << residu << "\n"; // std::cout << "residu = " << residu << "\n";
iteration += 1; // iteration += 1;
} while ((fabs(residu) > 1.e-16) or (fabs(old_guess - guess) > 1.e-10)); } while ((fabs(residu) > 1.e-16) or (fabs(old_guess - guess) > 1.e-10));
// std::cout << " nb Newton iterations " << iteration << "\n"; // std::cout << " nb Newton iterations " << iteration << "\n";
return guess; return guess;
...@@ -500,6 +508,10 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const ...@@ -500,6 +508,10 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const
{ {
TinyVector<2> roots(0., 0.); TinyVector<2> roots(0., 0.);
double delta = pow(P.coefficients()[1], 2) - 4 * P.coefficients()[2] * P.coefficients()[0]; double delta = pow(P.coefficients()[1], 2) - 4 * P.coefficients()[2] * P.coefficients()[0];
// if (fabs(delta) > epsilon) {
// std::cout << "Find roots is only for symmetric matrix \n";
// exit(0);
// }
Assert(delta > -epsilon, "Find roots is only for symmetric matrix"); Assert(delta > -epsilon, "Find roots is only for symmetric matrix");
if (fabs(delta) < epsilon) { if (fabs(delta) < epsilon) {
roots[0] = -P.coefficients()[1] / (2. * P.coefficients()[2]); roots[0] = -P.coefficients()[1] / (2. * P.coefficients()[2]);
...@@ -512,15 +524,17 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const ...@@ -512,15 +524,17 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const
return roots; return roots;
} }
void TinyMatrix<3>
EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA) EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A)
{ {
TinyMatrix<3> expA;
auto [eigenvalues, eigenvectors] = findEigen(A); auto [eigenvalues, eigenvectors] = findEigen(A);
TinyMatrix<3> D = zero; TinyMatrix<3> D = zero;
for (size_t i = 0; i < 3; ++i) { for (size_t i = 0; i < 3; ++i) {
D(i, i) = std::exp(eigenvalues[i]); D(i, i) = std::exp(eigenvalues[i]);
} }
expA = eigenvectors * D * transpose(eigenvectors); expA = eigenvectors * D * transpose(eigenvectors);
return expA;
} }
EigenvalueSolver::EigenvalueSolver() {} EigenvalueSolver::EigenvalueSolver() {}
...@@ -100,7 +100,7 @@ class EigenvalueSolver ...@@ -100,7 +100,7 @@ class EigenvalueSolver
TinyVector<2> _findLastRoots(Polynomial<2> P, const double epsilon) const; TinyVector<2> _findLastRoots(Polynomial<2> P, const double epsilon) const;
void computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA); TinyMatrix<3> computeExpForTinyMatrix(const TinyMatrix<3>& A);
EigenvalueSolver(); EigenvalueSolver();
~EigenvalueSolver() = default; ~EigenvalueSolver() = default;
......
...@@ -170,46 +170,13 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS ...@@ -170,46 +170,13 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
return vector; return vector;
} }
R3x3
matrix_exp(const R3x3 A) const
{
double normA = frobeniusNorm(A);
R3x3 expA = identity;
if (normA <= 1.e-12) {
return expA;
}
SmallMatrix<double> expB{3};
SmallMatrix<double> B(3, 3);
// B.fill(0.);
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
B(i, i) = A(i, j);
}
}
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;
}
double double
_compute_chi(const R3x3 S, const double yield) const _compute_chi(const R3x3 S, const double yield) const
{ {
const R3x3 Id = identity; const R3x3 Id = identity;
const R3x3 Sd = S - 1. / 3 * trace(S) * Id; const R3x3 Sd = S - 1. / 3 * trace(S) * Id;
// auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestB); // auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestB);
R3x3 B; // exit(0);
EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, B);
std::cout << B << "\n";
exit(0);
double chi = 0; double chi = 0;
const double normS = frobeniusNorm(Sd); const double normS = frobeniusNorm(Sd);
if ((normS - std::sqrt(2. / 3) * yield) > 0.) { if ((normS - std::sqrt(2. / 3) * yield) > 0.) {
...@@ -612,7 +579,9 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS ...@@ -612,7 +579,9 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
CellValue<double> new_yield = copy(yield.cellValues()); CellValue<double> new_yield = copy(yield.cellValues());
CellValue<R3x3> sigma_s = copy(sigma.cellValues()); CellValue<R3x3> sigma_s = copy(sigma.cellValues());
CellValue<double> chi{mesh.connectivity()}; CellValue<double> chi{mesh.connectivity()};
CellValue<double> rationorm{mesh.connectivity()};
chi.fill(0.); chi.fill(0.);
rationorm.fill(0.);
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j]; const auto& cell_nodes = cell_to_node_matrix[j];
...@@ -626,18 +595,26 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS ...@@ -626,18 +595,26 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
momentum_fluxes += Fjr(j, R); momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]); energy_fluxes += dot(Fjr(j, R), ur[r]);
} }
R3x3 I = identity;
R3x3 gradv3d = to3D(gradv); R3x3 gradv3d = to3D(gradv);
R3x3 sigmas = sigma[j] - 1. / 3 * trace(sigma[j]) * I;
const R3x3 elastic_fluxes = gradv3d * Fe[j]; const R3x3 elastic_fluxes = gradv3d * Fe[j];
const double dt_over_Mj = dt / (rho[j] * Vj[j]); const double dt_over_Mj = dt / (rho[j] * Vj[j]);
const double dt_over_Vj = dt / Vj[j]; const double dt_over_Vj = dt / Vj[j];
new_u[j] += dt_over_Mj * momentum_fluxes; new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes; new_E[j] += dt_over_Mj * energy_fluxes;
new_Fe[j] += dt_over_Vj * elastic_fluxes; new_Fe[j] += dt_over_Vj * elastic_fluxes;
const double fenorm0 = frobeniusNorm(new_Fe[j]);
chi[j] = _compute_chi(sigma[j], yield[j]); chi[j] = _compute_chi(sigma[j], yield[j]);
const R3x3 M = dt * chi[j] * zeta[j] * new_Fe[j]; const R3x3 M = -dt * chi[j] * zeta[j] * sigmas;
new_Fe[j] = new_Fe[j] * matrix_exp(M);
new_Fe[j] = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
const double fenorm1 = frobeniusNorm(new_Fe[j]);
rationorm[j] = fenorm1 / fenorm0;
}); });
std::cout << "sum_chi " << sum(chi) << "\n"; std::cout << "sum_chi " << sum(chi) << " min norm ratio " << min(rationorm) << " max norm ratio " << max(rationorm)
<< "\n";
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj(); CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for( parallel_for(
......
...@@ -190,7 +190,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]") ...@@ -190,7 +190,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
REQUIRE((B(i, j) - TestA2(i, j)) / frobeniusNorm(TestA2) == Catch::Approx(0).margin(1E-8)); REQUIRE((B(i, j) - TestA2(i, j)) / frobeniusNorm(TestA2) == Catch::Approx(0).margin(1E-8));
} }
} }
EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, expA2); expA2 = EigenvalueSolver{}.computeExpForTinyMatrix(TestA2);
for (size_t i = 0; i < 3; ++i) { for (size_t i = 0; i < 3; ++i) {
Diag(i, i) = std::exp(eigenvalues2[i]); Diag(i, i) = std::exp(eigenvalues2[i]);
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment