[cig-commits] r17087 - in short/3D/PyLith/branches/v1.5-stable: libsrc/materials unittests/libtests/materials/data

[willic3 at geodynamics.org]

Author: willic3
Date: 2010-08-12 20:49:12 -0700 (Thu, 12 Aug 2010)
New Revision: 17087

Modified:
   short/3D/PyLith/branches/v1.5-stable/libsrc/materials/DruckerPrager3D.cc
   short/3D/PyLith/branches/v1.5-stable/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
Log:
Added Drucker-Prager fixes to v1.5-stable that should have been done
originally (they were already in trunk).



Modified: short/3D/PyLith/branches/v1.5-stable/libsrc/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/libsrc/materials/DruckerPrager3D.cc	2010-08-12 22:31:38 UTC (rev 17086)
+++ short/3D/PyLith/branches/v1.5-stable/libsrc/materials/DruckerPrager3D.cc	2010-08-13 03:49:12 UTC (rev 17087)
@@ -540,10 +540,11 @@
 				      strainPPTpdt[5]/ae + devStressInitial[5]};
     const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
     const double yieldFunction = 3.0* alphaYield * trialMeanStress +
-      pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
-							  trialDevStress) -
+      sqrt(0.5 *
+	   pylith::materials::ElasticMaterial::scalarProduct3D( trialDevStress,
+								trialDevStress)) -
       beta;
-    PetscLogFlops(74);
+    PetscLogFlops(76);
 
     // If yield function is greater than zero, compute elastoplastic stress.
     if (yieldFunction >= 0.0) {
@@ -802,9 +803,11 @@
 				    strainPPTpdt[5]/ae + devStressInitial[5]};
   const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
   const double yieldFunction = 3.0* alphaYield * trialMeanStress +
-    pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
-							trialDevStress) - beta;
-  PetscLogFlops(74);
+    sqrt(0.5 *
+	 pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
+							     trialDevStress)) -
+    beta;
+  PetscLogFlops(76);
   
   // If yield function is greater than zero, compute elastoplastic stress and
   // corresponding tangent matrix.
@@ -1059,9 +1062,11 @@
 				    strainPPTpdt[5]/ae + devStressInitial[5]};
   const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
   const double yieldFunction = 3.0* alphaYield * trialMeanStress +
-    pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
-							trialDevStress) - beta;
-  PetscLogFlops(74);
+    sqrt(0.5 *
+	 pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
+							     trialDevStress)) -
+    beta;
+  PetscLogFlops(76);
 
   // If yield function is greater than zero, compute plastic strains.
   // Otherwise, plastic strains remain the same.

Modified: short/3D/PyLith/branches/v1.5-stable/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py	2010-08-12 22:31:38 UTC (rev 17086)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py	2010-08-13 03:49:12 UTC (rev 17087)
@@ -93,6 +93,8 @@
     # Second case has different values for friction angle and dilatation angle.
     frictionAngleB = math.radians(25.0)
     dilatationAngleB = math.radians(25.0)
+    # frictionAngleB = 0.0
+    # dilatationAngleB = 0.0
     cohesionB = 1.0e4
     strainB = [4.1e-4, 4.2e-4, 4.3e-4, 4.4e-4, 4.5e-4, 4.6e-4]
     initialStressB = [5.1e4, 5.2e4, 5.3e4, 5.4e4, 5.5e4, 5.6e4]
@@ -256,7 +258,9 @@
     trialDevStress = strainPPTpdt/ae + devStressInitial
     trialMeanStress = meanStrainPPTpdt/am + meanStressInitial
     yieldFunction = 3.0 * alphaYieldV * trialMeanStress + \
-                    self._scalarProduct(trialDevStress,  trialDevStress) - betaV
+                    math.sqrt(0.5 * self._scalarProduct(trialDevStress,
+                                                        trialDevStress)) - \
+                                                        betaV
 
     # If yield function is greater than zero, compute elastoplastic stress.
     if (yieldFunction >= 0.0):