From 9f15f3e31442c7956e4c52808649503424ff036a Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Thu, 25 May 2023 18:42:33 +0200
Subject: [PATCH] Add an epsilon to avoid collapsing stress tensor issue
---
src/scheme/HypoelasticSolver.cpp | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp
index f41632e76..c4076e353 100644
--- a/src/scheme/HypoelasticSolver.cpp
+++ b/src/scheme/HypoelasticSolver.cpp
@@ -491,10 +491,11 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
+ const double eps = 1.e-16;
const Rdxd I = identity;
const double norms = l2Norm(sigmad[j]);
- const double yieldf = norms - std::sqrt(2 / 3) * yieldy[j];
- const Rdxd N = 1. / norms * sigmad[j];
+ 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;
--
GitLab