[cig-commits] r16560 - in short/3D/PyLith/trunk: doc/userguide/materials libsrc/materials unittests/libtests/materials unittests/libtests/materials/data
willic3 at geodynamics.org
willic3 at geodynamics.org
Sun Apr 18 16:56:05 PDT 2010
Author: willic3
Date: 2010-04-18 16:56:03 -0700 (Sun, 18 Apr 2010)
New Revision: 16560
Added:
short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.hh
Modified:
short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx
short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh
Log:
Finished fixing Drucker-Prager tangent matrix and checked in
Drucker-Prager libtests, which now pass.
Modified: short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx 2010-04-17 22:02:34 UTC (rev 16559)
+++ short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx 2010-04-18 23:56:03 UTC (rev 16560)
@@ -4403,12 +4403,16 @@
\end_layout
\begin_layout Standard
-To compute the elastoplastic tangent matrix we proceed in a manner analogous
- to the viscoelastic case, using vector representations of the stress and
- strain tensors, and then computing the derivative of the stress vector
- with respect to the strain vector:
+To compute the elastoplastic tangent matrix we begin by writing
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:105"
+
+\end_inset
+
+ as a single equation in terms of stress and strain vectors:
\begin_inset Formula \begin{equation}
-C_{ij}^{EP}=\frac{\partial\phantom{}^{t+\Delta t}\sigma_{i}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}\frac{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}{\partial\phantom{}^{t+\Delta t}e_{l}}\frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}+R_{i}\frac{\partial\phantom{}^{t+\Delta t}P}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\:,\label{eq:117}\end{equation}
+^{t+\Delta t}\sigma_{i}=\frac{1}{a_{E}}\left(^{t+\Delta t}e_{i}^{\prime}-\Delta e_{i}^{P}\right)+S_{i}^{I}+\frac{R_{i}}{a_{m}}\left(^{t+\Delta t}\theta^{\prime}-\Delta\theta^{P}\right)+P^{I}\label{eq:117}\end{equation}
\end_inset
@@ -4419,170 +4423,117 @@
\end_inset
-First considering the deviatoric portion, we have
+The elastoplastic tangent matrix is then given by
\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}{\partial\phantom{}^{t+\Delta t}e_{l}}=\delta_{kl}\:,\label{eq:119}\end{equation}
+C_{ij}^{EP}=\frac{\partial\phantom{}^{t+\Delta t}\sigma_{i}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{a_{E}}\left(\frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)+\frac{R_{i}}{a_{m}}\left(\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\partial\Delta\theta^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)\:.\label{eq:119}\end{equation}
\end_inset
-and
-\begin_inset Formula \begin{align}
-\frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}} & =\frac{1}{3}\left[\begin{array}{ccc}
-2 & -1 & -1\\
--1 & 2 & -1\\
--1 & -1 & 2\end{array}\right]\:;\; l,j=1,2,3\nonumber \\
-\frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}} & =\delta_{lj}\:;\;\mathrm{otherwise.}\label{eq:120}\end{align}
-\end_inset
+\end_layout
-From equation
+\begin_layout Standard
+From
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:105"
+reference "eq:16"
\end_inset
-, we have
-\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=\frac{1}{a_{E}}\left(\delta_{ik}-\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}\right)\:,\label{eq:121}\end{equation}
-
-\end_inset
-
-and from equation
+ and
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:116"
+reference "eq:106"
\end_inset
- we have
+, we have
\begin_inset Formula \begin{equation}
-\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}\frac{1}{\sqrt{2}\,\phantom{}^{t+\Delta t}d}\left(^{t+\Delta t}e_{i}^{\prime}+a_{E}S_{i}^{I}\right)+\frac{\lambda}{\sqrt{2}}\left[\frac{-1}{^{t+\Delta t}d^{2}}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}\left(^{t+\Delta t}e_{i}^{\prime}+a_{E}S_{i}^{I}\right)+\frac{\delta_{ik}}{^{t+\Delta t}d}\right]\:.\label{eq:122}\end{equation}
+\frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{3}\left[\begin{array}{cccccc}
+2 & -1 & -1 & 0 & 0 & 0\\
+-1 & 2 & -1 & 0 & 0 & 0\\
+-1 & -1 & 2 & 0 & 0 & 0\\
+0 & 0 & 0 & 3 & 0 & 0\\
+0 & 0 & 0 & 0 & 3 & 0\\
+0 & 0 & 0 & 0 & 0 & 3\end{array}\right]\:,\label{eq:120}\end{equation}
\end_inset
-The derivative of
-\begin_inset Formula $^{t+\Delta t}d$
-\end_inset
+and from
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:15"
- is
-\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=\frac{a_{E}T_{k}^{I}+\phantom{}^{t+\Delta t}E_{k}}{\phantom{}^{t+\Delta t}d}\:,\label{eq:123}\end{equation}
-
\end_inset
-where
-\begin_inset Formula \begin{align}
-T_{k}^{I} & =S_{k}^{I}\;\mathrm{and}\;\phantom{}^{t+\Delta t}E_{k}=\phantom{}^{t+\Delta t}e_{k}^{\prime}\:;\; k=1,2,3\nonumber \\
-T_{k}^{I} & =2S_{k}^{I}\;\mathrm{and}\;\phantom{}^{t+\Delta t}E_{k}=2\phantom{}^{t+\Delta t}e_{k}^{\prime}\:;\; k=4,5,6\:.\label{eq:124}\end{align}
+ and
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:106"
\end_inset
-We then need to compute the derivative of
-\begin_inset Formula $^{t+\Delta t}\lambda$
-\end_inset
-
-:
+ we have
\begin_inset Formula \begin{equation}
-\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=\frac{2a_{E}a_{m}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{3\alpha_{f}}{a_{m}}\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{l}}\frac{\partial\phantom{}^{t+\Delta t}\epsilon_{l}}{\partial\phantom{}^{t+\Delta t}e_{m}}\frac{\partial\phantom{}^{t+\Delta t}e_{m}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}+\frac{1}{\sqrt{2}a_{E}}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}\right)\:.\label{eq:125}\end{equation}
+\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{R_{j}}{3}\:.\label{eq:121}\end{equation}
\end_inset
-The first derivative in
+From equation
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:125"
+reference "eq:116"
\end_inset
- is simply
+ we have
\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{l}}=\frac{R_{l}}{3}\:.\label{eq:126}\end{equation}
+\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{\sqrt{2}\,\phantom{}^{t+\Delta t}d}\left[\left(^{t+\Delta t}e_{i}^{\prime}+a_{E}S_{i}^{I}\right)\left(\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\lambda}{\phantom{}^{t+\Delta t}d}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)+\lambda\frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right]\:.\label{eq:122}\end{equation}
\end_inset
-The third derivative yields the same result as equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:119"
-
+The derivative of
+\begin_inset Formula $^{t+\Delta t}d$
\end_inset
-, and the fourth is given by equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:123"
+ is
+\begin_inset Formula \begin{equation}
+\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{a_{E}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\phantom{}^{t+\Delta t}d}\:,\label{eq:123}\end{equation}
\end_inset
-.
- The second derivative is analogous to equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:120"
-
-\end_inset
-
-:
+where
\begin_inset Formula \begin{align}
-\frac{\partial\phantom{}^{t+\Delta t}\epsilon_{l}}{\partial\phantom{}^{t+\Delta t}e_{m}} & =3\left[\begin{array}{ccc}
-\frac{1}{2} & -1 & -1\\
--1 & \frac{1}{2} & -1\\
--1 & -1 & \frac{1}{2}\end{array}\right]\:;\; l,m=1,2,3\nonumber \\
-\frac{\partial\phantom{}^{t+\Delta t}\epsilon_{l}}{\partial\phantom{}^{t+\Delta t}e_{m}} & =\delta_{lm}\:;\;\mathrm{otherwise.}\label{eq:127}\end{align}
+T_{j}^{I} & =S_{j}^{I}\;\mathrm{and}\;\phantom{}^{t+\Delta t}E_{j}=\phantom{}^{t+\Delta t}e_{j}^{\prime}\:;\; j=1,2,3\nonumber \\
+T_{j}^{I} & =2S_{j}^{I}\;\mathrm{and}\;\phantom{}^{t+\Delta t}E_{j}=2\phantom{}^{t+\Delta t}e_{j}^{\prime}\:;\; j=4,5,6\:.\label{eq:124}\end{align}
\end_inset
-
-\end_layout
-
-\begin_layout Standard
-For the volumetric portion we use equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:113"
-
+The derivative of
+\begin_inset Formula $^{t+\Delta t}\lambda$
\end_inset
- to obtain
-\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}P}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{a_{m}}\left(\frac{R_{j}}{3}-\alpha_{g}\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)\:.\label{eq:128}\end{equation}
+ is a function of derivatives already computed:
+\begin_inset Formula \begin{align}
+\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}} & =\frac{2a_{E}a_{m}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{3\alpha_{f}}{a_{m}}\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}+\frac{1}{\sqrt{2}a_{E}}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)\nonumber \\
+ & =\frac{2a_{E}a_{m}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{\alpha_{f}R_{j}}{a_{m}}+\frac{a_{E}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\sqrt{2}a_{E}\phantom{}^{t+\Delta t}d}\right)\:.\label{eq:125}\end{align}
\end_inset
-From equation
+Finally, from
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:114"
+reference "eq:108"
\end_inset
-,
+, the derivative of the volumetric plastic strain increment is:
\begin_inset Formula \begin{equation}
-\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{2a_{E}a_{m}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{R_{j}\alpha_{f}}{a_{m}}+\frac{1}{\sqrt{2}a_{E}}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)\:.\label{eq:129}\end{equation}
+\frac{\partial\Delta\theta^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\alpha_{g}\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\:.\label{eq:126}\end{equation}
\end_inset
-From equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:111"
-\end_inset
-
- we have
-\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{a_{e}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\phantom{}^{t+\Delta t}d}\:,\label{eq:130}\end{equation}
-
-\end_inset
-
-so that
-\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}P}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{a_{m}}\left[\frac{R_{j}}{3}-\frac{2a_{E}a_{m}\alpha_{g}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{R_{j}\alpha_{f}}{a_{m}}+\frac{a_{E}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\sqrt{2}a_{E}\phantom{}^{t+\Delta t}d}\right)\right]\:.\label{eq:131}\end{equation}
-
-\end_inset
-
-
\end_layout
\begin_layout Section
Modified: short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.cc 2010-04-17 22:02:34 UTC (rev 16559)
+++ short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.cc 2010-04-18 23:56:03 UTC (rev 16560)
@@ -211,13 +211,11 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- const double alphaYield =
- 2.0 * sin(frictionAngle)/(sqrt(3.0) * (3.0 - sin(frictionAngle)));
- const double beta =
- 6.0 * cohesion *
- cos(frictionAngle)/(sqrt(3.0) * (3.0 - sin(frictionAngle)));
- const double alphaFlow =
- 2.0 * sin(dilatationAngle)/(sqrt(3.0) * (3.0 - sin(dilatationAngle)));
+ const double denomFriction = sqrt(3.0) * (3.0 - sin(frictionAngle));
+ const double denomDilatation = sqrt(3.0) * (3.0 - sin(dilatationAngle));
+ const double alphaYield = 2.0 * sin(frictionAngle)/denomFriction;
+ const double beta = 6.0 * cohesion * cos(frictionAngle)/denomFriction;
+ const double alphaFlow = 2.0 * sin(dilatationAngle)/denomDilatation;
if (lambda <= 0.0) {
std::ostringstream msg;
@@ -819,41 +817,14 @@
pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
strainPPTpdt) +
strainPPTpdtProd);
- const double plasticMult = 2.0 * ae * am *
- (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+ const double plasticFac = 2.0 * ae * am/
(6.0 * alphaYield * alphaFlow * ae + am);
+ const double meanStrainFac = 3.0 * alphaYield/am;
+ const double dFac = 1.0/(sqrt(2.0) * ae);
+ const double plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
// Define some constants, vectors, and matrices.
- const double vec1[] = {strainPPTpdt[0] + ae * devStressInitial[0],
- strainPPTpdt[1] + ae * devStressInitial[1],
- strainPPTpdt[2] + ae * devStressInitial[2],
- strainPPTpdt[3] + ae * devStressInitial[3],
- strainPPTpdt[4] + ae * devStressInitial[4],
- strainPPTpdt[5] + ae * devStressInitial[5]};
- const double dDdEPrime[] = {vec1[0]/d,
- vec1[1]/d,
- vec1[2]/d,
- 2.0 * vec1[3]/d,
- 2.0 * vec1[4]/d,
- 2.0 * vec1[5]/d};
- const double const1 = 2.0 * ae * am/
- (6.0 * alphaYield * alphaFlow * ae + am);
- const double const2 = 3.0 * alphaYield/am;
- const double const3 = 1.0/(sqrt(2.0) * ae);
- const double dLambdadEPrime[] = {
- const1 * (-1.5 * const2 + const3 * dDdEPrime[0]),
- const1 * (-1.5 * const2 + const3 * dDdEPrime[1]),
- const1 * (-1.5 * const2 + const3 * dDdEPrime[2]),
- const1 * ( const3 * dDdEPrime[3]),
- const1 * ( const3 * dDdEPrime[4]),
- const1 * ( const3 * dDdEPrime[5])};
- const double delta[6][6] = {
- {1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
- {0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
- {0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
- {0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
- {0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
- {0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
const double third = 1.0/3.0;
const double dEdEpsilon[6][6] = {
{ 2.0 * third, -third, -third, 0.0, 0.0, 0.0},
@@ -862,48 +833,46 @@
{ 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
{ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
{ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
+ const double vec1[] = {strainPPTpdt[0] + ae * devStressInitial[0],
+ strainPPTpdt[1] + ae * devStressInitial[1],
+ strainPPTpdt[2] + ae * devStressInitial[2],
+ strainPPTpdt[3] + ae * devStressInitial[3],
+ strainPPTpdt[4] + ae * devStressInitial[4],
+ strainPPTpdt[5] + ae * devStressInitial[5]};
+ const double dDdEpsilon[] = {vec1[0]/d,
+ vec1[1]/d,
+ vec1[2]/d,
+ 2.0 * vec1[3]/d,
+ 2.0 * vec1[4]/d,
+ 2.0 * vec1[5]/d};
+ const double dLambdadEpsilon[] = {
+ plasticFac * (alphaYield/am + dFac * dDdEpsilon[0]),
+ plasticFac * (alphaYield/am + dFac * dDdEpsilon[1]),
+ plasticFac * (alphaYield/am + dFac * dDdEpsilon[2]),
+ plasticFac * ( dFac * dDdEpsilon[3]),
+ plasticFac * ( dFac * dDdEpsilon[4]),
+ plasticFac * ( dFac * dDdEpsilon[5])};
- const double const4 = 1.0/(sqrt(2.0) *d);
- const double const5 = plasticMult/sqrt(2.0);
- const double const6 = -1.0/(d * d);
-
- const double dPressuredEpsilon[] = {
- (third - const1 * (alphaYield/am + const3 * dDdEPrime[0]))/am,
- (third - const1 * (alphaYield/am + const3 * dDdEPrime[1]))/am,
- (third - const1 * (alphaYield/am + const3 * dDdEPrime[2]))/am,
- (-const1 * const3 * dDdEPrime[3])/am,
- (-const1 * const3 * dDdEPrime[4])/am,
- (-const1 * const3 * dDdEPrime[5])/am};
-
- double dDeltaEdEPrime = 0.0;
- double dSdEPrime[tensorSize][tensorSize];
+ const double dFac2 = 1.0/(sqrt(2.0) * d);
+ double dDeltaEdEpsilon = 0.0;
+ // Compute elasticity matrix.
for (int iComp=0; iComp < tensorSize; ++iComp) {
- for (int kComp=0; kComp < tensorSize; ++kComp) {
- dDeltaEdEPrime = const4 * dLambdadEPrime[kComp] * vec1[iComp] +
- const5 * (const6 * dDdEPrime[kComp] * vec1[iComp] +
- delta[iComp][kComp]/d);
- dSdEPrime[iComp][kComp] = (delta[iComp][kComp] - dDeltaEdEPrime)/ae;
- } // for
- } // for
-
- // Matrix multiplication.
- double sum = 0.0;
- double pressureTerm = 0.0;
- double elasticMatrix[tensorSize][tensorSize];
- for (int iComp=0; iComp < tensorSize; ++iComp) {
for (int jComp=0; jComp < tensorSize; ++jComp) {
- pressureTerm = diag[iComp] * dPressuredEpsilon[jComp];
- sum = 0.0;
- for (int kComp=0; kComp < tensorSize; ++kComp) {
- sum += dSdEPrime[iComp][kComp] * dEdEpsilon[jComp][kComp];
- } // for
- elasticMatrix[iComp][jComp] = sum + pressureTerm;
+ int iCount = jComp + tensorSize * iComp;
+ dDeltaEdEpsilon = dFac2 * (vec1[iComp] *
+ (dLambdadEpsilon[jComp] -
+ plasticMult * dDdEpsilon[jComp]/d) +
+ plasticMult * dEdEpsilon[iComp][jComp]);
+ elasticConsts[iCount] = (dEdEpsilon[iComp][jComp] -
+ dDeltaEdEpsilon)/ae +
+ diag[iComp] * (third * diag[jComp] -
+ alphaFlow * dLambdadEpsilon[jComp])/am;
} // for
} // for
-
- PetscLogFlops(161 + tensorSize * tensorSize * (12 + 2 * tensorSize));
+ PetscLogFlops(109 + tensorSize * tensorSize * 15);
+
} else {
// No plastic strain.
const double lambda2mu = lambda + mu2;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2010-04-17 22:02:34 UTC (rev 16559)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2010-04-18 23:56:03 UTC (rev 16560)
@@ -33,6 +33,7 @@
TestMaxwellPlaneStrain.cc \
TestGenMaxwellIsotropic3D.cc \
TestPowerLaw3D.cc \
+ TestDruckerPrager3D.cc \
TestEffectiveStress.cc \
test_materials.cc
@@ -50,6 +51,7 @@
TestMaxwellIsotropic3D.hh \
TestMaxwellPlaneStrain.hh \
TestPowerLaw3D.hh \
+ TestDruckerPrager3D.hh \
TestEffectiveStress.hh
# Source files associated with testing data
@@ -68,7 +70,9 @@
data/GenMaxwellIsotropic3DElasticData.cc \
data/GenMaxwellIsotropic3DTimeDepData.cc \
data/PowerLaw3DElasticData.cc \
- data/PowerLaw3DTimeDepData.cc
+ data/PowerLaw3DTimeDepData.cc \
+ data/DruckerPrager3DElasticData.cc \
+ data/DruckerPrager3DTimeDepData.cc
noinst_HEADERS += \
@@ -87,6 +91,8 @@
data/MaxwellPlaneStrainTimeDepData.hh \
data/PowerLaw3DElasticData.hh \
data/PowerLaw3DTimeDepData.hh \
+ data/DruckerPrager3DElasticData.hh \
+ data/DruckerPrager3DTimeDepData.hh \
data/header.hh
AM_CPPFLAGS = \
Added: short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,215 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestDruckerPrager3D.hh" // Implementation of class methods
+
+#include "data/DruckerPrager3DElasticData.hh" // USES DruckerPrager3DElasticData
+#include "data/DruckerPrager3DTimeDepData.hh" // USES DruckerPrager3DTimeDepData
+
+#include "pylith/materials/DruckerPrager3D.hh" // USES DruckerPrager3D
+
+#include <cstring> // USES memcpy()
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestDruckerPrager3D );
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::materials::TestDruckerPrager3D::setUp(void)
+{ // setUp
+ _material = new DruckerPrager3D();
+ _matElastic = new DruckerPrager3D();
+ _data = new DruckerPrager3DElasticData();
+ _dataElastic = new DruckerPrager3DElasticData();
+ setupNormalizer();
+} // setUp
+
+// ----------------------------------------------------------------------
+// Test timeStep()
+void
+pylith::materials::TestDruckerPrager3D::testTimeStep(void)
+{ // testTimeStep
+ DruckerPrager3D material;
+
+ CPPUNIT_ASSERT_EQUAL(false, material._needNewJacobian);
+
+ const double dt1 = 1.0;
+ material.timeStep(dt1);
+ CPPUNIT_ASSERT_EQUAL(dt1, material.Material::timeStep());
+ CPPUNIT_ASSERT_EQUAL(true, material.needNewJacobian());
+
+ const double dt2 = 2.0;
+ material.timeStep(dt2);
+ CPPUNIT_ASSERT_EQUAL(dt2, material.Material::timeStep());
+ CPPUNIT_ASSERT_EQUAL(true, material.needNewJacobian());
+} // testTimeStep
+
+// ----------------------------------------------------------------------
+// Test useElasticBehavior()
+void
+pylith::materials::TestDruckerPrager3D::testUseElasticBehavior(void)
+{ // testUseElasticBehavior
+ DruckerPrager3D material;
+
+ material.useElasticBehavior(true);
+ // Some compilers/operating systems (cygwin) don't allow comparing
+ // pointers. Use first test to determine if we can compare pointers.
+ if (&pylith::materials::DruckerPrager3D::_calcStressElastic ==
+ material._calcStressFn) {
+ CPPUNIT_ASSERT(&pylith::materials::DruckerPrager3D::_calcStressElastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(&pylith::materials::DruckerPrager3D::_calcElasticConstsElastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(&pylith::materials::DruckerPrager3D::_updateStateVarsElastic ==
+ material._updateStateVarsFn);
+
+ material.useElasticBehavior(false);
+ CPPUNIT_ASSERT(&pylith::materials::DruckerPrager3D::_calcStressElastoplastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(&pylith::materials::DruckerPrager3D::_calcElasticConstsElastoplastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(&pylith::materials::DruckerPrager3D::_updateStateVarsElastoplastic ==
+ material._updateStateVarsFn);
+ } // if
+} // testUseElasticBehavior
+
+// ----------------------------------------------------------------------
+// Test usesHasStateVars()
+void
+pylith::materials::TestDruckerPrager3D::testHasStateVars(void)
+{ // testHasStateVars
+ DruckerPrager3D material;
+ CPPUNIT_ASSERT_EQUAL(true, material.hasStateVars());
+} // testHasStateVars
+
+// ----------------------------------------------------------------------
+// Test _calcStressElastic()
+void
+pylith::materials::TestDruckerPrager3D::test_calcStressElastic(void)
+{ // test_calcStressElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_calcStress();
+} // test_calcStressElastic
+
+// ----------------------------------------------------------------------
+// Test calcElasticConstsElastic()
+void
+pylith::materials::TestDruckerPrager3D::test_calcElasticConstsElastic(void)
+{ // test_calcElasticConstsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_calcElasticConsts();
+} // test_calcElasticConstsElastic
+
+// ----------------------------------------------------------------------
+// Test _updateStateVarsElastic()
+void
+pylith::materials::TestDruckerPrager3D::test_updateStateVarsElastic(void)
+{ // test_updateStateVarsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_updateStateVars();
+} // test_updateStateVarsElastic
+
+// ----------------------------------------------------------------------
+// Test _calcStressTimeDep()
+void
+pylith::materials::TestDruckerPrager3D::test_calcStressTimeDep(void)
+{ // test_calcStressTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new DruckerPrager3DTimeDepData();
+
+ double dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcStress();
+} // test_calcStressTimeDep
+
+// ----------------------------------------------------------------------
+// Test _calcElasticConstsTimeDep()
+void
+pylith::materials::TestDruckerPrager3D::test_calcElasticConstsTimeDep(void)
+{ // test_calcElasticConstsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new DruckerPrager3DTimeDepData();
+
+ double dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcElasticConsts();
+} // test_calcElasticConstsTimeDep
+
+// ----------------------------------------------------------------------
+// Test _updateStateVarsTimeDep()
+void
+pylith::materials::TestDruckerPrager3D::test_updateStateVarsTimeDep(void)
+{ // test_updateStateVarsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new DruckerPrager3DTimeDepData();
+
+ double dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_updateStateVars();
+
+} // test_updateStateVarsTimeDep
+
+// ----------------------------------------------------------------------
+// Test _stableTimeStepImplicit()
+void
+pylith::materials::TestDruckerPrager3D::test_stableTimeStepImplicit(void)
+{ // test_stableTimeStepImplicit
+ CPPUNIT_ASSERT(0 != _matElastic);
+
+ delete _dataElastic; _dataElastic = new DruckerPrager3DTimeDepData();
+
+ TestElasticMaterial::test_stableTimeStepImplicit();
+} // test_stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Test hasProperty().
+void
+pylith::materials::TestDruckerPrager3D::testHasProperty(void)
+{ // testHasProperty
+ DruckerPrager3D material;
+
+ CPPUNIT_ASSERT(material.hasProperty("mu"));
+ CPPUNIT_ASSERT(material.hasProperty("lambda"));
+ CPPUNIT_ASSERT(material.hasProperty("density"));
+ CPPUNIT_ASSERT(material.hasProperty("alpha_yield"));
+ CPPUNIT_ASSERT(material.hasProperty("beta"));
+ CPPUNIT_ASSERT(material.hasProperty("alpha_flow"));
+ CPPUNIT_ASSERT(!material.hasProperty("aaa"));
+} // testHasProperty
+
+// ----------------------------------------------------------------------
+// Test hasStateVar().
+void
+pylith::materials::TestDruckerPrager3D::testHasStateVar(void)
+{ // testHasStateVar
+ DruckerPrager3D material;
+
+ CPPUNIT_ASSERT(material.hasStateVar("plastic_strain"));
+} // testHasStateVar
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,117 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/materials/TestDruckerPrager3D.hh
+ *
+ * @brief C++ TestDruckerPrager3D object
+ *
+ * C++ unit testing for DruckerPrager3D.
+ */
+
+#if !defined(pylith_materials_testdruckerprager3d_hh)
+#define pylith_materials_testdruckerprager3d_hh
+
+#include "TestElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class TestDruckerPrager3D;
+ class DruckerPrager3DElasticData;
+ class DruckerPrager3DTimeDepData;
+ } // materials
+} // pylith
+
+/// C++ unit testing for DruckerPrager3D
+class pylith::materials::TestDruckerPrager3D : public TestElasticMaterial
+{ // class TestDruckerPrager3D
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestDruckerPrager3D );
+
+ CPPUNIT_TEST( testDimension );
+ CPPUNIT_TEST( testTensorSize );
+ CPPUNIT_TEST( testDBToProperties );
+ CPPUNIT_TEST( testNonDimProperties );
+ CPPUNIT_TEST( testDimProperties );
+ CPPUNIT_TEST( testDBToStateVars );
+ CPPUNIT_TEST( testNonDimStateVars );
+ CPPUNIT_TEST( testDimStateVars );
+ CPPUNIT_TEST( test_calcDensity );
+ CPPUNIT_TEST( test_stableTimeStepImplicit );
+
+ // Need to test Drucker-Prager elastoplastic specific behavior.
+ CPPUNIT_TEST( testTimeStep );
+ CPPUNIT_TEST( testUseElasticBehavior );
+ CPPUNIT_TEST( testHasStateVars );
+
+ CPPUNIT_TEST( test_calcStressElastic );
+ CPPUNIT_TEST( test_calcStressTimeDep );
+ CPPUNIT_TEST( test_calcElasticConstsElastic );
+ CPPUNIT_TEST( test_calcElasticConstsTimeDep );
+ CPPUNIT_TEST( test_updateStateVarsElastic );
+ CPPUNIT_TEST( test_updateStateVarsTimeDep );
+
+ CPPUNIT_TEST( testHasProperty );
+ CPPUNIT_TEST( testHasStateVar );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+ /// Test timeStep()
+ void testTimeStep(void);
+
+ /// Test useElasticBehavior()
+ void testUseElasticBehavior(void);
+
+ /// Test hasStateVars()
+ void testHasStateVars(void);
+
+ /// Test _calcStressElastic()
+ void test_calcStressElastic(void);
+
+ /// Test _calcElasticConstsElastic()
+ void test_calcElasticConstsElastic(void);
+
+ /// Test _updateStateVarsElastic()
+ void test_updateStateVarsElastic(void);
+
+ /// Test _calcStressTimeDep()
+ void test_calcStressTimeDep(void);
+
+ /// Test _calcElasticConstsTimeDep()
+ void test_calcElasticConstsTimeDep(void);
+
+ /// Test _updateStatevarsTimeDep()
+ void test_updateStateVarsTimeDep(void);
+
+ /// Test _stableTimeStepImplicit()
+ void test_stableTimeStepImplicit(void);
+
+ /// Test hasProperty()
+ void testHasProperty(void);
+
+ /// Test hasStateVar()
+ void testHasStateVar(void);
+
+}; // class TestDruckerPrager3D
+
+#endif // pylith_materials_testdruckerprager3d_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElastic.py (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElastic.py 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,245 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/DruckerPrager3DElastic.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ DruckerPrager3D object with elastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+import math
+
+# ----------------------------------------------------------------------
+dimension = 3
+numElasticConsts = 36
+tensorSize = 6
+
+# DruckerPrager3DElastic class
+class DruckerPrager3DElastic(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ DruckerPrager3D object with elastic behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="druckerprager3delastic"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ # import pdb
+ # pdb.set_trace()
+ numLocs = 2
+
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "friction_angle", "cohesion",
+ "dilatation_angle"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+
+ self.dbStateVarValues = ["plastic-strain-xx",
+ "plastic-strain-yy",
+ "plastic-strain-zz",
+ "plastic-strain-xy",
+ "plastic-strain-yz",
+ "plastic-strain-xz"
+ ]
+ self.numStateVarValues = numpy.array([6], dtype=numpy.int32)
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ # First case has different values for friction angle and dilatation angle.
+ frictionAngleA = math.radians(30.0)
+ dilatationAngleA = math.radians(20.0)
+ cohesionA = 3.0e5
+ strainA = [1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4]
+ initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
+ initialStrainA = [3.1e-4, 3.2e-4, 3.3e-4, 3.4e-4, 3.5e-4, 3.6e-4]
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+
+ denomFrictionA = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleA))
+ denomDilatationA = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleA))
+ alphaYieldA = 2.0 * math.sin(frictionAngleA)/denomFrictionA
+ betaA = 6.0 * cohesionA * math.cos(frictionAngleA)/denomFrictionA
+ alphaFlowA = 2.0 * math.sin(dilatationAngleA)/denomDilatationA
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ # Second case has same values for friction angle and dilatation angle.
+ frictionAngleB = math.radians(25.0)
+ dilatationAngleB = math.radians(25.0)
+ cohesionB = 1.0e5
+ 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]
+ initialStrainB = [6.1e-4, 6.2e-4, 6.3e-4, 6.4e-4, 6.5e-4, 6.6e-4]
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ denomFrictionB = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleB))
+ denomDilatationB = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleB))
+ alphaYieldB = 2.0 * math.sin(frictionAngleB)/denomFrictionB
+ betaB = 6.0 * cohesionB * math.cos(frictionAngleB)/denomFrictionB
+ alphaFlowB = 2.0 * math.sin(dilatationAngleB)/denomDilatationB
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.densityScale = 1.0e+3
+ self.timeScale = 1.0
+
+ self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
+ frictionAngleA, cohesionA, \
+ dilatationAngleA],
+ [densityB, vsB, vpB, \
+ frictionAngleB, cohesionB, \
+ dilatationAngleB] ],
+ dtype=numpy.float64)
+ self.properties = numpy.array([ [densityA, muA, lambdaA, \
+ alphaYieldA, betaA, \
+ alphaFlowA],
+ [densityB, muB, lambdaB, \
+ alphaYieldB, betaB, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ self.dbStateVars = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.stateVars = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, \
+ alphaYieldA, betaA/mu0, \
+ alphaFlowA],
+ [densityB/density0, muB/mu0, lambdaB/mu0, \
+ alphaYieldB, betaB/mu0, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ self.stateVarsNondim = self.stateVars # no scaling
+
+ self.initialStress = numpy.array([initialStressA,
+ initialStressB],
+ dtype=numpy.float64)
+ self.initialStrain = numpy.array([initialStrainA,
+ initialStrainB],
+ dtype=numpy.float64)
+
+ self.density = numpy.array([densityA,
+ densityB],
+ dtype=numpy.float64)
+
+ self.strain = numpy.array([strainA,
+ strainB],
+ dtype=numpy.float64)
+
+ self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts), \
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA, \
+ initialStressA, initialStrainA)
+ (self.elasticConsts[1,:], self.stress[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB, \
+ initialStressB, initialStrainB)
+
+ self.dtStableImplicit = 1.0e10
+
+ plasticStrainUpdated = numpy.zeros((numLocs, tensorSize),
+ dtype=numpy.float64)
+
+ self.stateVarsUpdated = numpy.array( [plasticStrainUpdated[0,:],
+ plasticStrainUpdated[1,:]],
+ dtype=numpy.float64)
+
+ return
+
+
+ def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV):
+ """
+ Compute stress and derivative of elasticity matrix.
+ """
+ C1111 = lambdaV + 2.0*muV
+ C1122 = lambdaV
+ C1133 = lambdaV
+ C1112 = 0.0
+ C1123 = 0.0
+ C1113 = 0.0
+ C2211 = lambdaV
+ C2222 = lambdaV + 2.0*muV
+ C2233 = lambdaV
+ C2212 = 0.0
+ C2223 = 0.0
+ C2213 = 0.0
+ C3311 = lambdaV
+ C3322 = lambdaV
+ C3333 = lambdaV + 2.0*muV
+ C3312 = 0.0
+ C3323 = 0.0
+ C3313 = 0.0
+ C1211 = 0.0
+ C1222 = 0.0
+ C1233 = 0.0
+ C1212 = 2.0*muV
+ C1223 = 0.0
+ C1213 = 0.0
+ C2311 = 0.0
+ C2322 = 0.0
+ C2333 = 0.0
+ C2312 = 0.0
+ C2323 = 2.0*muV
+ C2313 = 0.0
+ C1311 = 0.0
+ C1322 = 0.0
+ C1333 = 0.0
+ C1312 = 0.0
+ C1323 = 0.0
+ C1313 = 2.0*muV
+ elasticConsts = numpy.array([C1111, C1122, C1133, C1112, C1123, C1113,
+ C2211, C2222, C2233, C2212, C2223, C2213,
+ C3311, C3322, C3333, C3312, C3323, C3313,
+ C1211, C1222, C1233, C1212, C1223, C1213,
+ C2311, C2322, C2333, C2312, C2323, C2313,
+ C1311, C1322, C1333, C1312, C1323, C1313],
+ dtype=numpy.float64)
+
+ strain = numpy.reshape(strainV, (6,1))
+ initialStress = numpy.reshape(initialStressV, (tensorSize,1))
+ initialStrain = numpy.reshape(initialStrainV, (tensorSize,1))
+ elastic = numpy.array([ [C1111, C1122, C1133, C1112, C1123, C1113],
+ [C2211, C2222, C2233, C2212, C2223, C2213],
+ [C3311, C3322, C3333, C3312, C3323, C3313],
+ [C1211, C1222, C1233, C1212, C1223, C1213],
+ [C2311, C2322, C2333, C2312, C2323, C2313],
+ [C1311, C1322, C1333, C1312, C1323, C1313] ],
+ dtype=numpy.float64)
+ stress = numpy.dot(elastic, strain-initialStrain) + initialStress
+ return (elasticConsts, numpy.ravel(stress))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = DruckerPrager3DElastic()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.cc 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,358 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerprager3delastic.
+
+#include "DruckerPrager3DElasticData.hh"
+
+const int pylith::materials::DruckerPrager3DElasticData::_dimension = 3;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numLocs = 2;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numProperties = 6;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numStateVars = 1;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numDBProperties = 6;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numDBStateVars = 6;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numPropsQuadPt = 6;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numVarsQuadPt = 6;
+
+const double pylith::materials::DruckerPrager3DElasticData::_lengthScale = 1.00000000e+03;
+
+const double pylith::materials::DruckerPrager3DElasticData::_timeScale = 1.00000000e+00;
+
+const double pylith::materials::DruckerPrager3DElasticData::_pressureScale = 2.25000000e+10;
+
+const double pylith::materials::DruckerPrager3DElasticData::_densityScale = 1.00000000e+03;
+
+const double pylith::materials::DruckerPrager3DElasticData::_dtStableImplicit = 1.00000000e+10;
+
+const int pylith::materials::DruckerPrager3DElasticData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::DruckerPrager3DElasticData::_numStateVarValues[] = {
+6,
+};
+
+const char* pylith::materials::DruckerPrager3DElasticData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"friction_angle",
+"cohesion",
+"dilatation_angle",
+};
+
+const char* pylith::materials::DruckerPrager3DElasticData::_dbStateVarValues[] = {
+"plastic-strain-xx",
+"plastic-strain-yy",
+"plastic-strain-zz",
+"plastic-strain-xy",
+"plastic-strain-yz",
+"plastic-strain-xz",
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 5.23598776e-01,
+ 3.00000000e+05,
+ 3.49065850e-01,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 4.36332313e-01,
+ 1.00000000e+05,
+ 4.36332313e-01,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_dbStateVars[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 2.30940108e-01,
+ 3.60000000e+05,
+ 1.48583084e-01,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 1.89338478e-01,
+ 1.21811303e+05,
+ 1.89338478e-01,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_stateVars[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 2.30940108e-01,
+ 1.60000000e-05,
+ 1.48583084e-01,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 1.89338478e-01,
+ 5.41383567e-06,
+ 1.89338478e-01,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_stateVarsNondim[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_strain[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.30000000e-04,
+ 1.40000000e-04,
+ 1.50000000e-04,
+ 1.60000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.30000000e-04,
+ 4.40000000e-04,
+ 4.50000000e-04,
+ 4.60000000e-04,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_stress[] = {
+ -2.24790000e+07,
+ -2.24780000e+07,
+ -2.24770000e+07,
+ -8.97600000e+06,
+ -8.97500000e+06,
+ -8.97400000e+06,
+ -2.82900000e+06,
+ -2.82800000e+06,
+ -2.82700000e+06,
+ -1.09800000e+06,
+ -1.09700000e+06,
+ -1.09600000e+06,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_elasticConsts[] = {
+ 6.75000000e+10,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.25000000e+10,
+ 6.75000000e+10,
+ 2.25000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 6.75000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 8.64000000e+09,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.88000000e+09,
+ 8.64000000e+09,
+ 2.88000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 8.64000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.76000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.76000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.76000000e+09,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.30000000e+04,
+ 2.40000000e+04,
+ 2.50000000e+04,
+ 2.60000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.30000000e+04,
+ 5.40000000e+04,
+ 5.50000000e+04,
+ 5.60000000e+04,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_initialStrain[] = {
+ 3.10000000e-04,
+ 3.20000000e-04,
+ 3.30000000e-04,
+ 3.40000000e-04,
+ 3.50000000e-04,
+ 3.60000000e-04,
+ 6.10000000e-04,
+ 6.20000000e-04,
+ 6.30000000e-04,
+ 6.40000000e-04,
+ 6.50000000e-04,
+ 6.60000000e-04,
+};
+
+const double pylith::materials::DruckerPrager3DElasticData::_stateVarsUpdated[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+pylith::materials::DruckerPrager3DElasticData::DruckerPrager3DElasticData(void)
+{ // constructor
+ dimension = _dimension;
+ numLocs = _numLocs;
+ numProperties = _numProperties;
+ numStateVars = _numStateVars;
+ numDBProperties = _numDBProperties;
+ numDBStateVars = _numDBStateVars;
+ numPropsQuadPt = _numPropsQuadPt;
+ numVarsQuadPt = _numVarsQuadPt;
+ lengthScale = _lengthScale;
+ timeScale = _timeScale;
+ pressureScale = _pressureScale;
+ densityScale = _densityScale;
+ dtStableImplicit = _dtStableImplicit;
+ numPropertyValues = const_cast<int*>(_numPropertyValues);
+ numStateVarValues = const_cast<int*>(_numStateVarValues);
+ dbPropertyValues = const_cast<char**>(_dbPropertyValues);
+ dbStateVarValues = const_cast<char**>(_dbStateVarValues);
+ dbProperties = const_cast<double*>(_dbProperties);
+ dbStateVars = const_cast<double*>(_dbStateVars);
+ properties = const_cast<double*>(_properties);
+ stateVars = const_cast<double*>(_stateVars);
+ propertiesNondim = const_cast<double*>(_propertiesNondim);
+ stateVarsNondim = const_cast<double*>(_stateVarsNondim);
+ density = const_cast<double*>(_density);
+ strain = const_cast<double*>(_strain);
+ stress = const_cast<double*>(_stress);
+ elasticConsts = const_cast<double*>(_elasticConsts);
+ initialStress = const_cast<double*>(_initialStress);
+ initialStrain = const_cast<double*>(_initialStrain);
+ stateVarsUpdated = const_cast<double*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::DruckerPrager3DElasticData::~DruckerPrager3DElasticData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DElasticData.hh 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,104 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerprager3delastic.
+
+#if !defined(pylith_materials_druckerprager3delasticdata_hh)
+#define pylith_materials_druckerprager3delasticdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class DruckerPrager3DElasticData;
+ } // pylith
+} // materials
+
+class pylith::materials::DruckerPrager3DElasticData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ DruckerPrager3DElasticData(void);
+
+ /// Destructor
+ ~DruckerPrager3DElasticData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numLocs;
+
+ static const int _numProperties;
+
+ static const int _numStateVars;
+
+ static const int _numDBProperties;
+
+ static const int _numDBStateVars;
+
+ static const int _numPropsQuadPt;
+
+ static const int _numVarsQuadPt;
+
+ static const double _lengthScale;
+
+ static const double _timeScale;
+
+ static const double _pressureScale;
+
+ static const double _densityScale;
+
+ static const double _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const double _dbProperties[];
+
+ static const double _dbStateVars[];
+
+ static const double _properties[];
+
+ static const double _stateVars[];
+
+ static const double _propertiesNondim[];
+
+ static const double _stateVarsNondim[];
+
+ static const double _density[];
+
+ static const double _strain[];
+
+ static const double _stress[];
+
+ static const double _elasticConsts[];
+
+ static const double _initialStress[];
+
+ static const double _initialStrain[];
+
+ static const double _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_druckerprager3delasticdata_hh
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,348 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ DruckerPrager3D object with time dependent behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+import math
+
+# ----------------------------------------------------------------------
+dimension = 3
+numElasticConsts = 36
+tensorSize = 6
+
+# DruckerPrager3DTimeDep class
+class DruckerPrager3DTimeDep(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ DruckerPrager3D object with time dependent behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="druckerprager3dtimedep"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ # import pdb
+ # pdb.set_trace()
+ numLocs = 2
+
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "friction_angle", "cohesion",
+ "dilatation_angle"]
+ self.propertyValues = ["density", "mu", "lambda",
+ "alpha_yield", "beta", "alpha_flow"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+
+ self.dbStateVarValues = ["plastic-strain-xx",
+ "plastic-strain-yy",
+ "plastic-strain-zz",
+ "plastic-strain-xy",
+ "plastic-strain-yz",
+ "plastic-strain-xz"
+ ]
+ self.stateVarValues = ["plastic-strain"]
+ self.numStateVarValues = numpy.array([6], dtype=numpy.int32)
+
+ self.dt = 2.0e5
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ # First case has same values for friction angle and dilatation angle.
+ frictionAngleA = math.radians(30.0)
+ dilatationAngleA = math.radians(20.0)
+ cohesionA = 3.0e5
+ strainA = [1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4]
+ initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
+ initialStrainA = [3.6e-5, 3.5e-5, 3.4e-5, 3.3e-5, 3.2e-5, 3.1e-5]
+ # initialStressA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ # initialStrainA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+
+ denomFrictionA = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleA))
+ denomDilatationA = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleA))
+ alphaYieldA = 2.0 * math.sin(frictionAngleA)/denomFrictionA
+ betaA = 6.0 * cohesionA * math.cos(frictionAngleA)/denomFrictionA
+ alphaFlowA = 2.0 * math.sin(dilatationAngleA)/denomDilatationA
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ # Second case has different values for friction angle and dilatation angle.
+ frictionAngleB = math.radians(25.0)
+ dilatationAngleB = math.radians(25.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]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.3e-5, 6.6e-5, 6.5e-5, 6.4e-5]
+ # initialStressB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ # initialStrainB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ denomFrictionB = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleB))
+ denomDilatationB = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleB))
+ alphaYieldB = 2.0 * math.sin(frictionAngleB)/denomFrictionB
+ betaB = 6.0 * cohesionB * math.cos(frictionAngleB)/denomFrictionB
+ alphaFlowB = 2.0 * math.sin(dilatationAngleB)/denomDilatationB
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.densityScale = 1.0e+3
+ self.timeScale = 1.0
+
+ self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
+ frictionAngleA, cohesionA, \
+ dilatationAngleA],
+ [densityB, vsB, vpB, \
+ frictionAngleB, cohesionB, \
+ dilatationAngleB] ],
+ dtype=numpy.float64)
+ self.properties = numpy.array([ [densityA, muA, lambdaA, \
+ alphaYieldA, betaA, \
+ alphaFlowA],
+ [densityB, muB, lambdaB, \
+ alphaYieldB, betaB, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ self.dbStateVars = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, \
+ alphaYieldA, betaA/mu0, \
+ alphaFlowA],
+ [densityB/density0, muB/mu0, lambdaB/mu0, \
+ alphaYieldB, betaB/mu0, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ self.initialStress = numpy.array([initialStressA,
+ initialStressB],
+ dtype=numpy.float64)
+ self.initialStrain = numpy.array([initialStrainA,
+ initialStrainB],
+ dtype=numpy.float64)
+
+ self.density = numpy.array([densityA,
+ densityB],
+ dtype=numpy.float64)
+
+ # Define state variables
+ plasStrainA = [4.1e-5, 4.2e-5, 4.3e-5, 4.4e-5, 4.5e-5, 4.6e-5]
+ plasStrainB = [1.1e-5, 1.2e-5, 1.3e-5, 1.4e-5, 1.5e-5, 1.6e-5]
+ # plasStrainA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ # plasStrainB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ self.stateVars = numpy.array([[plasStrainA], [plasStrainB] ],
+ dtype=numpy.float64)
+ self.stateVarsNondim = self.stateVars # No scaling.
+ self.strain = numpy.array([strainA,
+ strainB],
+ dtype=numpy.float64)
+
+ self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize),
+ dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts), \
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA,
+ alphaYieldA, betaA, alphaFlowA,
+ plasStrainA,
+ initialStressA, initialStrainA)
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB,
+ alphaYieldB, betaB, alphaFlowB,
+ plasStrainB,
+ initialStressB, initialStrainB)
+
+ self.dtStableImplicit = 1.0e10
+
+ return
+
+
+ def _scalarProduct(self, tensor1, tensor2):
+ """
+ Compute the scalar product of two tensors stored in vector form.
+ """
+ scalarProduct = tensor1[0] * tensor2[0] + \
+ tensor1[1] * tensor2[1] + \
+ tensor1[2] * tensor2[2] + \
+ 2.0 * (tensor1[3] * tensor2[3] + \
+ tensor1[4] * tensor2[4] + \
+ tensor1[5] * tensor2[5])
+ return scalarProduct
+
+
+ def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
+ muV, lambdaV, alphaYieldV, betaV, alphaFlowV,
+ plasStrainT, initialStress, initialStrain):
+ """
+ Function to compute a particular stress component as a function of a
+ strain component.
+ """
+ strainTest = numpy.array(strainTpdt, dtype=numpy.float64)
+ strainTest[strainComp] = strainVal
+ stressTpdt, visStrainTpdt = self._computeStress(strainTest, muV, lambdaV,
+ alphaYieldV, betaV,
+ alphaFlowV,
+ plasStrainT,
+ initialStress,
+ initialStrain)
+ return stressTpdt[stressComp]
+
+
+ def _computeStress(self, strainTpdt, muV, lambdaV, alphaYieldV, betaV,
+ alphaFlowV, plasStrainT, initialStress, initialStrain):
+ """
+ Function to compute stresses and plastic strains.
+ """
+
+ # Constants
+ mu2 = 2.0 * muV
+ lamPlusMu = lambdaV + muV
+ bulkModulus = lambdaV + mu2/3.0
+ ae = 1.0/mu2
+ am = 1.0/(3.0 * bulkModulus)
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0], dtype=numpy.float64)
+
+ # Values from previous time step
+ meanPlasStrainT = (plasStrainT[0] + plasStrainT[1] + plasStrainT[2])/3.0
+ devPlasStrainT = plasStrainT - meanPlasStrainT * diag
+
+ # Initial stress values
+ meanStressInitial = (initialStress[0] + initialStress[1] +
+ initialStress[2])/3.0
+ devStressInitial = initialStress - meanStressInitial * diag
+
+ # Initial strain values
+ meanStrainInitial = (initialStrain[0] + initialStrain[1] +
+ initialStrain[2])/3.0
+ devStrainInitial = initialStrain - meanStrainInitial * diag
+
+ # Values for current time step
+ meanStrainTpdt = (strainTpdt[0] + strainTpdt[1] + strainTpdt[2])/3.0
+ meanStrainPPTpdt = meanStrainTpdt - meanPlasStrainT - meanStrainInitial
+ strainPPTpdt = strainTpdt - meanStrainTpdt * diag - \
+ devPlasStrainT - devStrainInitial
+
+ # Compute trial elastic stresses and yield function to see if yield should
+ # occur.
+ trialDevStress = strainPPTpdt/ae + devStressInitial
+ trialMeanStress = meanStrainPPTpdt/am + meanStressInitial
+ yieldFunction = 3.0 * alphaYieldV * trialMeanStress + \
+ self._scalarProduct(trialDevStress, trialDevStress) - betaV
+
+ # If yield function is greater than zero, compute elastoplastic stress.
+ if (yieldFunction >= 0.0):
+ devStressInitialProd = self._scalarProduct(devStressInitial,
+ devStressInitial)
+ strainPPTpdtProd = self._scalarProduct(strainPPTpdt, strainPPTpdt)
+ d = math.sqrt(ae * ae * devStressInitialProd + \
+ 2.0 * ae * self._scalarProduct(devStressInitial, \
+ strainPPTpdt) + \
+ strainPPTpdtProd)
+ plasticMult = 2.0 * ae * am * \
+ (3.0 * alphaYieldV * meanStrainPPTpdt/am + \
+ d/(math.sqrt(2.0) * ae) - betaV)/ \
+ (6.0 * alphaYieldV * alphaFlowV * ae + am)
+ deltaMeanPlasticStrain = plasticMult * alphaFlowV
+ meanStressTpdt = (meanStrainPPTpdt - deltaMeanPlasticStrain)/am + \
+ meanStressInitial
+
+ deltaDevPlasticStrain = 0.0
+ stressTpdt = numpy.zeros( (tensorSize), dtype=numpy.float64)
+ plasStrainTpdt = numpy.zeros( (tensorSize), dtype=numpy.float64)
+
+ for iComp in range(tensorSize):
+ deltaDevPlasticStrain = plasticMult * (strainPPTpdt[iComp] + \
+ ae * devStressInitial[iComp])/ \
+ (math.sqrt(2.0) * d)
+ devStressTpdt = (strainPPTpdt[iComp] - deltaDevPlasticStrain)/ae + \
+ devStressInitial[iComp]
+ stressTpdt[iComp] = devStressTpdt + diag[iComp] * meanStressTpdt
+ plasStrainTpdt[iComp] = plasStrainT[iComp] + deltaDevPlasticStrain + \
+ diag[iComp] * deltaMeanPlasticStrain
+
+ return stressTpdt, plasStrainTpdt
+
+
+ def _calcStress(self, strainV, muV, lambdaV, alphaYieldV, betaV, alphaFlowV,
+ plasStrainV, initialStressV, initialStrainV):
+ """
+ Compute stress, updated state variables and derivative of elasticity matrix.
+ """
+ import scipy.misc
+
+ # Define some numpy arrays
+ strainTpdt = numpy.array(strainV, dtype=numpy.float64)
+ plasStrainT = numpy.array(plasStrainV, dtype=numpy.float64)
+ initialStress = numpy.array(initialStressV, dtype=numpy.float64)
+ initialStrain = numpy.array(initialStrainV, dtype=numpy.float64)
+ stressTpdt, plasStrainTpdt = self._computeStress(strainTpdt, muV, lambdaV,
+ alphaYieldV, betaV,
+ alphaFlowV, plasStrainT,
+ initialStress,
+ initialStrain)
+ stateVarsUpdated = numpy.array( [plasStrainTpdt], dtype=numpy.float64)
+ # Compute components of tangent constitutive matrix using numerical
+ # derivatives.
+ derivDx = 1.0e-12
+ derivOrder = 3
+ elasticConstsList = []
+
+ for stressComp in range(tensorSize):
+ for strainComp in range(tensorSize):
+ dStressDStrain = scipy.misc.derivative(self._calcStressComponent,
+ strainTpdt[strainComp],
+ dx=derivDx,
+ args=(strainComp,
+ stressComp,
+ strainTpdt, muV, lambdaV,
+ alphaYieldV, betaV,
+ alphaFlowV,
+ plasStrainT,
+ initialStress,
+ initialStrain),
+ order=derivOrder)
+ elasticConstsList.append(dStressDStrain)
+
+ elasticConsts = numpy.array(elasticConstsList, dtype=numpy.float64)
+
+ return (elasticConsts, numpy.ravel(stressTpdt),
+ numpy.ravel(stateVarsUpdated))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = DruckerPrager3DTimeDep()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,358 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerprager3dtimedep.
+
+#include "DruckerPrager3DTimeDepData.hh"
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_dimension = 3;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numLocs = 2;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numProperties = 6;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numStateVars = 1;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numDBProperties = 6;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numDBStateVars = 6;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numPropsQuadPt = 6;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numVarsQuadPt = 6;
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_lengthScale = 1.00000000e+03;
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_timeScale = 1.00000000e+00;
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_pressureScale = 2.25000000e+10;
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_densityScale = 1.00000000e+03;
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_dtStableImplicit = 1.00000000e+10;
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::DruckerPrager3DTimeDepData::_numStateVarValues[] = {
+6,
+};
+
+const char* pylith::materials::DruckerPrager3DTimeDepData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"friction_angle",
+"cohesion",
+"dilatation_angle",
+};
+
+const char* pylith::materials::DruckerPrager3DTimeDepData::_dbStateVarValues[] = {
+"plastic-strain-xx",
+"plastic-strain-yy",
+"plastic-strain-zz",
+"plastic-strain-xy",
+"plastic-strain-yz",
+"plastic-strain-xz",
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 5.23598776e-01,
+ 3.00000000e+05,
+ 3.49065850e-01,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 4.36332313e-01,
+ 1.00000000e+04,
+ 4.36332313e-01,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_dbStateVars[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 2.30940108e-01,
+ 3.60000000e+05,
+ 1.48583084e-01,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 1.89338478e-01,
+ 1.21811303e+04,
+ 1.89338478e-01,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_stateVars[] = {
+ 4.10000000e-05,
+ 4.20000000e-05,
+ 4.30000000e-05,
+ 4.40000000e-05,
+ 4.50000000e-05,
+ 4.60000000e-05,
+ 1.10000000e-05,
+ 1.20000000e-05,
+ 1.30000000e-05,
+ 1.40000000e-05,
+ 1.50000000e-05,
+ 1.60000000e-05,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 2.30940108e-01,
+ 1.60000000e-05,
+ 1.48583084e-01,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 1.89338478e-01,
+ 5.41383567e-07,
+ 1.89338478e-01,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_stateVarsNondim[] = {
+ 4.10000000e-05,
+ 4.20000000e-05,
+ 4.30000000e-05,
+ 4.40000000e-05,
+ 4.50000000e-05,
+ 4.60000000e-05,
+ 1.10000000e-05,
+ 1.20000000e-05,
+ 1.30000000e-05,
+ 1.40000000e-05,
+ 1.50000000e-05,
+ 1.60000000e-05,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_strain[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.30000000e-04,
+ 1.40000000e-04,
+ 1.50000000e-04,
+ 1.60000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.30000000e-04,
+ 4.40000000e-04,
+ 4.50000000e-04,
+ 4.60000000e-04,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_stress[] = {
+ 5.55066469e+05,
+ 5.54377286e+05,
+ 5.53688104e+05,
+ -4.36889895e+03,
+ -5.05808169e+03,
+ -5.74726442e+03,
+ 9.73972396e+05,
+ 9.67659189e+05,
+ 9.61345981e+05,
+ -2.85301193e+05,
+ -2.93159178e+05,
+ -3.01017164e+05,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_elasticConsts[] = {
+ 2.69552119e+10,
+ 2.61171422e+10,
+ 2.52103070e+10,
+ -1.14973070e+10,
+ -1.33109780e+10,
+ -1.51246488e+10,
+ 2.56403618e+10,
+ 2.47114216e+10,
+ 2.39200125e+10,
+ -1.09057219e+10,
+ -1.26260719e+10,
+ -1.43464216e+10,
+ 2.42567462e+10,
+ 2.34432321e+10,
+ 2.25609525e+10,
+ -1.03141367e+10,
+ -1.19411658e+10,
+ -1.35681945e+10,
+ -8.77107957e+09,
+ -8.47528712e+09,
+ -8.17949468e+09,
+ 3.68143839e+09,
+ 4.34178877e+09,
+ 4.93337366e+09,
+ -1.01546951e+10,
+ -9.81224227e+09,
+ -9.46978915e+09,
+ 4.34178907e+09,
+ 4.95792951e+09,
+ 5.71160063e+09,
+ -1.15383109e+10,
+ -1.11491971e+10,
+ -1.07600836e+10,
+ 4.93337427e+09,
+ 5.71160093e+09,
+ 6.42106241e+09,
+ 2.65083753e+09,
+ 3.40097380e+09,
+ 3.37872031e+09,
+ -2.01128150e+09,
+ -2.06667779e+09,
+ -2.12207396e+09,
+ 3.40097345e+09,
+ 2.60654697e+09,
+ 3.35689774e+09,
+ -1.99182658e+09,
+ -2.04668706e+09,
+ -2.10154743e+09,
+ 3.37872049e+09,
+ 3.35689820e+09,
+ 2.56268622e+09,
+ -1.97237171e+09,
+ -2.02669634e+09,
+ -2.08102085e+09,
+ -1.00564075e+09,
+ -9.95913200e+08,
+ -9.86186147e+08,
+ 1.06799416e+08,
+ 9.03403590e+08,
+ 9.27618763e+08,
+ -1.03333886e+09,
+ -1.02334327e+09,
+ -1.01334814e+09,
+ 9.03403910e+08,
+ 1.55897025e+08,
+ 9.53168172e+08,
+ -1.06103672e+09,
+ -1.05077369e+09,
+ -1.04051051e+09,
+ 9.27618908e+08,
+ 9.53168172e+08,
+ 2.06328375e+08,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.30000000e+04,
+ 2.40000000e+04,
+ 2.50000000e+04,
+ 2.60000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.30000000e+04,
+ 5.40000000e+04,
+ 5.50000000e+04,
+ 5.60000000e+04,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_initialStrain[] = {
+ 3.60000000e-05,
+ 3.50000000e-05,
+ 3.40000000e-05,
+ 3.30000000e-05,
+ 3.20000000e-05,
+ 3.10000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.30000000e-05,
+ 6.60000000e-05,
+ 6.50000000e-05,
+ 6.40000000e-05,
+};
+
+const double pylith::materials::DruckerPrager3DTimeDepData::_stateVarsUpdated[] = {
+ 6.92302201e-05,
+ 8.02677575e-05,
+ 9.13052948e-05,
+ 1.07630420e-04,
+ 1.18667957e-04,
+ 1.29705495e-04,
+ 2.84142902e-04,
+ 2.94412556e-04,
+ 3.04682210e-04,
+ 4.32906457e-04,
+ 4.45444302e-04,
+ 4.57982146e-04,
+};
+
+pylith::materials::DruckerPrager3DTimeDepData::DruckerPrager3DTimeDepData(void)
+{ // constructor
+ dimension = _dimension;
+ numLocs = _numLocs;
+ numProperties = _numProperties;
+ numStateVars = _numStateVars;
+ numDBProperties = _numDBProperties;
+ numDBStateVars = _numDBStateVars;
+ numPropsQuadPt = _numPropsQuadPt;
+ numVarsQuadPt = _numVarsQuadPt;
+ lengthScale = _lengthScale;
+ timeScale = _timeScale;
+ pressureScale = _pressureScale;
+ densityScale = _densityScale;
+ dtStableImplicit = _dtStableImplicit;
+ numPropertyValues = const_cast<int*>(_numPropertyValues);
+ numStateVarValues = const_cast<int*>(_numStateVarValues);
+ dbPropertyValues = const_cast<char**>(_dbPropertyValues);
+ dbStateVarValues = const_cast<char**>(_dbStateVarValues);
+ dbProperties = const_cast<double*>(_dbProperties);
+ dbStateVars = const_cast<double*>(_dbStateVars);
+ properties = const_cast<double*>(_properties);
+ stateVars = const_cast<double*>(_stateVars);
+ propertiesNondim = const_cast<double*>(_propertiesNondim);
+ stateVarsNondim = const_cast<double*>(_stateVarsNondim);
+ density = const_cast<double*>(_density);
+ strain = const_cast<double*>(_strain);
+ stress = const_cast<double*>(_stress);
+ elasticConsts = const_cast<double*>(_elasticConsts);
+ initialStress = const_cast<double*>(_initialStress);
+ initialStrain = const_cast<double*>(_initialStrain);
+ stateVarsUpdated = const_cast<double*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::DruckerPrager3DTimeDepData::~DruckerPrager3DTimeDepData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.hh 2010-04-18 23:56:03 UTC (rev 16560)
@@ -0,0 +1,104 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerprager3dtimedep.
+
+#if !defined(pylith_materials_druckerprager3dtimedepdata_hh)
+#define pylith_materials_druckerprager3dtimedepdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class DruckerPrager3DTimeDepData;
+ } // pylith
+} // materials
+
+class pylith::materials::DruckerPrager3DTimeDepData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ DruckerPrager3DTimeDepData(void);
+
+ /// Destructor
+ ~DruckerPrager3DTimeDepData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numLocs;
+
+ static const int _numProperties;
+
+ static const int _numStateVars;
+
+ static const int _numDBProperties;
+
+ static const int _numDBStateVars;
+
+ static const int _numPropsQuadPt;
+
+ static const int _numVarsQuadPt;
+
+ static const double _lengthScale;
+
+ static const double _timeScale;
+
+ static const double _pressureScale;
+
+ static const double _densityScale;
+
+ static const double _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const double _dbProperties[];
+
+ static const double _dbStateVars[];
+
+ static const double _properties[];
+
+ static const double _stateVars[];
+
+ static const double _propertiesNondim[];
+
+ static const double _stateVarsNondim[];
+
+ static const double _density[];
+
+ static const double _strain[];
+
+ static const double _stress[];
+
+ static const double _elasticConsts[];
+
+ static const double _initialStress[];
+
+ static const double _initialStrain[];
+
+ static const double _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_druckerprager3dtimedepdata_hh
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh 2010-04-17 22:02:34 UTC (rev 16559)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh 2010-04-18 23:56:03 UTC (rev 16560)
@@ -40,6 +40,11 @@
--data.object=GenMaxwellIsotropic3DElasticData \
--data.parent=ElasticMaterialData
+ python DruckerPrager3DElastic.py \
+ --data.namespace=pylith,materials \
+ --data.object=DruckerPrager3DElasticData \
+ --data.parent=ElasticMaterialData
+
# 2-D ----------------------------------------------------------------
python ElasticPlaneStrain.py \
@@ -94,6 +99,11 @@
--data.object=PowerLaw3DTimeDepData \
--data.parent=ElasticMaterialData
+ python DruckerPrager3DTimeDep.py \
+ --data.namespace=pylith,materials \
+ --data.object=DruckerPrager3DTimeDepData \
+ --data.parent=ElasticMaterialData
+
fi
More information about the CIG-COMMITS
mailing list