diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp
index f41632e761082a55ac9cf83f4c313af005366010..c4076e353bb8ed3d1e8e9346e3879592014c57f0 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;