From 3609391a8339565c1b49431386f6be725a4d7013 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Thu, 20 Feb 2025 17:41:32 +0100
Subject: [PATCH] ?
---
src/algebra/EigenvalueSolver.cpp | 17 ++++++++++-------
src/mesh/MeshFlatNodeBoundary.cpp | 12 ++++++------
src/scheme/HyperplasticSolverO2.cpp | 29 +++++++++++++++++++++++++++++
3 files changed, 45 insertions(+), 13 deletions(-)
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 378cebb94..4d915dc10 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -488,18 +488,21 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
Polynomial<2> Q = derivative(P);
Polynomial<1> R = derivative(Q);
// size_t iteration = 0;
- double residu = 0.;
+ double residu = 0.;
+ size_t iteration = 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";
+ // std::cout << "guess = " << guess << " old_guess " << old_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";
+ // std::cout << "residu = " << residu << " delta " << fabs(old_guess - guess) << "\n";
+ iteration += 1;
+ } while (((fabs(residu) > 1.e-16) or (fabs(old_guess - guess) > 1.e-10)) and (iteration < 10000));
+ if (iteration == 100) {
+ std::cout << " nb Newton iterations reached, residu " << residu << "\n";
+ }
return guess;
}
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index 647298e97..5621a6230 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -24,12 +24,12 @@ MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType::
}
});
- if (parallel::allReduceOr(is_bad)) {
- std::ostringstream ost;
- ost << "invalid boundary \"" << rang::fgB::yellow << this->m_ref_node_list.refId() << rang::style::reset
- << "\": boundary is not flat!";
- throw NormalError(ost.str());
- }
+ // if (parallel::allReduceOr(is_bad)) {
+ // std::ostringstream ost;
+ // ost << "invalid boundary \"" << rang::fgB::yellow << this->m_ref_node_list.refId() << rang::style::reset
+ // << "\": boundary is not flat!";
+ // throw NormalError(ost.str());
+ // }
}
template <>
diff --git a/src/scheme/HyperplasticSolverO2.cpp b/src/scheme/HyperplasticSolverO2.cpp
index 41543fca0..c4efdd128 100644
--- a/src/scheme/HyperplasticSolverO2.cpp
+++ b/src/scheme/HyperplasticSolverO2.cpp
@@ -769,6 +769,35 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
CellValue<double> rationorm{mesh.connectivity()};
chi.fill(0.);
rationorm.fill(0.);
+ // for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+ // const auto& cell_nodes = cell_to_node_matrix[j];
+
+ // Rd momentum_fluxes = zero;
+ // double energy_fluxes = 0;
+ // Rdxd gradv = zero;
+ // for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ // const NodeId r = cell_nodes[R];
+ // gradv += tensorProduct(ur[r], Cjr(j, R));
+ // momentum_fluxes += Fjr(j, R);
+ // energy_fluxes += dot(Fjr(j, R), ur[r]);
+ // }
+ // R3x3 I = identity;
+ // R3x3 gradv3d = to3D(gradv);
+ // R3x3 sigmas = sigma[j] - 1. / 3 * trace(sigma[j]) * I;
+ // const R3x3 elastic_fluxes = gradv3d * Fe[j];
+ // const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ // const double dt_over_Vj = dt / Vj[j];
+ // new_u[j] += dt_over_Mj * momentum_fluxes;
+ // new_E[j] += dt_over_Mj * energy_fluxes;
+ // new_Fe[j] += dt_over_Vj * elastic_fluxes;
+ // const double fenorm0 = frobeniusNorm(new_Fe[j]);
+ // chi[j] = _compute_chi(sigma[j], yield[j]);
+ // const R3x3 M = -dt * chi[j] * zeta[j] * sigmas;
+
+ // new_Fe[j] = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
+ // const double fenorm1 = frobeniusNorm(new_Fe[j]);
+ // rationorm[j] = fenorm1 / fenorm0;
+ // }
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
--
GitLab