From 6ee4681a48757cea73e8e0d1fd8a7ba1e2d89b8a Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Thu, 31 Aug 2023 10:43:29 +0200
Subject: [PATCH] bugs fix
---
src/scheme/HypoelasticSolver.cpp | 15 ++++++++++++---
1 file changed, 12 insertions(+), 3 deletions(-)
diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp
index c4076e353..6036a30f6 100644
--- a/src/scheme/HypoelasticSolver.cpp
+++ b/src/scheme/HypoelasticSolver.cpp
@@ -476,7 +476,12 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_sigmad = copy(sigmad.cellValues());
- auto chi = [&](const double a) -> double { return (a < 0 ? 0 : 1); };
+ auto chi = [&](const double a, const double b) -> double {
+ if ((a >= 0) and (b > 0))
+ return 1;
+ else
+ return 0;
+ };
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
@@ -497,8 +502,9 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
const double yieldf = norms - std::sqrt(2. / 3.) * yieldy[j];
const Rdxd N = 1. / (norms + eps) * sigmad[j];
const Rdxd D = 0.5 * (gradv + transpose(gradv));
- const Rdxd Dd = D - 1 / 3 * trace(D) * I;
- const Rdxd Dp = chi(yieldf) * scalarProduct(N, D) * N;
+ const Rdxd Dd = D - 1. / Dimension * trace(D) * I;
+ const double scalarND = scalarProduct(N, D);
+ const Rdxd Dp = chi(yieldf, scalarND) * scalarND * N;
const Rdxd gradva = 0.5 * (gradv - transpose(gradv));
const Rdxd jaumann = sigmad[j] * gradva - gradva * sigmad[j];
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
@@ -506,6 +512,9 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
new_sigmad[j] += dt_over_Vj * (2 * mu[j] * (Dd - Dp) - jaumann);
+ if (trace(new_sigmad[j]) >= 1.e-6) {
+ std::cout << j << " trace " << trace(new_sigmad[j]) << "\n";
+ }
});
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
--
GitLab