[cig-commits] [commit] willic3/fix-plasticity: Modified Drucker-Prager for the case where there is no plastic yield and tensile yield is not allowed. (5265cab)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/pylith

On branch  : willic3/fix-plasticity
Link       : https://github.com/geodynamics/pylith/compare/5a98e61eb36798668e6db4cd9329b6f3430a112b...5265cabb42a53c87bf4721228408b095415f7e13

>---------------------------------------------------------------

commit 5265cabb42a53c87bf4721228408b095415f7e13
Author: Charles Williams <C.Williams at gns.cri.nz>
Date:   Tue Nov 18 11:51:48 2014 +1300

    Modified Drucker-Prager for the case where there is no plastic yield and
    tensile yield is not allowed.


>---------------------------------------------------------------

5265cabb42a53c87bf4721228408b095415f7e13
 libsrc/pylith/materials/DruckerPrager3D.cc         | 50 ++++++++++------------
 .../pylith/materials/DruckerPragerPlaneStrain.cc   | 50 ++++++++++------------
 2 files changed, 44 insertions(+), 56 deletions(-)

diff --git a/libsrc/pylith/materials/DruckerPrager3D.cc b/libsrc/pylith/materials/DruckerPrager3D.cc
index 5dce6d0..746a423 100644
--- a/libsrc/pylith/materials/DruckerPrager3D.cc
+++ b/libsrc/pylith/materials/DruckerPrager3D.cc
@@ -71,7 +71,7 @@ namespace pylith {
 	"vp" ,
 	"friction-angle",
 	"cohesion",
-	"dilatation-angle",
+	"dilatation-angle"
       };
 
       /// Number of state variables.
@@ -1254,25 +1254,6 @@ pylith::materials::DruckerPrager3D::_updateStateVarsElastoplastic(
   const PylithScalar yieldFunction =
     3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
 
-  if (!_allowTensileYield) {
-    const PylithScalar feasibleFunction =
-      3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
-    // If mean stress is too negative to project back to the yield surface,
-    // throw an exception.
-    if (feasibleFunction < 0.0) {
-      std::ostringstream msg;
-      msg << "Infeasible stress state - cannot project back to yield surface.\n"
-	  << "  alphaYield:       " << alphaYield << "\n"
-	  << "  alphaFlow:        " << alphaFlow << "\n"
-	  << "  beta:             " << beta << "\n"
-	  << "  trialMeanStress:  " << trialMeanStress << "\n"
-	  << "  stressInvar2:     " << stressInvar2 << "\n"
-	  << "  yieldFunction:    " << yieldFunction << "\n"
-	  << "  feasibleFunction: " << feasibleFunction << "\n";
-      throw std::runtime_error(msg.str());
-    } // if
-  } // if
-
 #if 0 // DEBUGGING
   std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
   std::cout << "  alphaYield:       " << alphaYield << std::endl;
@@ -1299,16 +1280,29 @@ pylith::materials::DruckerPrager3D::_updateStateVarsElastoplastic(
     const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
 
     PylithScalar plasticMult = 0.0;
-    if (_allowTensileYield) {
-      plasticMult = std::min(PylithScalar(sqrt(2.0)*d), 
-			     plasticFac * (meanStrainFac *
-					   (meanStrainPPTpdt/am +
-					    meanStressInitial) +
-					   dFac * d - beta));
-    } else {
-      plasticMult = plasticFac *
+    const PylithScalar plasticMultNormal = plasticFac *
 	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
 	 dFac * d - beta);
+    const PylithScalar plasticMultTensile = sqrt(2.0) * d;
+    
+    if (_allowTensileYield) {
+      plasticMult = std::min(plasticMultNormal, plasticMultTensile);
+    } else if (plasticMultNormal <= plasticMultTensile) {
+      plasticMult = plasticMultNormal;
+    } else {
+      std::ostringstream msg;
+      const PylithScalar stressInvar2Comp = 0.5 *
+	(sqrt(2.0) * d - plasticMultNormal)/ae;
+      msg << "Infeasible stress state. Cannot project back to yield surface.\n"
+	  << "  alphaYield:         " << alphaYield << "\n"
+	  << "  alphaFlow:          " << alphaFlow << "\n"
+	  << "  beta:               " << beta << "\n"
+	  << "  d:                  " << d << "\n"
+	  << "  plasticMultNormal:  " << plasticMultNormal << "\n"
+	  << "  plasticMultTensile: " << plasticMultTensile << "\n"
+	  << "  yieldFunction:      " << yieldFunction << "\n"
+	  << "  stressInvar2Comp:   " << stressInvar2Comp << "\n";
+      throw std::runtime_error(msg.str());
     } // if/else
 
     const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
diff --git a/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc b/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
index 81f4363..df46f26 100644
--- a/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
+++ b/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
@@ -1222,25 +1222,6 @@ pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastoplastic(
   const PylithScalar yieldFunction =
     3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
 
-  if (!_allowTensileYield) {
-    const PylithScalar feasibleFunction =
-      3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
-    // If mean stress is too negative to project back to the yield surface,
-    // throw an exception.
-    if (feasibleFunction < 0.0) {
-      std::ostringstream msg;
-      msg << "Infeasible stress state - cannot project back to yield surface.\n"
-	  << "  alphaYield:       " << alphaYield << "\n"
-	  << "  alphaFlow:        " << alphaFlow << "\n"
-	  << "  beta:             " << beta << "\n"
-	  << "  trialMeanStress:  " << trialMeanStress << "\n"
-	  << "  stressInvar2:     " << stressInvar2 << "\n"
-	  << "  yieldFunction:    " << yieldFunction << "\n"
-	  << "  feasibleFunction: " << feasibleFunction << "\n";
-      throw std::runtime_error(msg.str());
-    } // if
-  } // if
-
 #if 0 // DEBUGGING
   std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
   std::cout << "  alphaYield:       " << alphaYield << std::endl;
@@ -1268,16 +1249,28 @@ pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastoplastic(
     const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
 
     PylithScalar plasticMult =  0.0;
+    const PylithScalar plasticMultNormal = plasticFac *
+      (meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
+       dFac * d - beta);
+  const PylithScalar plasticMultTensile = sqrt(2.0) * d;
     if (_allowTensileYield) {
-      plasticMult = std::min(PylithScalar(sqrt(2.0) * d),
-			     plasticFac * (meanStrainFac *
-					   (meanStrainPPTpdt/am +
-					    meanStressInitial) +
-					   dFac * d - beta));
+      plasticMult = std::min(plasticMultNormal, plasticMultTensile);
+    } else if (plasticMultNormal <= plasticMultTensile) {
+      plasticMult = plasticMultNormal;
     } else {
-      plasticMult = plasticFac *
-	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
-	 dFac * d - beta);
+      std::ostringstream msg;
+      const PylithScalar stressInvar2Comp = 0.5 *
+	(sqrt(2.0) * d - plasticMultNormal)/ae;
+      msg << "Infeasible stress state. Cannot project back to yield surface.\n"
+	  << "  alphaYield:         " << alphaYield << "\n"
+	  << "  alphaFlow:          " << alphaFlow << "\n"
+	  << "  beta:               " << beta << "\n"
+	  << "  d:                  " << d << "\n"
+	  << "  plasticMultNormal:  " << plasticMultNormal << "\n"
+	  << "  plasticMultTensile: " << plasticMultTensile << "\n"
+	  << "  yieldFunction:      " << yieldFunction << "\n"
+	  << "  stressInvar2Comp:   " << stressInvar2Comp << "\n";
+      throw std::runtime_error(msg.str());
     } // if/else
 
     const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
@@ -1292,7 +1285,8 @@ pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastoplastic(
       } // for
     } else {
       for (int iComp=0; iComp < tensorSizePS; ++iComp) {
-	stateVars[s_plasticStrain+iComp] += diag[iComp] * deltaMeanPlasticStrain;
+	stateVars[s_plasticStrain+iComp] +=
+	  diag[iComp] * deltaMeanPlasticStrain;
       } // for
     } // if/else