clean implementation of criterion based on the interface force

NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp

@@ -410,20 +410,34 @@ void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, cons
         // If Coalescence has not yet occured at least once
         // Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface
 
-        // Get normal traction on the interface from tau stress
-        double sigmax = 0.;
+        // Get effective traction stress on the interface from tau stress
+        static SVector3 tractionForce; // traction force vector
+        double tauEff = 0.; double tauN = 0.; double tauTanSq = 0.; // tranction force components
+        double beta = 1.; // ratio => should be moved into law creation as a parameter
         static SVector3 Norm;
-        static SVector3 tractionForce;
         Norm = q1->getConstRefToReferenceOutwardNormal();
         if (Norm.norm() > 0){ // is it usefulll ???
           Norm.normalize();
+          // Get normal part
           for (int i=0; i<3; i++){
             tractionForce(i) = 0.;
             for (int j=0; j<3; j++){
               tractionForce(i) += kCor(i,j)*Norm(j);
             }
+           tauN += tractionForce(i)*Norm(i);
+          }
+          // Deduce tangent part
+          for (int i=0; i<3; i++){
+            double temp = tractionForce(i)-tauN*Norm(i);
+            tauTanSq += temp*temp;
+          }
+          // Compute effective stress
+          if (tauN > 0.){
+            tauEff = sqrt(tauN*tauN + beta*tauTanSq);
+          }
+          else{
+            tauEff = sqrt(beta*tauTanSq);
           }
-          sigmax = tractionForce.norm();
         }
         else{
           Msg::Fatal("mlawNonLocalPorousCoupledLaw::updatePlasticState: zero normal");
@@ -447,7 +461,7 @@ void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, cons
           q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
         }
 
-        if ((sigmax - Cft*R > _tol) and (yieldThomason > _tol)){
+        if ((tauEff - Cft*R > _tol) and (yieldThomason > _tol)){
           // If coalescence occurs, update Cft offset value and others variables
 
           // Update onset coalescence variables