[cig-commits] r18413 - in short/3D/PyLith/trunk: libsrc/pylith/materials unittests/libtests/materials/data

willic3 at geodynamics.org willic3 at geodynamics.org
Sat May 21 10:00:32 PDT 2011


Author: willic3
Date: 2011-05-21 10:00:31 -0700 (Sat, 21 May 2011)
New Revision: 18413

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.hh
   short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
   short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
   short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc
Log:
Update MaxwellPlaneStrain to use initial stress in the out-of-plane
direction as a state variable, and updated unit tests to test this.



Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc	2011-05-21 16:24:29 UTC (rev 18412)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc	2011-05-21 17:00:31 UTC (rev 18413)
@@ -24,6 +24,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES double_array
+#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -64,17 +65,19 @@
       const char* dbProperties[] = {"density", "vs", "vp" , "viscosity"};
 
       /// Number of state variables.
-      const int numStateVars = 2;
+      const int numStateVars = 3;
 
       /// State variables.
       const Metadata::ParamDescription stateVars[] = {
+	{ "stress_zz_initial", 1, pylith::topology::FieldBase::SCALAR },
 	{ "total_strain", tensorSize, pylith::topology::FieldBase::TENSOR },
-	{ "viscous_strain", 4, pylith::topology::FieldBase::TENSOR },
+	{ "viscous_strain", 4, pylith::topology::FieldBase::OTHER },
       };
 
       /// Values expected in state variables spatial database
-      const int numDBStateVars = tensorSize + 4;
-      const char* dbStateVars[] = { "total-strain-xx",
+      const int numDBStateVars = 1 + tensorSize + 4;
+      const char* dbStateVars[] = { "stress-zz-initial",
+				    "total-strain-xx",
 				    "total-strain-yy",
 				    "total-strain-xy",
 				    "viscous-strain-xx",
@@ -111,14 +114,20 @@
   pylith::materials::MaxwellPlaneStrain::db_vp + 1;
 
 // Indices of state variables
-const int pylith::materials::MaxwellPlaneStrain::s_totalStrain = 0;
+const int pylith::materials::MaxwellPlaneStrain::s_stressZZInitial = 0;
 
+const int pylith::materials::MaxwellPlaneStrain::s_totalStrain =
+  pylith::materials::MaxwellPlaneStrain::s_stressZZInitial + 1;
+
 const int pylith::materials::MaxwellPlaneStrain::s_viscousStrain = 
   pylith::materials::MaxwellPlaneStrain::s_totalStrain + 3;
 
 // Indices of database values (order must match dbStateVars)
-const int pylith::materials::MaxwellPlaneStrain::db_totalStrain = 0;
+const int pylith::materials::MaxwellPlaneStrain::db_stressZZInitial = 0;
 
+const int pylith::materials::MaxwellPlaneStrain::db_totalStrain =
+  pylith::materials::MaxwellPlaneStrain::db_stressZZInitial + 1;
+
 const int pylith::materials::MaxwellPlaneStrain::db_viscousStrain = 
   pylith::materials::MaxwellPlaneStrain::db_totalStrain + 3;
 
@@ -200,8 +209,8 @@
     throw std::runtime_error(msg.str());
   } // if
 
-  const double mu = density * vs*vs;
-  const double lambda = density * vp*vp - 2.0*mu;
+  const double mu = density * vs * vs;
+  const double lambda = density * vp * vp - 2.0 * mu;
 
   if (lambda <= 0.0) {
     std::ostringstream msg;
@@ -285,7 +294,7 @@
   const int numDBValues = dbValues.size();
   assert(_MaxwellPlaneStrain::numDBStateVars == numDBValues);
 
-  const int totalSize = _tensorSize + 4;
+  const int totalSize = 1 + _tensorSize + 4;
   assert(totalSize == _numVarsQuadPt);
   assert(totalSize == numDBValues);
   memcpy(stateValues, &dbValues[0], totalSize*sizeof(double));
@@ -294,6 +303,38 @@
 } // _dbToStateVars
 
 // ----------------------------------------------------------------------
+// Nondimensionalize state variables.
+void
+pylith::materials::MaxwellPlaneStrain::_nondimStateVars(double* const values,
+							const int nvalues) const
+{ // _nondimStateVars
+  assert(0 != _normalizer);
+  assert(0 != values);
+  assert(nvalues == _numVarsQuadPt);
+
+  const double pressureScale = _normalizer->pressureScale();
+  _normalizer->nondimensionalize(&values[s_stressZZInitial], 1, pressureScale);
+
+  PetscLogFlops(1);
+} // _nondimStateVars
+
+// ----------------------------------------------------------------------
+// Dimensionalize state variables.
+void
+pylith::materials::MaxwellPlaneStrain::_dimStateVars(double* const values,
+						     const int nvalues) const
+{ // _dimStateVars
+  assert(0 != _normalizer);
+  assert(0 != values);
+  assert(nvalues == _numVarsQuadPt);
+
+  const double pressureScale = _normalizer->pressureScale();
+  _normalizer->dimensionalize(&values[s_stressZZInitial], 1, pressureScale);
+
+  PetscLogFlops(1);
+} // _dimStateVars
+
+// ----------------------------------------------------------------------
 // Compute density at location from properties.
 void
 pylith::materials::MaxwellPlaneStrain::_calcDensity(double* const density,
@@ -377,13 +418,6 @@
 					 const int initialStrainSize,
 					 const bool computeStateVars)
 { // _calcStressViscoelastic
-  // Note that there is a problem with the way things are presently set up.
-  // Since the stress tensor is required to have a dimension of 3 for a 2D
-  // problem, the ZZ component of stress is never computed, even though it is
-  // part of the solution.  The only way around this at present is to compute
-  // this component as part of postprocessing, which will require knowledge of
-  // the current stress, strain, and viscous strain values as well as the
-  // initial values.
   assert(0 != stress);
   assert(_MaxwellPlaneStrain::tensorSize == stressSize);
   assert(0 != properties);
@@ -402,16 +436,15 @@
   const double mu = properties[p_mu];
   const double lambda = properties[p_lambda];
   const double maxwellTime = properties[p_maxwellTime];
+  const double stressZZInitial = stateVars[s_stressZZInitial];
 
   const double mu2 = 2.0 * mu;
   const double bulkModulus = lambda + mu2 / 3.0;
 
   // Initial stress and strain values
   const double meanStrainInitial = (initialStrain[0] + initialStrain[1]) / 3.0;
-  // :BUG:? CHARLES - Doesn't sigma_zz contribute to
-  // meanStressInitial? Isn't sigma_zz nonzero and found from sigma_xx
-  // and sigma_yy?
-  const double meanStressInitial = (initialStress[0] + initialStress[1]) / 3.0;
+  const double meanStressInitial = (initialStress[0] + initialStress[1] +
+				    stressZZInitial) / 3.0;
   const double devStrainInitial[] = {initialStrain[0] - meanStrainInitial,
 				     initialStrain[1] - meanStrainInitial,
 				     initialStrain[2]};
@@ -444,7 +477,7 @@
   stress[2] = mu2 * (_viscousStrain[3] - devStrainInitial[2]) +
     devStressInitial[2];
 
-  PetscLogFlops(28);
+  PetscLogFlops(30);
 } // _calcStressViscoelastic
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.hh	2011-05-21 16:24:29 UTC (rev 18412)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.hh	2011-05-21 17:00:31 UTC (rev 18413)
@@ -102,9 +102,21 @@
   void _dbToStateVars(double* const stateValues,
 		      const double_array& dbValues);
 
-  // Note: We do not need to dimensionalize or nondimensionalize state
-  // variables because there are strains, which are dimensionless.
+  /** Nondimensionalize state variables..
+   *
+   * @param values Array of state variables.
+   * @param nvalues Number of values.
+   */
+  void _nondimStateVars(double* const values,
+                        const int nvalues) const;
 
+  /** Dimensionalize state variables.
+   *
+   * @param values Array of state variables.
+   * @param nvalues Number of values.
+   */
+  void _dimStateVars(double* const values,
+                     const int nvalues) const;
 
   /** Compute density from properties.
    *
@@ -487,8 +499,10 @@
   static const int db_vp;
   static const int db_viscosity;
 
+  static const int s_stressZZInitial;
   static const int s_totalStrain;
   static const int s_viscousStrain;
+  static const int db_stressZZInitial;
   static const int db_totalStrain;
   static const int db_viscousStrain;
 

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py	2011-05-21 16:24:29 UTC (rev 18412)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py	2011-05-21 17:00:31 UTC (rev 18413)
@@ -45,8 +45,8 @@
     """
     ElasticMaterialApp.__init__(self, name)
 
-    # import pdb
-    # pdb.set_trace()
+    import pdb
+    pdb.set_trace()
 
     numLocs = 2
 
@@ -56,7 +56,8 @@
     self.dbPropertyValues = ["density", "vs", "vp", "viscosity"]
     self.numPropertyValues = numpy.array([1, 1, 1, 1], dtype=numpy.int32)
 
-    self.dbStateVarValues = ["total-strain-xx",
+    self.dbStateVarValues = ["stress-zz-initial",
+                             "total-strain-xx",
                              "total-strain-yy",
                              "total-strain-xy",
                              "viscous-strain-xx",
@@ -64,7 +65,7 @@
                              "viscous-strain-zz",
                              "viscous-strain-xy"
                              ]
-    self.numStateVarValues = numpy.array([3, 4], dtype=numpy.int32)
+    self.numStateVarValues = numpy.array([1, 3, 4], dtype=numpy.int32)
     
     densityA = 2500.0
     vsA = 3000.0
@@ -101,9 +102,9 @@
                                      dtype=numpy.float64)
 
     # TEMPORARY, need to determine how to use initial state variables
-    self.dbStateVars = numpy.zeros( (numLocs, tensorSize+4),
+    self.dbStateVars = numpy.zeros( (numLocs, tensorSize+5),
                                     dtype=numpy.float64)
-    self.stateVars = numpy.zeros( (numLocs, tensorSize+4),
+    self.stateVars = numpy.zeros( (numLocs, tensorSize+5),
                                   dtype=numpy.float64)
 
     mu0 = self.pressureScale
@@ -134,7 +135,7 @@
     self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
     self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts), \
                                       dtype=numpy.float64)
-    self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + 4), \
+    self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + 5), \
                                          dtype=numpy.float64)
 
     (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
@@ -179,8 +180,10 @@
                                 dtype=numpy.float64)
     viscousStrainVec = numpy.ravel(viscousStrain)
     strainVec = numpy.ravel(strain)
+    stressZZInitial = numpy.zeros(1, dtype=numpy.float64)
     
-    stateVarsUpdated = numpy.concatenate((strainVec, viscousStrainVec))
+    stateVarsUpdated = numpy.concatenate((stressZZInitial, strainVec,
+                                          viscousStrainVec))
     return (elasticConsts, numpy.ravel(stress), numpy.ravel(stateVarsUpdated))
   
 

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc	2011-05-21 16:24:29 UTC (rev 18412)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc	2011-05-21 17:00:31 UTC (rev 18413)
@@ -27,15 +27,15 @@
 
 const int pylith::materials::MaxwellPlaneStrainElasticData::_numProperties = 4;
 
-const int pylith::materials::MaxwellPlaneStrainElasticData::_numStateVars = 2;
+const int pylith::materials::MaxwellPlaneStrainElasticData::_numStateVars = 3;
 
 const int pylith::materials::MaxwellPlaneStrainElasticData::_numDBProperties = 4;
 
-const int pylith::materials::MaxwellPlaneStrainElasticData::_numDBStateVars = 7;
+const int pylith::materials::MaxwellPlaneStrainElasticData::_numDBStateVars = 8;
 
 const int pylith::materials::MaxwellPlaneStrainElasticData::_numPropsQuadPt = 4;
 
-const int pylith::materials::MaxwellPlaneStrainElasticData::_numVarsQuadPt = 7;
+const int pylith::materials::MaxwellPlaneStrainElasticData::_numVarsQuadPt = 8;
 
 const double pylith::materials::MaxwellPlaneStrainElasticData::_lengthScale =   1.00000000e+03;
 
@@ -55,6 +55,7 @@
 };
 
 const int pylith::materials::MaxwellPlaneStrainElasticData::_numStateVarValues[] = {
+1,
 3,
 4,
 };
@@ -67,6 +68,7 @@
 };
 
 const char* pylith::materials::MaxwellPlaneStrainElasticData::_dbStateVarValues[] = {
+"stress-zz-initial",
 "total-strain-xx",
 "total-strain-yy",
 "total-strain-xy",
@@ -102,6 +104,8 @@
   0.00000000e+00,
   0.00000000e+00,
   0.00000000e+00,
+  0.00000000e+00,
+  0.00000000e+00,
 };
 
 const double pylith::materials::MaxwellPlaneStrainElasticData::_properties[] = {
@@ -130,6 +134,8 @@
   0.00000000e+00,
   0.00000000e+00,
   0.00000000e+00,
+  0.00000000e+00,
+  0.00000000e+00,
 };
 
 const double pylith::materials::MaxwellPlaneStrainElasticData::_propertiesNondim[] = {
@@ -158,6 +164,8 @@
   0.00000000e+00,
   0.00000000e+00,
   0.00000000e+00,
+  0.00000000e+00,
+  0.00000000e+00,
 };
 
 const double pylith::materials::MaxwellPlaneStrainElasticData::_density[] = {
@@ -223,6 +231,7 @@
 };
 
 const double pylith::materials::MaxwellPlaneStrainElasticData::_stateVarsUpdated[] = {
+  0.00000000e+00,
   1.10000000e-04,
   1.20000000e-04,
   1.40000000e-04,
@@ -230,6 +239,7 @@
   4.33333333e-05,
  -7.66666667e-05,
   1.40000000e-04,
+  0.00000000e+00,
   4.10000000e-04,
   4.20000000e-04,
   4.40000000e-04,

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py	2011-05-21 16:24:29 UTC (rev 18412)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py	2011-05-21 17:00:31 UTC (rev 18413)
@@ -45,8 +45,8 @@
     """
     ElasticMaterialApp.__init__(self, name)
 
-    # import pdb
-    # pdb.set_trace()
+    import pdb
+    pdb.set_trace()
 
     numLocs = 2
 
@@ -57,7 +57,8 @@
     self.propertyValues = ["density", "mu", "lambda", "maxwellTime"]
     self.numPropertyValues = numpy.array([1, 1, 1, 1], dtype=numpy.int32)
 
-    self.dbStateVarValues = ["total-strain-xx",
+    self.dbStateVarValues = ["stress-zz-initial",
+                             "total-strain-xx",
                              "total-strain-yy",
                              "total-strain-xy",
                              "viscous-strain-xx",
@@ -65,8 +66,8 @@
                              "viscous-strain-zz",
                              "viscous-strain-xy"
                              ]
-    self.stateVarValues = ["total-strain", "viscous-strain"]
-    self.numStateVarValues = numpy.array([3, 4], dtype=numpy.int32)
+    self.stateVarValues = ["stress-zz-initial","total-strain", "viscous-strain"]
+    self.numStateVarValues = numpy.array([1, 3, 4], dtype=numpy.int32)
 
     self.dt = 2.0e5
 
@@ -81,6 +82,7 @@
     lambdaA = vpA*vpA*densityA - 2.0*muA
     maxwellTimeA = viscosityA / muA
     meanStrainA = (strainA[0] + strainA[1])/3.0
+    stressInitialZZA = numpy.array([2.0e4], dtype=numpy.float64)
 
     densityB = 2000.0
     vsB = 1200.0
@@ -93,6 +95,7 @@
     lambdaB = vpB*vpB*densityB - 2.0*muB
     maxwellTimeB = viscosityB / muB
     meanStrainB = (strainB[0] + strainB[1])/3.0
+    stressInitialZZB = numpy.array([5.0e4], dtype=numpy.float64)
 
     diag = numpy.array([1.0, 1.0, 1.0, 0.0],
                        dtype=numpy.float64)
@@ -113,8 +116,11 @@
                                      dtype=numpy.float64)
 
     # TEMPORARY, need to determine how to use initial state variables
-    self.dbStateVars = numpy.zeros( (numLocs, tensorSize+4),
+    # At present, only the first (stressInitialZZ) is being used.
+    self.dbStateVars = numpy.zeros( (numLocs, tensorSize+5),
                                     dtype=numpy.float64)
+    self.dbStateVars[0, 0] = stressInitialZZA
+    self.dbStateVars[1, 0] = stressInitialZZB
 
     mu0 = self.pressureScale
     density0 = self.densityScale
@@ -156,17 +162,26 @@
     viscousStrainVecB = numpy.ravel(viscousStrainB)
     strainVecA = numpy.array(strainA, dtype=numpy.float64)
     strainVecB = numpy.array(strainB, dtype=numpy.float64)
-    stateVarsA = numpy.concatenate((strainVecA, viscousStrainVecA))
-    stateVarsB = numpy.concatenate((strainVecB, viscousStrainVecB))
-    self.stateVars = numpy.array([ [stateVarsA],
-                                   [stateVarsB] ],
+    stateVarsA = numpy.concatenate((stressInitialZZA, strainVecA,
+                                    viscousStrainVecA))
+    stateVarsB = numpy.concatenate((stressInitialZZB, strainVecB,
+                                    viscousStrainVecB))
+    self.stateVars = numpy.array([ stateVarsA,
+                                   stateVarsB ],
                                  dtype=numpy.float64)
-    self.stateVarsNondim = self.stateVars # no scaling
+    stressInitialZZANondim = stressInitialZZA/mu0
+    stressInitialZZBNondim = stressInitialZZB/mu0
+    self.stateVarsNondim = numpy.zeros( (numLocs, tensorSize + 5),
+                                        dtype=numpy.float64)
+    self.stateVarsNondim[:] = self.stateVars[:] # no scaling
+    # Scale stressInitialZZ
+    self.stateVarsNondim[0, 0] = stressInitialZZANondim
+    self.stateVarsNondim[1, 0] = stressInitialZZBNondim
     
     self.strain = numpy.array([totalStrainA, totalStrainB],
                                dtype=numpy.float64)
     self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
-    self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + 4),
+    self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + 5),
                                          dtype=numpy.float64)
     self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts),
                                       dtype=numpy.float64)
@@ -175,20 +190,26 @@
                      self._calcStress(strainA, 
                                       muA, lambdaA, maxwellTimeA,
                                       totalStrainA, viscousStrainA,
-                                      initialStressA, initialStrainA)
+                                      initialStressA, initialStrainA,
+                                      stateVarsA)
     (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
                               self._calcStress(strainB, 
                                                muB, lambdaB, maxwellTimeB, 
                                                totalStrainB, viscousStrainB,
-                                               initialStressB, initialStrainB)
+                                               initialStressB, initialStrainB,
+                                               stateVarsB)
 
+    self.stateVarsUpdated[0, 0] = stressInitialZZA
+    self.stateVarsUpdated[1, 0] = stressInitialZZB
+
     self.dtStableImplicit = 0.2*min(maxwellTimeA, maxwellTimeB)
     return
 
 
   def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
-                         muV, lambdaV, maxwellTimeV, dqV, totalStrainT,
-                         viscousStrainT, initialStress, initialStrain):
+                           muV, lambdaV, maxwellTimeV, dqV, totalStrainT,
+                           viscousStrainT, initialStress, initialStrain,
+                           stateVars):
     """
     Function to compute a particular stress component as a function of a
     strain component.
@@ -203,13 +224,14 @@
                                                         totalStrainT,
                                                         viscousStrainT,
                                                         initialStress,
-                                                        initialStrain)
+                                                        initialStrain,
+                                                        stateVars)
     return stressTpdt[stressComp]
 
 
   def _computeStress(self, strainTpdt, muV, lambdaV, maxwellTimeV, dqV,
                      strainT, viscousStrainT,
-                     initialStress, initialStrain):
+                     initialStress, initialStrain, stateVars):
     """
     Function to compute stresses and viscous strains for a given strain.
     """
@@ -222,7 +244,7 @@
     meanStrainInitial = \
     (initialStrain[0] + initialStrain[1]) / 3.0
     meanStressInitial = \
-    (initialStress[0] + initialStress[1]) / 3.0
+    (initialStress[0] + initialStress[1] + stateVars[0]) / 3.0
 
     devStrainInitial = initialStrain - numpy.array(diag) * meanStrainInitial
     devStressInitial = initialStress - numpy.array(diag) * meanStressInitial
@@ -263,9 +285,7 @@
     devStressTpdt22 = elasFac * \
                       (viscousStrainTpdt[1] - devStrainInitial[1]) + \
                       devStressInitial[1]
-    # This is somewhat messed up, since we can't specify initial stresses or
-    # strains in the out-of-plane direction.
-    devStressTpdt33 = elasFac * viscousStrainTpdt[2]
+    devStressTpdt33 = elasFac * viscousStrainTpdt[2] + stateVars[0]
     devStressTpdt12 = elasFac * \
                       (viscousStrainTpdt[3] - devStrainInitial[2]) + \
                       devStressInitial[2]
@@ -303,7 +323,7 @@
 
   
   def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, totalStrainV,
-                  viscousStrainV, initialStressV, initialStrainV):
+                  viscousStrainV, initialStressV, initialStrainV, stateVars):
     """
     Compute stress, derivative of elasticity matrix, and updated state
     variables. This assumes behavior is always viscoelastic.
@@ -316,6 +336,7 @@
     viscousStrainT = numpy.array(viscousStrainV, dtype=numpy.float64)
     initialStress = numpy.array(initialStressV, dtype=numpy.float64)
     initialStrain = numpy.array(initialStrainV, dtype=numpy.float64)
+    stressZZInitial = numpy.array([stateVars[0]], dtype=numpy.float64)
     
     # Compute current stress and viscous strain.
     dqV = self._computeViscousFactor(maxwellTimeV)
@@ -327,12 +348,14 @@
                                                         strainT,
                                                         viscousStrainT,
                                                         initialStress,
-                                                        initialStrain)
+                                                        initialStrain,
+                                                        stateVars)
                                                         
     # Form updated state variables
     strainTpdtVec = numpy.ravel(strainTpdt)
     viscousStrainTpdtVec = numpy.ravel(viscousStrainTpdt)
-    stateVarsUpdated = numpy.concatenate((strainTpdtVec, viscousStrainTpdtVec))
+    stateVarsUpdated = numpy.concatenate((stressZZInitial, strainTpdtVec,
+                                          viscousStrainTpdtVec))
 
     # Compute components of tangent constitutive matrix using numerical
     # derivatives.
@@ -356,7 +379,8 @@
                                                      strainT,
                                                      viscousStrainT,
                                                      initialStress,
-                                                     initialStrain),
+                                                     initialStrain,
+                                                     stateVars),
                                                order=derivOrder)
         elasticConstsList.append(dStressDStrain)
 

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc	2011-05-21 16:24:29 UTC (rev 18412)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc	2011-05-21 17:00:31 UTC (rev 18413)
@@ -27,15 +27,15 @@
 
 const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numProperties = 4;
 
-const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numStateVars = 2;
+const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numStateVars = 3;
 
 const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numDBProperties = 4;
 
-const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numDBStateVars = 7;
+const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numDBStateVars = 8;
 
 const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numPropsQuadPt = 4;
 
-const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numVarsQuadPt = 7;
+const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numVarsQuadPt = 8;
 
 const double pylith::materials::MaxwellPlaneStrainTimeDepData::_lengthScale =   1.00000000e+03;
 
@@ -55,6 +55,7 @@
 };
 
 const int pylith::materials::MaxwellPlaneStrainTimeDepData::_numStateVarValues[] = {
+1,
 3,
 4,
 };
@@ -67,6 +68,7 @@
 };
 
 const char* pylith::materials::MaxwellPlaneStrainTimeDepData::_dbStateVarValues[] = {
+"stress-zz-initial",
 "total-strain-xx",
 "total-strain-yy",
 "total-strain-xy",
@@ -88,6 +90,7 @@
 };
 
 const double pylith::materials::MaxwellPlaneStrainTimeDepData::_dbStateVars[] = {
+  2.00000000e+04,
   0.00000000e+00,
   0.00000000e+00,
   0.00000000e+00,
@@ -95,6 +98,7 @@
   0.00000000e+00,
   0.00000000e+00,
   0.00000000e+00,
+  5.00000000e+04,
   0.00000000e+00,
   0.00000000e+00,
   0.00000000e+00,
@@ -116,6 +120,7 @@
 };
 
 const double pylith::materials::MaxwellPlaneStrainTimeDepData::_stateVars[] = {
+  2.00000000e+04,
   1.10000000e-04,
   1.20000000e-04,
   1.40000000e-04,
@@ -123,6 +128,7 @@
   4.33333333e-05,
  -7.66666667e-05,
   1.40000000e-04,
+  5.00000000e+04,
   4.10000000e-04,
   4.20000000e-04,
   4.40000000e-04,
@@ -144,6 +150,7 @@
 };
 
 const double pylith::materials::MaxwellPlaneStrainTimeDepData::_stateVarsNondim[] = {
+  8.88888889e-07,
   1.10000000e-04,
   1.20000000e-04,
   1.40000000e-04,
@@ -151,6 +158,7 @@
   4.33333333e-05,
  -7.66666667e-05,
   1.40000000e-04,
+  2.22222222e-06,
   4.10000000e-04,
   4.20000000e-04,
   4.40000000e-04,
@@ -185,13 +193,13 @@
 
 const double pylith::materials::MaxwellPlaneStrainTimeDepData::_elasticConsts[] = {
   6.74326010e+10,
+  2.25336994e+10,
+  0.00000000e+00,
   2.25336989e+10,
+  6.74326010e+10,
   0.00000000e+00,
-  2.25336994e+10,
-  6.74326005e+10,
   0.00000000e+00,
   0.00000000e+00,
-  0.00000000e+00,
   4.48989011e+10,
   8.63988977e+09,
   2.88005546e+09,
@@ -223,6 +231,7 @@
 };
 
 const double pylith::materials::MaxwellPlaneStrainTimeDepData::_stateVarsUpdated[] = {
+  2.00000000e+04,
   1.20000000e-04,
   1.30000000e-04,
   1.50000000e-04,
@@ -230,6 +239,7 @@
   4.64646160e-05,
  -8.29741309e-05,
   1.49348949e-04,
+  5.00000000e+04,
   4.20000000e-04,
   4.30000000e-04,
   4.50000000e-04,



More information about the CIG-COMMITS mailing list