[cig-commits] [commit] baagaard/add-release-2.1.0, baagaard/dynrup-new-lagrange, baagaard/feature-output-station-names, baagaard/feature-progress-monitor, baagaard/fix-custom-faultpc, baagaard/fix-faults-intersect, knepley/feature-petsc-fe, knepley/upgrade-petsc-master, maint, master, willic3/fix-plasticity: Modified manual, code, and unit tests for error in Drucker-Prager elastoplasticity. We were missing the initial pressure term. I haven't tested yet. (b59a515)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Wed Nov 5 15:46:22 PST 2014


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

On branches: baagaard/add-release-2.1.0,baagaard/dynrup-new-lagrange,baagaard/feature-output-station-names,baagaard/feature-progress-monitor,baagaard/fix-custom-faultpc,baagaard/fix-faults-intersect,knepley/feature-petsc-fe,knepley/upgrade-petsc-master,maint,master,willic3/fix-plasticity
Link       : https://github.com/geodynamics/pylith/compare/f33c75b19fd60eedb2a3405db76a1fee333bb1d7...5b6d812b1612809fea3bd331c4e5af98c25a536a

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

commit b59a515bf037254a708f031fa2772ca388b4ea65
Author: Charles Williams <C.Williams at gns.cri.nz>
Date:   Mon Oct 6 17:01:55 2014 +1300

    Modified manual, code, and unit tests for error in Drucker-Prager
    elastoplasticity. We were missing the initial pressure term. I haven't
    tested yet.


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

b59a515bf037254a708f031fa2772ca388b4ea65
 doc/userguide/materials/materials.lyx              | 93 +++++++++++----------
 libsrc/pylith/materials/DruckerPrager3D.cc         | 24 ++++--
 .../pylith/materials/DruckerPragerPlaneStrain.cc   | 24 ++++--
 .../materials/data/DruckerPrager3DTimeDep.py       |  2 +-
 .../materials/data/DruckerPrager3DTimeDepData.cc   | 96 +++++++++++-----------
 .../data/DruckerPragerPlaneStrainTimeDep.py        |  2 +-
 .../data/DruckerPragerPlaneStrainTimeDepData.cc    | 32 ++++----
 7 files changed, 149 insertions(+), 124 deletions(-)

diff --git a/doc/userguide/materials/materials.lyx b/doc/userguide/materials/materials.lyx
index 2cedeaa..a0c8bd9 100644
--- a/doc/userguide/materials/materials.lyx
+++ b/doc/userguide/materials/materials.lyx
@@ -1,5 +1,5 @@
-#LyX 2.0 created this file. For more info see http://www.lyx.org/
-\lyxformat 413
+#LyX 2.1 created this file. For more info see http://www.lyx.org/
+\lyxformat 474
 \begin_document
 \begin_header
 \textclass book
@@ -10,18 +10,18 @@
 \maintain_unincluded_children false
 \language english
 \language_package default
-\inputencoding latin1
+\inputencoding iso8859-1
 \fontencoding global
 \font_roman default
 \font_sans default
 \font_typewriter default
+\font_math auto
 \font_default_family default
 \use_non_tex_fonts false
 \font_sc false
 \font_osf false
 \font_sf_scale 100
 \font_tt_scale 100
-
 \graphics default
 \default_output_format default
 \output_sync 0
@@ -32,15 +32,24 @@
 \use_hyperref false
 \papersize default
 \use_geometry true
-\use_amsmath 1
-\use_esint 0
-\use_mhchem 1
-\use_mathdots 1
+\use_package amsmath 1
+\use_package amssymb 1
+\use_package cancel 1
+\use_package esint 0
+\use_package mathdots 1
+\use_package mathtools 1
+\use_package mhchem 1
+\use_package stackrel 1
+\use_package stmaryrd 1
+\use_package undertilde 1
 \cite_engine basic
+\cite_engine_type default
+\biblio_style plain
 \use_bibtopic false
 \use_indices false
 \paperorientation portrait
 \suppress_date false
+\justification true
 \use_refstyle 0
 \index Index
 \shortcut idx
@@ -502,7 +511,7 @@ status open
 
 \begin_layout Plain Layout
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -529,7 +538,7 @@ cell_data_fields
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="6" columns="4">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="center" valignment="top" width="1.5in">
 <column alignment="center" valignment="top" width="1.8in">
 <column alignment="center" valignment="top" width="1.5in">
@@ -850,7 +859,7 @@ status open
 
 \begin_layout Plain Layout
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -867,7 +876,7 @@ Order of components in tensor state-variables for material models.
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="5" columns="4">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="center" valignment="top" width="1.25in">
 <column alignment="center" valignment="top" width="0.5in">
 <column alignment="center" valignment="top" width="1.25in">
@@ -1445,7 +1454,7 @@ status open
 \begin_layout Plain Layout
 \noindent
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 Values in spatial database used as parameters in the elastic material constituti
@@ -1457,7 +1466,7 @@ ve models.
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="4" columns="2">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -1917,7 +1926,7 @@ status open
 \begin_layout Plain Layout
 \noindent
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -1934,7 +1943,7 @@ Available viscoelastic materials for PyLith.
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="7" columns="2">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="left" valignment="top" width="2.85in">
 <column alignment="center" valignment="top" width="2.47in">
 <row>
@@ -2111,7 +2120,7 @@ status open
 \end_inset
 
 
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -2284,7 +2293,7 @@ status open
 \begin_layout Plain Layout
 \noindent
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -2301,10 +2310,10 @@ Mathematical notation used in this section.
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="3" columns="3">
-<features tabularvalignment="middle">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
+<features rotate="0" tabularvalignment="middle">
+<column alignment="center" valignment="top">
+<column alignment="center" valignment="top">
+<column alignment="center" valignment="top">
 <row>
 <cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
 \begin_inset Text
@@ -2967,7 +2976,7 @@ status open
 \begin_layout Plain Layout
 \noindent
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -2985,7 +2994,7 @@ Values in spatial database used as parameters in the linear Maxwell viscoelastic
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="5" columns="2">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -3137,7 +3146,7 @@ status open
 \begin_layout Plain Layout
 \noindent
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -3155,7 +3164,7 @@ Values in spatial database used as parameters in the generalized linear
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="10" columns="2">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -4480,7 +4489,7 @@ status open
 \begin_layout Plain Layout
 \noindent
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -4498,7 +4507,7 @@ Values in spatial database used as parameters in the nonlinear power-law
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="7" columns="2">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -4943,7 +4952,7 @@ sideways false
 status open
 
 \begin_layout Plain Layout
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -4969,11 +4978,11 @@ fit_mohr_coulomb
 \align center
 \begin_inset Tabular
 <lyxtabular version="3" rows="4" columns="4">
-<features tabularvalignment="middle">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
+<features rotate="0" tabularvalignment="middle">
+<column alignment="center" valignment="top">
+<column alignment="center" valignment="top">
+<column alignment="center" valignment="top">
+<column alignment="center" valignment="top">
 <row>
 <cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
 \begin_inset Text
@@ -5346,7 +5355,7 @@ reference "eq:105"
 
 \end_inset
 
- and taking the scalar product of both sides:
+ and taking the scalar inner product of both sides:
 \begin_inset Formula 
 \begin{equation}
 \lambda=\sqrt{2}\,\phantom{}^{t+\Delta t}d-2a_{E}\sqrt{^{t+\Delta t}J_{2}^{\prime}}\:,\label{eq:110}
@@ -5403,7 +5412,7 @@ We then use the yield condition (
  to obtain:
 \begin_inset Formula 
 \begin{equation}
-\lambda=\frac{2a_{E}a_{m}\left(\frac{3\alpha_{f}}{a_{m}}\phantom{}^{t+\Delta t}\theta^{\prime}+\frac{^{t+\Delta t}d}{\sqrt{2}a_{E}}-\beta\bar{c}\right)}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\:.\label{eq:114}
+\lambda=\frac{2a_{E}a_{m}\left(\frac{3\alpha_{f}}{a_{m}}\phantom{}^{t+\Delta t}\theta^{\prime}+\frac{^{t+\Delta t}d}{\sqrt{2}a_{E}}-\beta+3\alpha_{f}P^{I}\right)}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\:.\label{eq:114}
 \end{equation}
 
 \end_inset
@@ -5689,7 +5698,7 @@ status open
 
 \begin_layout Plain Layout
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 \begin_inset CommandInset label
@@ -5707,9 +5716,9 @@ c model with perfect plasticity.
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="7" columns="2">
-<features tabularvalignment="middle">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
+<features rotate="0" tabularvalignment="middle">
+<column alignment="center" valignment="top">
+<column alignment="center" valignment="top">
 <row>
 <cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
 \begin_inset Text
@@ -6135,7 +6144,7 @@ status open
 \begin_layout Plain Layout
 \noindent
 \align center
-\begin_inset Caption
+\begin_inset Caption Standard
 
 \begin_layout Plain Layout
 Values in spatial database for initial state variables for 3D problems.
@@ -6156,7 +6165,7 @@ reference "tab:material-model-output"
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="3" columns="2">
-<features tabularvalignment="middle">
+<features rotate="0" tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
diff --git a/libsrc/pylith/materials/DruckerPrager3D.cc b/libsrc/pylith/materials/DruckerPrager3D.cc
index 8c96583..5dce6d0 100644
--- a/libsrc/pylith/materials/DruckerPrager3D.cc
+++ b/libsrc/pylith/materials/DruckerPrager3D.cc
@@ -651,12 +651,14 @@ pylith::materials::DruckerPrager3D::_calcStressElastoplastic(
 	plasticMult = 
 	  std::min(sqrt(2.0)*d,
 		   2.0 * ae * am *
-		   (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae)
+		   (3.0 * alphaYield * (meanStrainPPTpdt/am + meanStressInitial)
+		    + d/(sqrt(2.0) * ae)
 		    - beta)/ (6.0 * alphaYield * alphaFlow * ae + am));
       } else {
 	plasticMult = 
 	  2.0 * ae * am *
-	  (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+	  (3.0 * alphaYield * (meanStrainPPTpdt/am + meanStressInitial)
+	   + d/(sqrt(2.0) * ae) - beta)/
 	  (6.0 * alphaYield * alphaFlow * ae + am);
       } // if/else
 
@@ -968,7 +970,7 @@ pylith::materials::DruckerPrager3D::_calcElasticConstsElastoplastic(
 	   scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
     const PylithScalar plasticFac = 2.0 * ae * am/
       (6.0 * alphaYield * alphaFlow * ae + am);
-    const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
+    const PylithScalar meanStrainFac = 3.0 * alphaYield;
     const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
 
     PylithScalar plasticMult = 0.0;
@@ -976,13 +978,15 @@ pylith::materials::DruckerPrager3D::_calcElasticConstsElastoplastic(
     PylithScalar dFac2 = 0.0;
     if (_allowTensileYield) {
       const PylithScalar testMult = plasticFac *
-	(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
+	 dFac * d - beta);
       tensileYield = (sqrt(2.0)*d < testMult) ? true : false;
       plasticMult = (tensileYield) ? sqrt(2.0)*d : testMult;
       dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
     } else {
       plasticMult = plasticFac *
-	(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
+	 dFac * d - beta);
       dFac2 = 1.0/(sqrt(2.0) * d);
     } // if/else
 
@@ -1291,16 +1295,20 @@ pylith::materials::DruckerPrager3D::_updateStateVarsElastoplastic(
 	   scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
     const PylithScalar plasticFac = 2.0 * ae * am/
       (6.0 * alphaYield * alphaFlow * ae + am);
-    const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
+    const PylithScalar meanStrainFac = 3.0 * alphaYield;
     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 + dFac * d - beta));
+			     plasticFac * (meanStrainFac *
+					   (meanStrainPPTpdt/am +
+					    meanStressInitial) +
+					   dFac * d - beta));
     } else {
       plasticMult = plasticFac *
-	(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
+	 dFac * d - beta);
     } // if/else
 
     const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
diff --git a/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc b/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
index 5594730..81f4363 100644
--- a/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
+++ b/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
@@ -666,12 +666,14 @@ pylith::materials::DruckerPragerPlaneStrain::_calcStressElastoplastic(
 	plasticMult = 
 	  std::min(sqrt(2.0)*d,
 		   2.0 * ae * am *
-		   (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) 
+		   (3.0 * alphaYield * (meanStrainPPTpdt/am + meanStressInitial)
+		    + d/(sqrt(2.0) * ae) 
 		    - beta)/ (6.0 * alphaYield * alphaFlow * ae + am));
       } else {
 	plasticMult = 
 	  2.0 * ae * am *
-	  (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+	  (3.0 * alphaYield * (meanStrainPPTpdt/am + meanStressInitial) +
+	   d/(sqrt(2.0) * ae) - beta)/
 	  (6.0 * alphaYield * alphaFlow * ae + am);
       } // if/else
 
@@ -982,7 +984,7 @@ pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastoplastic(
 	   strainPPTpdtProd);
     const PylithScalar plasticFac = 2.0 * ae * am/
       (6.0 * alphaYield * alphaFlow * ae + am);
-    const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
+    const PylithScalar meanStrainFac = 3.0 * alphaYield;
     const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
 
     PylithScalar plasticMult = 0.0;
@@ -990,13 +992,15 @@ pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastoplastic(
     PylithScalar dFac2 = 0.0;
     if (_allowTensileYield) {
       const PylithScalar testMult = plasticFac *
-	(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
+	 dFac * d - beta);
       tensileYield = (sqrt(2.0) * d < testMult) ? true: false;
       plasticMult = tensileYield ? sqrt(2.0) * d : testMult;
       dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
     } else {
       plasticMult = plasticFac *
-	(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
+	 dFac * d - beta);
       dFac2 = 1.0/(sqrt(2.0) * d);
     } // if/else
 
@@ -1260,16 +1264,20 @@ pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastoplastic(
 	   strainPPTpdtProd);
     const PylithScalar plasticFac = 2.0 * ae * am/
       (6.0 * alphaYield * alphaFlow * ae + am);
-    const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
+    const PylithScalar meanStrainFac = 3.0 * alphaYield;
     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 + dFac * d - beta));
+			     plasticFac * (meanStrainFac *
+					   (meanStrainPPTpdt/am +
+					    meanStressInitial) +
+					   dFac * d - beta));
     } else {
       plasticMult = plasticFac *
-	(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+	(meanStrainFac * (meanStrainPPTpdt/am + meanStressInitial) +
+	 dFac * d - beta);
     } // if/else
 
     const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
diff --git a/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py b/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
index bad3394..fe7f3e4 100644
--- a/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
+++ b/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
@@ -278,7 +278,7 @@ class DruckerPrager3DTimeDep(ElasticMaterialApp):
                     strainPPTpdtProd)
       dFac = math.sqrt(2.0) * d
       testMult = 2.0 * ae * am * \
-                 (3.0 * alphaYieldV * meanStrainPPTpdt/am + \
+                 (3.0 * alphaYieldV * (meanStrainPPTpdt/am + meanStressInitial) + \
                   d/(math.sqrt(2.0) * ae) - betaV)/ \
                   (6.0 * alphaYieldV * alphaFlowV * ae + am)
       plasticMult = min(testMult, dFac)
diff --git a/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc b/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
index a3698f9..7bae235 100644
--- a/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
+++ b/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
@@ -191,12 +191,12 @@ const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_strain[] = {
 };
 
 const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_stress[] = {
- -1.63946744e+07,
- -6.17719414e+06,
- -5.86690637e+06,
- -2.02684646e+06,
- -2.02615846e+06,
- -2.02547046e+06,
+ -1.63921277e+07,
+ -6.18945621e+06,
+ -5.87961817e+06,
+ -2.02390884e+06,
+ -2.02322184e+06,
+ -2.02253484e+06,
   4.21392632e+06,
   4.21392632e+06,
   4.21392632e+06,
@@ -206,42 +206,42 @@ const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_stress[] = {
 };
 
 const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_elasticConsts[] = {
-  6.57004938e+10,
-  2.99833221e+10,
-  2.98388535e+10,
-  1.88737456e+09,
-  1.88673474e+09,
-  1.88609213e+09,
-  2.09815549e+10,
-  3.02171782e+10,
- -1.40253222e+09,
-  8.61894060e+09,
-  8.61601438e+09,
-  8.61308910e+09,
-  2.05637179e+10,
- -1.67589867e+09,
-  2.86086951e+10,
-  8.82336730e+09,
-  8.82037124e+09,
-  8.81737750e+09,
-  2.72937212e+09,
-  6.09515538e+09,
-  6.19736873e+09,
-  2.96246322e+10,
- -1.33489072e+09,
- -1.33443763e+09,
-  2.72844546e+09,
-  6.09308621e+09,
-  6.19526533e+09,
- -1.33489084e+09,
-  2.96255386e+10,
- -1.33398466e+09,
-  2.72751914e+09,
-  6.09101739e+09,
-  6.19316124e+09,
- -1.33443763e+09,
- -1.33398443e+09,
-  2.96264446e+10,
+  6.56928355e+10,
+  2.99876509e+10,
+  2.98421849e+10,
+  1.90042052e+09,
+  1.89977698e+09,
+  1.89913157e+09,
+  2.09858837e+10,
+  3.01923398e+10,
+ -1.38202077e+09,
+  8.61270912e+09,
+  8.60978523e+09,
+  8.60686414e+09,
+  2.05670479e+10,
+ -1.65538862e+09,
+  2.85848565e+10,
+  8.81655049e+09,
+  8.81355722e+09,
+  8.81056674e+09,
+  2.73589557e+09,
+  6.09203975e+09,
+  6.19396020e+09,
+  2.95835844e+10,
+ -1.33106776e+09,
+ -1.33061595e+09,
+  2.73496681e+09,
+  6.08997175e+09,
+  6.19185797e+09,
+ -1.33106823e+09,
+  2.95844882e+10,
+ -1.33016438e+09,
+  2.73403805e+09,
+  6.08790386e+09,
+  6.18975528e+09,
+ -1.33061595e+09,
+ -1.33016414e+09,
+  2.95853915e+10,
   4.98723239e+09,
   4.80000023e+09,
   4.61276807e+09,
@@ -311,12 +311,12 @@ const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_initialStrain
 };
 
 const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_stateVarsUpdated[] = {
- -7.89512481e-06,
-  9.60719814e-05,
-  1.00198920e-04,
-  2.35743657e-05,
-  2.45812990e-05,
-  2.55882324e-05,
+ -8.05139365e-06,
+  9.62447955e-05,
+  1.00381728e-04,
+  2.35090854e-05,
+  2.45160409e-05,
+  2.55229964e-05,
   5.53592834e-05,
   6.61856723e-05,
   7.70120612e-05,
diff --git a/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py b/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py
index 9cb4c66..3a0d1b2 100644
--- a/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py
+++ b/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py
@@ -291,7 +291,7 @@ class DruckerPragerPlaneStrainTimeDep(ElasticMaterialApp):
                     strainPPTpdtProd)
       dFac = math.sqrt(2.0) * d
       testMult = 2.0 * ae * am * \
-                 (3.0 * alphaYieldV * meanStrainPPTpdt/am + \
+                 (3.0 * alphaYieldV * (meanStrainPPTpdt/am + meanStressInitial) + \
                   d/(math.sqrt(2.0) * ae) - betaV)/ \
                   (6.0 * alphaYieldV * alphaFlowV * ae + am)
       plasticMult = min(testMult, dFac)
diff --git a/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc b/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc
index 313ff9a..c9a488a 100644
--- a/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc
+++ b/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc
@@ -179,24 +179,24 @@ const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_stra
 };
 
 const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_stress[] = {
- -1.92922457e+07,
- -4.93459240e+06,
- -2.89163969e+06,
+ -1.92892468e+07,
+ -4.94969012e+06,
+ -2.88799503e+06,
   2.05547473e+06,
   2.05547473e+06,
   1.45519152e-11,
 };
 
 const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_elasticConsts[] = {
-  6.78169280e+10,
-  3.08737531e+10,
- -2.64315493e+09,
-  1.98734119e+10,
-  2.43162760e+10,
-  1.57343177e+10,
-  8.93897377e+08,
-  1.00826339e+10,
-  3.98038772e+10,
+  6.78100809e+10,
+  3.08704115e+10,
+ -2.62248144e+09,
+  1.98700731e+10,
+  2.42954488e+10,
+  1.57192741e+10,
+  9.04234592e+08,
+  1.00751121e+10,
+  3.97562357e+10,
   3.37901525e+09,
   3.25490593e+09,
  -1.82383531e+09,
@@ -228,10 +228,10 @@ const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_init
 
 const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_stateVarsUpdated[] = {
   2.30000000e+04,
-  3.64658670e-05,
-  4.84291274e-05,
-  4.75397462e-05,
-  4.27919930e-05,
+  3.62995494e-05,
+  4.86649558e-05,
+  4.76695884e-05,
+  4.27110006e-05,
   5.40000000e+04,
   2.05251755e-04,
   2.16078144e-04,



More information about the CIG-COMMITS mailing list