bugs fix

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();