[cig-commits] r18416 - in short/3D/PyLith/branches/multifields: libsrc/pylith/materials libsrc/pylith/topology unittests/libtests/materials/data
brad at geodynamics.org
brad at geodynamics.org
Sat May 21 12:47:53 PDT 2011
Author: brad
Date: 2011-05-21 12:47:53 -0700 (Sat, 21 May 2011)
New Revision: 18416
Modified:
short/3D/PyLith/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.cc
short/3D/PyLith/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.hh
short/3D/PyLith/branches/multifields/libsrc/pylith/topology/RefineVol8Face4Edges2.cc
short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc
short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc
Log:
Merge from trunk.
Modified: short/3D/PyLith/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.cc 2011-05-21 19:47:16 UTC (rev 18415)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.cc 2011-05-21 19:47:53 UTC (rev 18416)
@@ -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/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.hh
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.hh 2011-05-21 19:47:16 UTC (rev 18415)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/materials/MaxwellPlaneStrain.hh 2011-05-21 19:47:53 UTC (rev 18416)
@@ -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/branches/multifields/libsrc/pylith/topology/RefineVol8Face4Edges2.cc
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/topology/RefineVol8Face4Edges2.cc 2011-05-21 19:47:16 UTC (rev 18415)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/topology/RefineVol8Face4Edges2.cc 2011-05-21 19:47:53 UTC (rev 18416)
@@ -539,10 +539,12 @@
} // for
} // for
+#if 0 // debuggin
oldSendOverlap->view("OLD SEND OVERLAP");
oldRecvOverlap->view("OLD RECV OVERLAP");
newSendOverlap->view("NEW SEND OVERLAP");
newRecvOverlap->view("NEW RECV OVERLAP");
+#endif
} // overlapAddNewVertces
Modified: short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py 2011-05-21 19:47:16 UTC (rev 18415)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py 2011-05-21 19:47:53 UTC (rev 18416)
@@ -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/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc 2011-05-21 19:47:16 UTC (rev 18415)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc 2011-05-21 19:47:53 UTC (rev 18416)
@@ -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/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py 2011-05-21 19:47:16 UTC (rev 18415)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py 2011-05-21 19:47:53 UTC (rev 18416)
@@ -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/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc 2011-05-21 19:47:16 UTC (rev 18415)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDepData.cc 2011-05-21 19:47:53 UTC (rev 18416)
@@ -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