[cig-commits] r18497 - in short/3D/PyLith/trunk: libsrc/pylith/materials unittests/libtests/materials unittests/libtests/materials/data
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon May 30 22:27:02 PDT 2011
Author: willic3
Date: 2011-05-30 22:27:02 -0700 (Mon, 30 May 2011)
New Revision: 18497
Modified:
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc
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/GenMaxwellIsotropic3DElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
Log:
More cleaning up of materials.
Primary changes are use of initial strains and how unit tests are
performed for time-dependent cases (e.g. numerical calculation of
constitutive matrix, etc.).
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -23,6 +23,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
@@ -375,7 +376,7 @@
assert(_numVarsQuadPt == numStateVars);
// It's unclear what to do for an elasto-plastic material, which has no
// inherent time scale. For now, just set dtStable to a large value.
- const double dtStable = 1.0e10;
+ const double dtStable = pylith::PYLITH_MAXDOUBLE;
PetscLogFlops(0);
return dtStable;
} // _stableTimeStepImplicit
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -223,7 +223,6 @@
assert(0 != initialStrain);
assert(_ElasticIsotropic3D::tensorSize == initialStrainSize);
- const double density = properties[p_density];
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
@@ -277,7 +276,6 @@
assert(0 != initialStrain);
assert(_ElasticIsotropic3D::tensorSize == initialStrainSize);
- const double density = properties[p_density];
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -77,7 +77,7 @@
};
/// Number of state variables.
- const int numStateVars = 1+numMaxwellModels;
+ const int numStateVars = 1 + numMaxwellModels;
/// State variables. :KLUDGE: Not generalized over number of models.
const Metadata::ParamDescription stateVars[numStateVars] = {
@@ -88,7 +88,7 @@
};
// Values expected in state variables spatial database
- const int numDBStateVars = tensorSize + numMaxwellModels*tensorSize;
+ const int numDBStateVars = tensorSize + numMaxwellModels * tensorSize;
const char* dbStateVars[numDBStateVars] = {"total-strain-xx",
"total-strain-yy",
"total-strain-zz",
@@ -310,7 +310,7 @@
propValues[p_maxwellTime + imodel] = maxwellTime;
} // for
- PetscLogFlops(6+2*numMaxwellModels);
+ PetscLogFlops(6 + 3 * numMaxwellModels);
} // _dbToProperties
// ----------------------------------------------------------------------
@@ -505,25 +505,41 @@
const double mu2 = 2.0 * mu;
const double bulkModulus = lambda + mu2/3.0;
- // :TODO: Need to determine how to incorporate initial strain and
- // state variables
- const double e11 = totalStrain[0] - initialStrain[0];
- const double e22 = totalStrain[1] - initialStrain[1];
- const double e33 = totalStrain[2] - initialStrain[2];
+ // Initial stress and strain values
+ const double meanStrainInitial =
+ (initialStrain[0] + initialStrain[1] + initialStrain[2])/3.0;
+ const double meanStressInitial =
+ (initialStress[0] + initialStress[1] + initialStress[2])/3.0;
+ const double devStrainInitial[] = {initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ initialStrain[2] - meanStrainInitial,
+ initialStrain[3],
+ initialStrain[4],
+ initialStrain[5]};
+ const double devStressInitial[] = {initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ initialStress[2] - meanStressInitial,
+ initialStress[3],
+ initialStress[4],
+ initialStress[5]};
+ // :TODO: Need to determine how to incorporate state variables
+ // Mean stress and strain for time t + dt
+ const double meanStrainTpdt = (totalStrain[0] +
+ totalStrain[1] +
+ totalStrain[2]) / 3.0;
- const double e123 = e11 + e22 + e33;
- const double meanStrainTpdt = e123 / 3.0;
- const double meanStressTpdt = bulkModulus * e123;
+ const double meanStressTpdt = 3.0 * bulkModulus *
+ (meanStrainTpdt - meanStrainInitial) + meanStressInitial;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
-
+
double visFrac = 0.0;
for (int imodel=0; imodel < numMaxwellModels; ++imodel)
visFrac += muRatio[imodel];
assert(visFrac <= 1.0);
const double elasFrac = 1.0 - visFrac;
- PetscLogFlops(8 + numMaxwellModels);
+ PetscLogFlops(23 + numMaxwellModels);
// Get viscous strains
if (computeStateVars) {
@@ -546,18 +562,20 @@
double devStrainTpdt = 0.0;
double devStressTpdt = 0.0;
for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
- devStressTpdt = elasFrac*devStrainTpdt;
+ devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt -
+ devStrainInitial[iComp];
+ devStressTpdt = elasFrac * devStrainTpdt;
for (int model=0; model < numMaxwellModels; ++model) {
- devStressTpdt += muRatio[model] * _viscousStrain[model*tensorSize+iComp];
+ devStressTpdt += muRatio[model] *
+ _viscousStrain[model * tensorSize+iComp];
} // for
- devStressTpdt = mu2*devStressTpdt;
- stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
- initialStress[iComp];
+ devStressTpdt = mu2 * devStressTpdt + devStressInitial[iComp];
+ stress[iComp] = diag[iComp] * (meanStressTpdt + meanStressInitial) +
+ devStressTpdt;
} // for
- PetscLogFlops((7 + 2*numMaxwellModels) * tensorSize);
+ PetscLogFlops((9 + 3 * numMaxwellModels) * tensorSize);
} // _calcStressViscoelastic
// ----------------------------------------------------------------------
@@ -769,10 +787,14 @@
const int tensorSize = _tensorSize;
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
- const double meanStrainTpdt = (e11 + e22 + e33) / 3.0;
+ const double strainTpdt[] = {totalStrain[0] - initialStrain[0],
+ totalStrain[1] - initialStrain[1],
+ totalStrain[2] - initialStrain[2],
+ totalStrain[3] - initialStrain[3],
+ totalStrain[4] - initialStrain[4],
+ totalStrain[5] - initialStrain[5]};
+ const double meanStrainTpdt =
+ (strainTpdt[0] + strainTpdt[1] + strainTpdt[2])/3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
@@ -784,7 +806,7 @@
double devStrain = 0.0;
double shearRatio = 0.0;
for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStrain = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ devStrain = strainTpdt[iComp] - diag[iComp] * meanStrainTpdt;
// Maxwell model 1
stateVars[s_viscousStrain1+iComp] = devStrain;
// Maxwell model 2
@@ -792,7 +814,7 @@
// Maxwell model 3
stateVars[s_viscousStrain3+iComp] = devStrain;
} // for
- PetscLogFlops(3 + 2 * tensorSize);
+ PetscLogFlops(6 + 2 * tensorSize);
_needNewJacobian = true;
} // _updateStateVarsElastic
@@ -916,16 +938,8 @@
};
// :TODO: Need to account for initial values for state variables
- // and the initial strain??
-
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
- const double e12 = totalStrain[3];
- const double e23 = totalStrain[4];
- const double e13 = totalStrain[5];
-
- const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
+ const double meanStrainTpdt =
+ (totalStrain[0] + totalStrain[1] + totalStrain[2])/3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
const double meanStrainT =
@@ -948,14 +962,15 @@
double deltaStrain = 0.0;
for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
- devStrainT = stateVars[s_totalStrain+iComp] - diag[iComp]*meanStrainT;
+ devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ devStrainT = stateVars[s_totalStrain+iComp] - diag[iComp] * meanStrainT;
deltaStrain = devStrainTpdt - devStrainT;
// Maxwell model 1
int imodel = 0;
if (0.0 != muRatio[imodel]) {
- _viscousStrain[imodel*tensorSize+iComp] = exp(-_dt/maxwellTime[imodel])*
+ _viscousStrain[imodel * tensorSize+iComp] =
+ exp(-_dt/maxwellTime[imodel]) *
stateVars[s_viscousStrain1 + iComp] + dq[imodel] * deltaStrain;
PetscLogFlops(6);
} // if
@@ -963,7 +978,8 @@
// Maxwell model 2
imodel = 1;
if (0.0 != muRatio[imodel]) {
- _viscousStrain[imodel*tensorSize+iComp] = exp(-_dt/maxwellTime[imodel])*
+ _viscousStrain[imodel*tensorSize+iComp] =
+ exp(-_dt/maxwellTime[imodel]) *
stateVars[s_viscousStrain2 + iComp] + dq[imodel] * deltaStrain;
PetscLogFlops(6);
} // if
@@ -971,14 +987,15 @@
// Maxwell model 3
imodel = 2;
if (0.0 != muRatio[imodel]) {
- _viscousStrain[imodel*tensorSize+iComp] = exp(-_dt/maxwellTime[imodel])*
+ _viscousStrain[imodel*tensorSize+iComp] =
+ exp(-_dt/maxwellTime[imodel]) *
stateVars[s_viscousStrain3 + iComp] + dq[imodel] * deltaStrain;
PetscLogFlops(6);
} // if
} // for
- PetscLogFlops(5*tensorSize);
+ PetscLogFlops(5 * tensorSize);
} // _computeStateVars
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -305,7 +305,7 @@
propValues[p_maxwellTime + imodel] = maxwellTime;
} // for
- PetscLogFlops(6 + 2 * numMaxwellModels);
+ PetscLogFlops(6 + 3 * numMaxwellModels);
} // _dbToProperties
// ----------------------------------------------------------------------
@@ -944,7 +944,7 @@
} // for
- PetscLogFlops(numMaxwellModels * (4 * (5 + 18)));
+ PetscLogFlops(5 * 4);
} // _computeStateVars
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellIsotropic3D.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellIsotropic3D.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -432,8 +432,7 @@
initialStress[4],
initialStress[5]};
- // :TODO: Need to determine how to incorporate initial strain and
- // state variables
+ // :TODO: Need to determine how to incorporate state variables
const double meanStrainTpdt = (totalStrain[0] +
totalStrain[1] +
totalStrain[2]) / 3.0;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -128,69 +128,10 @@
void
pylith::materials::TestGenMaxwellIsotropic3D::test_updateStateVarsElastic(void)
{ // testUpdateStateVarsElastic
- // :TODO: Use TestElasticMaterial::test_updateStateVars
- // instead. This requires moving the calculation of the expected
- // state vars below to the Python code (where it belongs) and
- // setting the stateVarsUpdate attribute in the Python object.
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
- GenMaxwellIsotropic3D material;
- material.useElasticBehavior(true);
-
- const bool computeStateVars = true;
-
- const int numLocs = _dataElastic->numLocs;
- const int numPropsQuadPt = _dataElastic->numPropsQuadPt;
- const int numVarsQuadPt = _dataElastic->numVarsQuadPt;
- const int tensorSize = material.tensorSize();
-
- double_array stress(tensorSize);
- double_array properties(numPropsQuadPt);
- double_array stateVars(numVarsQuadPt);
- double_array strain(tensorSize);
- double_array initialStress(tensorSize);
- double_array initialStrain(tensorSize);
-
- for (int iLoc=0; iLoc < numLocs; ++iLoc) {
- memcpy(&properties[0], &_dataElastic->properties[iLoc*numPropsQuadPt],
- numPropsQuadPt*sizeof(double));
- memcpy(&stateVars[0], &_dataElastic->stateVars[iLoc*numVarsQuadPt],
- numVarsQuadPt*sizeof(double));
- memcpy(&strain[0], &_dataElastic->strain[iLoc*tensorSize],
- tensorSize*sizeof(double));
- memcpy(&initialStress[0], &_dataElastic->initialStress[iLoc*tensorSize],
- tensorSize*sizeof(double));
- memcpy(&initialStrain[0], &_dataElastic->initialStrain[iLoc*tensorSize],
- tensorSize*sizeof(double));
-
- const double meanStrain = (strain[0] + strain[1] + strain[2]) / 3.0;
-
- // Compute expected state variables
- double_array stateVarsE(numVarsQuadPt);
- const int s_totalStrain = 0;
- const int s_viscousStrain = s_totalStrain + tensorSize;
-
- // State variable 'total_strain' should match 'strain'
- for (int i=0; i < tensorSize; ++i)
- stateVarsE[s_totalStrain+i] = strain[i];
-
- // State variable 'viscous_strain'
- const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
- const int numMaxwellModels = 3;
- for (int imodel=0; imodel < numMaxwellModels; ++imodel)
- for (int i=0; i < tensorSize; ++i)
- stateVarsE[s_viscousStrain+imodel*tensorSize+i] =
- strain[i] - diag[i]*meanStrain;
-
- material._updateStateVars(&stateVars[0], stateVars.size(),
- &properties[0], properties.size(),
- &strain[0], strain.size(),
- &initialStress[0], initialStress.size(),
- &initialStrain[0], initialStrain.size());
-
- const double tolerance = 1.0e-06;
- for (int i=0; i < numVarsQuadPt; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(stateVarsE[i], stateVars[i], tolerance);
- } // for
+ test_updateStateVars();
} // testUpdateStateVarsElastic
// ----------------------------------------------------------------------
@@ -228,105 +169,15 @@
void
pylith::materials::TestGenMaxwellIsotropic3D::test_updateStateVarsTimeDep(void)
{ // testUpdateStateVarsTimeDep
- // :TODO: Use TestElasticMaterial::test_updateStateVars
- // instead. This requires moving the calculation of the expected
- // state vars below to the Python code (where it belongs) and
- // setting the stateVarsUpdate attribute in the Python object.
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
- GenMaxwellIsotropic3D material;
- material.useElasticBehavior(false);
-
delete _dataElastic; _dataElastic = new GenMaxwellIsotropic3DTimeDepData();
- const double dt = 2.0e+5;
- material.timeStep(dt);
+ double dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_updateStateVars();
- const bool computeStateVars = true;
-
- const int numLocs = _dataElastic->numLocs;
- const int numPropsQuadPt = _dataElastic->numPropsQuadPt;
- const int numVarsQuadPt = _dataElastic->numVarsQuadPt;
- const int tensorSize = material.tensorSize();
-
- double_array stress(tensorSize);
- double_array properties(numPropsQuadPt);
- double_array stateVars(numVarsQuadPt);
- double_array strain(tensorSize);
- double_array initialStress(tensorSize);
- double_array initialStrain(tensorSize);
-
- for (int iLoc=0; iLoc < numLocs; ++iLoc) {
- memcpy(&properties[0], &_dataElastic->properties[iLoc*numPropsQuadPt],
- numPropsQuadPt*sizeof(double));
- memcpy(&stateVars[0], &_dataElastic->stateVars[iLoc*numVarsQuadPt],
- numVarsQuadPt*sizeof(double));
- memcpy(&strain[0], &_dataElastic->strain[iLoc*tensorSize],
- tensorSize*sizeof(double));
- memcpy(&initialStress[0], &_dataElastic->initialStress[iLoc*tensorSize],
- tensorSize*sizeof(double));
- memcpy(&initialStrain[0], &_dataElastic->initialStrain[iLoc*tensorSize],
- tensorSize*sizeof(double));
-
- // Compute expected state variables
- double_array stateVarsE(numVarsQuadPt);
- const int numMaxwellModels = 3;
- const int s_totalStrain = 0;
- const int s_viscousStrain = s_totalStrain + tensorSize;
- const int p_shearRatio = 3;
- const int p_maxwellTime = p_shearRatio + numMaxwellModels;
-
- // State variable 'total_strain' should match 'strain'
- for (int i=0; i < tensorSize; ++i)
- stateVarsE[s_totalStrain+i] = strain[i];
-
- // State variable 'viscous_strain'
- double_array maxwellTime(numMaxwellModels);
- double_array shearRatio(numMaxwellModels);
- double_array dq(numMaxwellModels);
- for (int i=0; i < numMaxwellModels; ++i) {
- shearRatio[i] = properties[p_shearRatio+i];
- maxwellTime[i] = properties[p_maxwellTime+i];
- dq[i] = maxwellTime[i] * (1.0 - exp(-dt/maxwellTime[i]))/dt;
- } // for
- double_array strainT(tensorSize);
- for (int i=0; i < tensorSize; ++i)
- strainT[i] = stateVars[s_totalStrain+i];
- const double meanStrainT =
- (stateVars[s_totalStrain+0] +
- stateVars[s_totalStrain+1] +
- stateVars[s_totalStrain+2]) / 3.0;
- const double meanStrainTpdt = (strain[0] + strain[1] + strain[2]) / 3.0;
-
- double devStrainTpdt = 0.0;
- double devStrainT = 0.0;
- double deltaStrain = 0.0;
- double visStrain = 0.0;
- const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
-
- for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStrainTpdt = strain[iComp] - diag[iComp]*meanStrainTpdt;
- devStrainT = strainT[iComp] - diag[iComp]*meanStrainT;
- deltaStrain = devStrainTpdt - devStrainT;
- for (int imodel=0; imodel < numMaxwellModels; ++imodel) {
- stateVarsE[s_viscousStrain+imodel*tensorSize+iComp] =
- exp(-dt/maxwellTime[imodel]) *
- stateVars[s_viscousStrain+imodel*tensorSize+iComp] + dq[imodel] * deltaStrain;
- } // for
- } // for
-
- material._updateStateVars(&stateVars[0], stateVars.size(),
- &properties[0], properties.size(),
- &strain[0], strain.size(),
- &initialStress[0], initialStress.size(),
- &initialStrain[0], initialStrain.size());
-
- const double tolerance = 1.0e-06;
- for (int i=0; i < numVarsQuadPt; ++i)
- if (stateVarsE[i] > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stateVars[i]/stateVarsE[i], tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(stateVarsE[i], stateVars[i], tolerance);
- } // for
} // testUpdateStateVarsTimeDep
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py 2011-05-31 05:27:02 UTC (rev 18497)
@@ -189,7 +189,7 @@
plasStrainB,
initialStressB, initialStrainB)
- self.dtStableImplicit = 1.0e10
+ self.dtStableImplicit = 1.0e30
return
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -45,7 +45,7 @@
const double pylith::materials::DruckerPrager3DTimeDepData::_densityScale = 1.00000000e+03;
-const double pylith::materials::DruckerPrager3DTimeDepData::_dtStableImplicit = 1.00000000e+10;
+const double pylith::materials::DruckerPrager3DTimeDepData::_dtStableImplicit = 1.00000000e+30;
const int pylith::materials::DruckerPrager3DTimeDepData::_numPropertyValues[] = {
1,
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElastic.py 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElastic.py 2011-05-31 05:27:02 UTC (rev 18497)
@@ -46,15 +46,19 @@
"""
ElasticMaterialApp.__init__(self, name)
+ # import pdb
+ # pdb.set_trace()
+
numLocs = 2
- self.dt = 2.0e5
+
self.dimension = dimension
self.numLocs = numLocs
self.dbPropertyValues = ["density", "vs", "vp",
"shear-ratio-1", "shear-ratio-2", "shear-ratio-3",
"viscosity-1", "viscosity-2", "viscosity-3"]
- self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1, 1, 1, 1],
+ dtype=numpy.int32)
self.dbStateVarValues = ["total-strain-xx",
"total-strain-yy",
@@ -91,7 +95,7 @@
viscosityA = [1.0e18, 1.0e17, 1.0e19]
strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.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]
+ initialStrainA = [3.1e-5, 3.2e-5, 3.3e-5, 3.4e-5, 3.5e-5, 3.6e-5]
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
@@ -102,12 +106,12 @@
viscosityB = [1.0e18, 1.0e19, 1.0e20]
strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-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]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.3e-5, 6.4e-5, 6.5e-5, 6.6e-5]
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
- maxwellTimeA = [0.0, 0.0, 0.0]
- maxwellTimeB = [0.0, 0.0, 0.0]
+ maxwellTimeA = [1.0e30, 1.0e30, 1.0e30]
+ maxwellTimeB = [1.0e30, 1.0e30, 1.0e30]
for i in xrange(numMaxwellModels):
if shearRatioA[i] != 0.0:
maxwellTimeA[i] = viscosityA[i]/(muA*shearRatioA[i])
@@ -127,10 +131,12 @@
self.properties = numpy.array([propA, propB], dtype=numpy.float64)
# TEMPORARY, need to determine how to use initial state variables
- self.dbStateVars = numpy.zeros( (numLocs, tensorSize+numMaxwellModels*tensorSize),
- dtype=numpy.float64)
- self.stateVars = numpy.zeros( (numLocs, tensorSize+numMaxwellModels*tensorSize),
- dtype=numpy.float64)
+ self.dbStateVars = numpy.zeros((numLocs,
+ tensorSize + numMaxwellModels * tensorSize),
+ dtype=numpy.float64)
+ self.stateVars = numpy.zeros((numLocs,
+ tensorSize + numMaxwellModels * tensorSize),
+ dtype=numpy.float64)
mu0 = self.pressureScale
density0 = self.densityScale
@@ -138,10 +144,12 @@
self.propertiesNondim = \
numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0,
shearRatioA[0], shearRatioA[1], shearRatioA[2],
- maxwellTimeA[0]/time0, maxwellTimeA[1]/time0, maxwellTimeA[2]/time0],
+ maxwellTimeA[0]/time0, maxwellTimeA[1]/time0,
+ maxwellTimeA[2]/time0],
[densityB/density0, muB/mu0, lambdaB/mu0,
shearRatioB[0], shearRatioB[1], shearRatioB[2],
- maxwellTimeB[0]/time0, maxwellTimeB[1]/time0, maxwellTimeB[2]/time0] ],
+ maxwellTimeB[0]/time0, maxwellTimeB[1]/time0,
+ maxwellTimeB[2]/time0] ],
dtype=numpy.float64)
self.stateVarsNondim = self.stateVars # no scaling
@@ -165,18 +173,25 @@
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.stateVarsUpdated = numpy.zeros(
+ (numLocs, tensorSize + numMaxwellModels * tensorSize),
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA, \
+ initialStressA, initialStrainA,
+ self.stateVars[0,:])
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB, \
+ initialStressB, initialStrainB,
+ self.stateVars[1,:])
self.dtStableImplicit = 0.2*min(min(maxwellTimeA), min(maxwellTimeB))
return
- def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV):
+ def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV,
+ stateVars):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -224,9 +239,9 @@
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))
+ strain = numpy.reshape(strainV, (tensorSize,1)) - initialStrain
elastic = numpy.array([ [C1111, C1122, C1133, C1112, C1123, C1113],
[C2211, C2222, C2233, C2212, C2223, C2213],
[C3311, C3322, C3333, C3312, C3323, C3313],
@@ -234,8 +249,14 @@
[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))
+ stress = numpy.dot(elastic, strain) + initialStress
+ meanStrain = (strain[0,0] + strain[1,0] + strain[2,0])/3.0
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0], dtype=numpy.float64)
+ strainVec = numpy.array(strainV, dtype=numpy.float64)
+ viscousStrain = numpy.ravel(strain) - diag * meanStrain
+ stateVarsUpdated = numpy.concatenate((strainVec, viscousStrain,
+ viscousStrain, viscousStrain))
+ return (elasticConsts, numpy.ravel(stress), numpy.ravel(stateVarsUpdated))
# MAIN /////////////////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -342,18 +342,18 @@
};
const double pylith::materials::GenMaxwellIsotropic3DElasticData::_stress[] = {
- -1.57290000e+07,
- -1.12280000e+07,
- -6.72700000e+06,
- 4.52400000e+06,
- 9.02500000e+06,
- 1.35260000e+07,
- -6.14100000e+06,
- -5.56400000e+06,
- -4.98700000e+06,
- -1.04040000e+06,
- -4.63400000e+05,
- 1.13600000e+05,
+ 1.62660000e+07,
+ 2.11720000e+07,
+ 2.60780000e+07,
+ 1.82940000e+07,
+ 2.32000000e+07,
+ 2.81060000e+07,
+ 1.84236000e+06,
+ 2.47120000e+06,
+ 3.10004000e+06,
+ 2.27736000e+06,
+ 2.90620000e+06,
+ 3.53504000e+06,
};
const double pylith::materials::GenMaxwellIsotropic3DElasticData::_elasticConsts[] = {
@@ -447,22 +447,71 @@
};
const double pylith::materials::GenMaxwellIsotropic3DElasticData::_initialStrain[] = {
- 3.10000000e-04,
- 3.20000000e-04,
+ 3.10000000e-05,
+ 3.20000000e-05,
+ 3.30000000e-05,
+ 3.40000000e-05,
+ 3.50000000e-05,
+ 3.60000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.30000000e-05,
+ 6.40000000e-05,
+ 6.50000000e-05,
+ 6.60000000e-05,
+};
+
+const double pylith::materials::GenMaxwellIsotropic3DElasticData::_stateVarsUpdated[] = {
+ 1.10000000e-04,
+ 2.20000000e-04,
3.30000000e-04,
+ 4.40000000e-04,
+ 5.50000000e-04,
+ 6.60000000e-04,
+ -1.09000000e-04,
+ 0.00000000e+00,
+ 1.09000000e-04,
+ 4.06000000e-04,
+ 5.15000000e-04,
+ 6.24000000e-04,
+ -1.09000000e-04,
+ 0.00000000e+00,
+ 1.09000000e-04,
+ 4.06000000e-04,
+ 5.15000000e-04,
+ 6.24000000e-04,
+ -1.09000000e-04,
+ 0.00000000e+00,
+ 1.09000000e-04,
+ 4.06000000e-04,
+ 5.15000000e-04,
+ 6.24000000e-04,
+ 1.20000000e-04,
+ 2.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,
+ 4.50000000e-04,
+ 5.60000000e-04,
+ 6.70000000e-04,
+ -1.09000000e-04,
+ 2.71050543e-20,
+ 1.09000000e-04,
+ 3.86000000e-04,
+ 4.95000000e-04,
+ 6.04000000e-04,
+ -1.09000000e-04,
+ 2.71050543e-20,
+ 1.09000000e-04,
+ 3.86000000e-04,
+ 4.95000000e-04,
+ 6.04000000e-04,
+ -1.09000000e-04,
+ 2.71050543e-20,
+ 1.09000000e-04,
+ 3.86000000e-04,
+ 4.95000000e-04,
+ 6.04000000e-04,
};
-const double* pylith::materials::GenMaxwellIsotropic3DElasticData::_stateVarsUpdated = 0;
-
pylith::materials::GenMaxwellIsotropic3DElasticData::GenMaxwellIsotropic3DElasticData(void)
{ // constructor
dimension = _dimension;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh 2011-05-31 05:27:02 UTC (rev 18497)
@@ -101,7 +101,7 @@
static const double _initialStrain[];
- static const double* _stateVarsUpdated;
+ static const double _stateVarsUpdated[];
};
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py 2011-05-31 05:27:02 UTC (rev 18497)
@@ -24,6 +24,7 @@
from ElasticMaterialApp import ElasticMaterialApp
import numpy
+import math
# ----------------------------------------------------------------------
dimension = 3
@@ -47,6 +48,9 @@
"""
ElasticMaterialApp.__init__(self, name)
+ # import pdb
+ # pdb.set_trace()
+
numLocs = 2
self.dt = 2.0e5
self.dimension = dimension
@@ -55,7 +59,8 @@
self.dbPropertyValues = ["density", "vs", "vp",
"shear-ratio-1", "shear-ratio-2", "shear-ratio-3",
"viscosity-1", "viscosity-2", "viscosity-3"]
- self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1, 1, 1, 1],
+ dtype=numpy.int32)
self.dbStateVarValues = ["total-strain-xx",
"total-strain-yy",
@@ -92,9 +97,10 @@
viscosityA = [1.0e18, 1.0e17, 1.0e19]
strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.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]
+ initialStrainA = [3.1e-5, 3.2e-5, 3.3e-5, 3.4e-5, 3.5e-5, 3.6e-5]
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
+ meanStrainA = (strainA[0] + strainA[1] + strainA[2])/3.0
densityB = 2000.0
vsB = 1200.0
@@ -103,72 +109,61 @@
viscosityB = [1.0e18, 1.0e19, 1.0e20]
strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-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]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.3e-5, 6.4e-5, 6.5e-5, 6.6e-5]
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
+ meanStrainB = (strainB[0] + strainB[1] + strainB[2])/3.0
- # TEMPORARY, need to determine how to use initial strain
- initialStrainA = [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]
-
- # Simplest approach for now is to assume this is the first step after the
- # elastic solution. In that case, both the total strain from the last step
- # (strainT) and the total viscous strain (visStrain) are defined by the
- # assigned elastic strain.
+ # Compute Maxwell time and viscous strain.
diag = numpy.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0], dtype=numpy.float64)
- strainTA = numpy.array(strainA)
- strainTB = numpy.array(strainB)
- meanStrainA = (strainA[0] + strainA[1] + strainA[2])/3.0
- meanStrainB = (strainB[0] + strainB[1] + strainB[2])/3.0
+ strainTA = numpy.array(strainA, dtype=numpy.float64)
+ strainTB = numpy.array(strainB, dtype=numpy.float64)
- maxwellTimeA = [0.0, 0.0, 0.0]
- maxwellTimeB = [0.0, 0.0, 0.0]
- visStrainA = numpy.zeros( (numMaxwellModels, tensorSize), dtype=numpy.float64)
- visStrainB = numpy.zeros( (numMaxwellModels, tensorSize), dtype=numpy.float64)
+ maxwellTimeA = [1.0e30, 1.0e30, 1.0e30]
+ maxwellTimeB = [1.0e30, 1.0e30, 1.0e30]
+ visStrainA = numpy.zeros( (numMaxwellModels, tensorSize),
+ dtype=numpy.float64)
+ visStrainB = numpy.zeros( (numMaxwellModels, tensorSize),
+ dtype=numpy.float64)
for imodel in xrange(numMaxwellModels):
if shearRatioA[imodel] != 0.0:
maxwellTimeA[imodel] = viscosityA[imodel]/(muA*shearRatioA[imodel])
- visStrainA[imodel,:] = strainA[:] - diag[:] * meanStrainA
+ visStrainA[imodel,:] = strainTA[:] - diag[:] * meanStrainA
if shearRatioB[imodel] != 0.0:
maxwellTimeB[imodel] = viscosityB[imodel]/(muB*shearRatioB[imodel])
- visStrainB[imodel,:] = strainB[:] - diag[:] * meanStrainB
+ visStrainB[imodel,:] = strainTB[:] - diag[:] * meanStrainB
- self.lengthScale = 1.0e+3
- self.pressureScale = muA
- self.timeScale = 1.0
- self.densityScale = 1.0e+3
-
- propA = [densityA, vsA, vpA] + shearRatioA + viscosityA
- propB = [densityB, vsB, vpB] + shearRatioB + viscosityB
- self.dbProperties = numpy.array([propA, propB], dtype=numpy.float64)
+ dbPropA = [densityA, vsA, vpA] + shearRatioA + viscosityA
+ dbPropB = [densityB, vsB, vpB] + shearRatioB + viscosityB
+ self.dbProperties = numpy.array([dbPropA, dbPropB], dtype=numpy.float64)
propA = [densityA, muA, lambdaA] + shearRatioA + maxwellTimeA
propB = [densityB, muB, lambdaB] + shearRatioB + maxwellTimeB
self.properties = numpy.array([propA, propB], dtype=numpy.float64)
# TEMPORARY, need to determine how to use initial state variables
- self.dbStateVars = numpy.zeros( (numLocs, tensorSize+numMaxwellModels*tensorSize),
- dtype=numpy.float64)
- self.stateVars = numpy.zeros( (numLocs, tensorSize+numMaxwellModels*tensorSize),
- dtype=numpy.float64)
- self.stateVars[0,0:tensorSize] = strainTA
- self.stateVars[0,tensorSize:(1+numMaxwellModels)*tensorSize] = visStrainA.ravel()
- self.stateVars[1,0:tensorSize] = strainTB
- self.stateVars[1,tensorSize:(1+numMaxwellModels)*tensorSize] = visStrainB.ravel()
+ self.dbStateVars = numpy.zeros(
+ (numLocs, tensorSize + numMaxwellModels * tensorSize),
+ dtype=numpy.float64)
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.timeScale = 1.0
+ self.densityScale = 1.0e+3
+
mu0 = self.pressureScale
density0 = self.densityScale
time0 = self.timeScale
self.propertiesNondim = \
numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0,
shearRatioA[0], shearRatioA[1], shearRatioA[2],
- maxwellTimeA[0]/time0, maxwellTimeA[1]/time0, maxwellTimeA[2]/time0],
+ maxwellTimeA[0]/time0, maxwellTimeA[1]/time0,
+ maxwellTimeA[2]/time0],
[densityB/density0, muB/mu0, lambdaB/mu0,
shearRatioB[0], shearRatioB[1], shearRatioB[2],
- maxwellTimeB[0]/time0, maxwellTimeB[1]/time0, maxwellTimeB[2]/time0] ],
+ maxwellTimeB[0]/time0, maxwellTimeB[1]/time0,
+ maxwellTimeB[2]/time0] ],
dtype=numpy.float64)
- self.stateVarsNondim = self.stateVars # no scaling
-
self.initialStress = numpy.array([initialStressA,
initialStressB],
dtype=numpy.float64)
@@ -180,143 +175,229 @@
densityB],
dtype=numpy.float64)
- self.strain = numpy.array([strainA,
- strainB],
- dtype=numpy.float64)
-
+ # Simplest approach for now is to assume this is the first step
+ # after the elastic solution. In that case, both the total strain
+ # from the last step (total_strain) and the total viscous strain
+ # (viscous_strain) are defined by the assigned elastic strain.
+ # Revised approach. For a better test, I am setting the total strain
+ # for the current time step to be equal to the strain from the previous
+ # time step plus a constant amount.
+ totalStrainA = strainTA + 1.0e-5
+ totalStrainB = strainTB + 1.0e-5
+ visStrainVecA = numpy.ravel(visStrainA)
+ visStrainVecB = numpy.ravel(visStrainB)
+ stateVarsA = numpy.concatenate((strainTA, visStrainVecA))
+ stateVarsB = numpy.concatenate((strainTB, visStrainVecB))
+ self.stateVars = numpy.array([stateVarsA, stateVarsB], dtype=numpy.float64)
+ self.stateVarsNondim = self.stateVars # no scaling
+
+ self.strain = numpy.array([totalStrainA, totalStrainB], dtype=numpy.float64)
self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.stateVarsUpdated = numpy.zeros(
+ (numLocs, tensorSize + numMaxwellModels * tensorSize),
+ dtype=numpy.float64)
+
self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts), \
dtype=numpy.float64)
- (self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, muA, lambdaA, shearRatioA, maxwellTimeA,
- strainTA, visStrainA,
- initialStressA, initialStrainA)
- (self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, muB, lambdaB, shearRatioB, maxwellTimeB,
- strainTB, visStrainB,
- initialStressB, initialStrainB)
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA,
+ shearRatioA, maxwellTimeA,
+ totalStrainA, visStrainA,
+ initialStressA, initialStrainA,
+ stateVarsA)
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB,
+ shearRatioB, maxwellTimeB,
+ totalStrainB, visStrainB,
+ initialStressB, initialStrainB,
+ stateVarsB)
self.dtStableImplicit = 0.2*min(min(maxwellTimeA), min(maxwellTimeB))
return
+ def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
+ muV, lambdaV, shearRatioV, maxwellTimeV, dqV,
+ totalStrainT, viscousStrainT,
+ initialStress, initialStrain,
+ stateVars):
+ """
+ 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, viscousStrainTpdt = self._computeStress(strainTest,
+ muV,
+ lambdaV,
+ shearRatioV,
+ maxwellTimeV,
+ dqV,
+ totalStrainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain,
+ stateVars)
+ return stressTpdt[stressComp]
+
+
+ def _computeStress(self, strainTpdt, muV, lambdaV, shearRatioV, maxwellTimeV,
+ dqV, strainT, viscousStrainT,
+ initialStress, initialStrain, stateVars):
+ """
+ Function to compute stresses and viscous strains for a given strain.
+ """
+
+ bulkModulus = lambdaV + 2.0 * muV / 3.0
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0], dtype=numpy.float64)
+
+ # Initial stresses and strains
+ meanStrainInitial = numpy.dot(diag, initialStrain) / 3.0
+ meanStressInitial = numpy.dot(diag, initialStress) / 3.0
+
+ devStrainInitial = initialStrain - diag * meanStrainInitial
+ devStressInitial = initialStress - diag * meanStressInitial
+
+ # Strains from previous time step
+ meanStrainT = numpy.dot(diag, strainT) / 3.0
+ devStrainT = strainT - diag * meanStrainT - devStrainInitial
+
+ # Strains and mean stress for current time step
+ meanStrainTpdt = numpy.dot(diag, strainTpdt) / 3.0
+ devStrainTpdt = strainTpdt - diag * meanStrainTpdt - devStrainInitial
+ meanStressTpdt = 3.0 * bulkModulus * \
+ (meanStrainTpdt - meanStrainInitial) + meanStressInitial
+
+ # Compute viscous factors
+ visFrac = 0.0
+ visFac = 0.0
+ expFac = numpy.zeros((numMaxwellModels), dtype=numpy.float64)
+ for model in range(numMaxwellModels):
+ visFrac += shearRatioV[model]
+ visFac += shearRatioV[model] * dqV[model]
+ expFac[model] = math.exp(-self.dt/maxwellTimeV[model])
+
+ elasFrac = 1.0 - visFrac
+ stressTpdt = numpy.zeros( (tensorSize), dtype=numpy.float64)
+ viscousStrainTpdt = numpy.zeros( (numMaxwellModels, tensorSize),
+ dtype=numpy.float64)
+ elasFac = 2.0 * muV
+
+ # Compute viscous strain and stress
+ for iComp in range(tensorSize):
+ deltaStrain = devStrainTpdt[iComp] - devStrainT[iComp]
+ devStressTpdt = elasFrac * devStrainTpdt[iComp]
+ for model in range(numMaxwellModels):
+ if shearRatioV[model] != 0.0:
+ viscousStrainTpdt[model, iComp] = expFac[model] * \
+ viscousStrainT[model, iComp] + \
+ dqV[model] * deltaStrain
+ devStressTpdt += shearRatioV[model] * viscousStrainTpdt[model, iComp]
+ devStressTpdt = elasFac * devStressTpdt + devStressInitial[iComp]
+ stressTpdt[iComp] = diag[iComp] * (meanStressTpdt + meanStressInitial) + \
+ devStressTpdt
+
+ return stressTpdt, viscousStrainTpdt
+
+
+ def _computeViscousFactor(self, maxwellTime):
+ """
+ Compute viscous strain factor for a given Maxwell time.
+ """
+
+ timeFrac = 1.0e-5
+ numTerms = 5
+ dq = 0.0
+ if maxwellTime < timeFrac*self.dt:
+ fSign = 1.0
+ factorial = 1.0
+ fraction = 1.0
+ dq = 1.0
+ for iTerm in range(2, numTerms + 1):
+ factorial *= iTerm
+ fSign *= -1.0
+ fraction *= self.dt/maxwellTime
+ dq += fSign*fraction/factorial
+ else:
+ dq = maxwellTime*(1.0-math.exp(-self.dt/maxwellTime))/self.dt
+
+ return dq
+
+
def _calcStress(self, strainV, muV, lambdaV, shearRatioV, maxwellTimeV,
- strainTV, visStrainV, initialStressV, initialStrainV):
+ totalStrainV, viscousStrainV, initialStressV, initialStrainV,
+ stateVars):
"""
Compute stress and derivative of elasticity matrix.
This assumes behavior is always viscoelastic.
"""
- import math
- # import pdb
+ import scipy.misc
- bulkModulus = lambdaV + 2.0 * muV/3.0
+ # Define some numpy arrays
+ strainTpdt = numpy.array(totalStrainV, dtype=numpy.float64)
+ strainT = numpy.array(strainV, dtype=numpy.float64)
+ viscousStrainT = numpy.array(viscousStrainV, dtype=numpy.float64)
+ initialStress = numpy.array(initialStressV, dtype=numpy.float64)
+ initialStrain = numpy.array(initialStrainV, dtype=numpy.float64)
- diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
- traceStrainT = strainTV[0] + strainTV[1] + strainTV[2]
- traceStrainTpdt = strainV[0] + strainV[1] + strainV[2]
- meanStrainT = traceStrainT/3.0
- meanStrainTpdt = traceStrainTpdt/3.0
- meanStressTpdt = bulkModulus * traceStrainTpdt
+ # Compute viscous factors for each Maxwell model
+ dqV = numpy.zeros(numMaxwellModels, dtype=numpy.float64)
+ for model in range(numMaxwellModels):
+ dqV[model] = self._computeViscousFactor(maxwellTimeV[model])
- timeFrac = 1.0e-10
- numTerms = 5
- visFrac = 0.0
- visFac = 0.0
- dq = [0.0, 0.0, 0.0]
- for imodel in range(numMaxwellModels):
- visFrac += shearRatioV[imodel]
- if shearRatioV[imodel] != 0.0:
- if self.dt < timeFrac*maxwellTimeV[imodel]:
- fSign = 1.0
- factorial = 1.0
- fraction = 1.0
- dq[imodel] = 1.0
- for iTerm in range(2, numTerms + 1):
- factorial *= iTerm
- fSign *= -1.0
- fraction *= self.dt/maxwellTimeV[imodel]
- dq[imodel] += fSign*fraction/factorial
- elif maxwellTimeV[imodel] < timeFrac*self.dt:
- dq[imodel] = maxwellTimeV[imodel]/dt
- else:
- dq[imodel] = maxwellTimeV[imodel] * \
- (1.0-math.exp(-self.dt/maxwellTimeV[imodel]))/self.dt
- visFac += shearRatioV[imodel] * dq[imodel]
+ # Compute current stress and viscous strain
+ stressTpdt, viscousStrainTpdt = self._computeStress(strainTpdt,
+ muV,
+ lambdaV,
+ shearRatioV,
+ maxwellTimeV,
+ dqV,
+ strainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain,
+ stateVars)
- elasFrac = 1.0 - visFrac
- shearFac = muV*(elasFrac + visFac)/3.0
+ # Form updated state variables
+ strainTpdtVec = numpy.ravel(strainTpdt)
+ viscousStrainTpdtVec = numpy.ravel(viscousStrainTpdt)
+ stateVarsUpdated = numpy.concatenate((strainTpdtVec, viscousStrainTpdtVec))
- C1111 = bulkModulus + 4.0*shearFac
- C1122 = bulkModulus - 2.0*shearFac
- C1133 = C1122
- C1112 = 0.0
- C1123 = 0.0
- C1113 = 0.0
- C2211 = C1122
- C2222 = C1111
- C2233 = C1122
- C2212 = 0.0
- C2223 = 0.0
- C2213 = 0.0
- C3311 = C1133
- C3322 = C2233
- C3333 = C1111
- C3312 = 0.0
- C3323 = 0.0
- C3313 = 0.0
- C1211 = 0.0
- C1222 = 0.0
- C1233 = 0.0
- C1212 = 6.0*shearFac
- C1223 = 0.0
- C1213 = 0.0
- C2311 = 0.0
- C2322 = 0.0
- C2333 = 0.0
- C2312 = 0.0
- C2323 = C1212
- C2313 = 0.0
- C1311 = 0.0
- C1322 = 0.0
- C1333 = 0.0
- C1312 = 0.0
- C1323 = 0.0
- C1313 = C1212
- 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)
+ # Compute components of tangent constitutive matrix using numerical
+ # derivatives.
+ derivDx = 1.0e-12
+ derivOrder = 3
+ elasticConstsList = []
- stressV = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ for stressComp in range(tensorSize):
+ for strainComp in range(tensorSize):
+ dStressDStrain = stressComp + strainComp
+ dStressDStrain = scipy.misc.derivative(self._calcStressComponent,
+ strainTpdt[strainComp],
+ dx=derivDx,
+ args=(strainComp,
+ stressComp,
+ strainTpdt,
+ muV,
+ lambdaV,
+ shearRatioV,
+ maxwellTimeV,
+ dqV,
+ strainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain,
+ stateVars),
+ order=derivOrder)
+ elasticConstsList.append(dStressDStrain)
- elasFac = 2.0*muV
- devStrainTpdt = 0.0
- devStrainT = 0.0
- deltaStrain = 0.0
- devStressTpdt = 0.0
- visStrain = 0.0
- for iComp in xrange(tensorSize):
- devStrainTpdt = strainV[iComp] - diag[iComp]*meanStrainTpdt
- devStrainT = strainTV[iComp] - diag[iComp]*meanStrainT
- deltaStrain = devStrainTpdt - devStrainT
- devStressTpdt = elasFrac*devStrainTpdt
- for imodel in range(numMaxwellModels):
- if shearRatioV[imodel] != 0.0:
- visStrain = math.exp(-self.dt/maxwellTimeV[imodel]) * \
- visStrainV[imodel,iComp] + \
- dq[imodel] * deltaStrain
- devStressTpdt += shearRatioV[imodel] * visStrain
- devStressTpdt = elasFac * devStressTpdt
- stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt + \
- initialStressV[iComp]
-
- stress = numpy.reshape(stressV, (6,1))
- return (elasticConsts, numpy.ravel(stress))
-
+ elasticConsts = numpy.array(elasticConstsList, dtype=numpy.float64)
+ return (elasticConsts, numpy.ravel(stressTpdt),
+ numpy.ravel(stateVarsUpdated))
+
+
# MAIN /////////////////////////////////////////////////////////////////
if __name__ == "__main__":
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc 2011-05-31 05:27:02 UTC (rev 18497)
@@ -327,108 +327,108 @@
};
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_strain[] = {
- 1.10000000e-04,
- 2.20000000e-04,
- 3.30000000e-04,
- 4.40000000e-04,
- 5.50000000e-04,
- 6.60000000e-04,
1.20000000e-04,
2.30000000e-04,
3.40000000e-04,
4.50000000e-04,
5.60000000e-04,
6.70000000e-04,
+ 1.30000000e-04,
+ 2.40000000e-04,
+ 3.50000000e-04,
+ 4.60000000e-04,
+ 5.70000000e-04,
+ 6.80000000e-04,
};
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_stress[] = {
- 1.98288741e+07,
- 2.47720000e+07,
- 2.97151259e+07,
- 1.97925037e+07,
- 2.47356296e+07,
- 2.96787555e+07,
- 2.72941620e+06,
- 3.36400000e+06,
- 3.99858380e+06,
- 2.64593371e+06,
- 3.28051751e+06,
- 3.91510131e+06,
+ 1.73848741e+07,
+ 2.23190000e+07,
+ 2.72531259e+07,
+ 1.99361456e+07,
+ 2.48702715e+07,
+ 2.98043974e+07,
+ 2.03492020e+06,
+ 2.66720000e+06,
+ 3.29947980e+06,
+ 2.55607698e+06,
+ 3.18835678e+06,
+ 3.82063657e+06,
};
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_elasticConsts[] = {
- 6.74761278e+10,
- 2.25119361e+10,
- 2.25119361e+10,
+ 6.74761292e+10,
+ 2.25119367e+10,
+ 2.25119367e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 2.25119361e+10,
- 6.74761278e+10,
- 2.25119361e+10,
+ 2.25119367e+10,
+ 6.74761292e+10,
+ 2.25119367e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 2.25119361e+10,
- 2.25119361e+10,
- 6.74761278e+10,
+ 2.25119386e+10,
+ 2.25119386e+10,
+ 6.74761292e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 4.49641917e+10,
+ 4.49641924e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 4.49641917e+10,
+ 4.49641924e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 4.49641917e+10,
- 8.63995090e+09,
- 2.88002455e+09,
- 2.88002455e+09,
+ 4.49641906e+10,
+ 8.63995089e+09,
+ 2.88002449e+09,
+ 2.88002449e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 2.88002455e+09,
- 8.63995090e+09,
- 2.88002455e+09,
+ 2.88002449e+09,
+ 8.63995100e+09,
+ 2.88002449e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 2.88002455e+09,
- 2.88002455e+09,
- 8.63995090e+09,
+ 2.88002472e+09,
+ 2.88002472e+09,
+ 8.63995077e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 5.75992635e+09,
+ 5.75992605e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 5.75992635e+09,
+ 5.75992628e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 5.75992635e+09,
+ 5.75992605e+09,
};
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStress[] = {
@@ -447,21 +447,70 @@
};
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStrain[] = {
- 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,
+ 3.10000000e-05,
+ 3.20000000e-05,
+ 3.30000000e-05,
+ 3.40000000e-05,
+ 3.50000000e-05,
+ 3.60000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.30000000e-05,
+ 6.40000000e-05,
+ 6.50000000e-05,
+ 6.60000000e-05,
};
-const double* pylith::materials::GenMaxwellIsotropic3DTimeDepData::_stateVarsUpdated = 0;
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_stateVarsUpdated[] = {
+ 1.20000000e-04,
+ 2.30000000e-04,
+ 3.40000000e-04,
+ 4.50000000e-04,
+ 5.60000000e-04,
+ 6.70000000e-04,
+ -1.09752778e-04,
+ -2.70745840e-20,
+ 1.09752778e-04,
+ 4.48999871e-04,
+ 5.58752650e-04,
+ 6.68505428e-04,
+ -1.09506112e-04,
+ -2.70441593e-20,
+ 1.09506112e-04,
+ 4.48001982e-04,
+ 5.57508094e-04,
+ 6.67014206e-04,
+ -1.09990100e-04,
+ -2.71038346e-20,
+ 1.09990100e-04,
+ 4.49959952e-04,
+ 5.59950052e-04,
+ 6.69940153e-04,
+ 1.30000000e-04,
+ 2.40000000e-04,
+ 3.50000000e-04,
+ 4.60000000e-04,
+ 5.70000000e-04,
+ 6.80000000e-04,
+ -1.09987329e-04,
+ 1.56113122e-24,
+ 1.09987329e-04,
+ 4.59947587e-04,
+ 5.69934916e-04,
+ 6.79922244e-04,
+ -1.09998733e-04,
+ 1.56123808e-25,
+ 1.09998733e-04,
+ 4.59994758e-04,
+ 5.69993491e-04,
+ 6.79992224e-04,
+ -1.09999873e-04,
+ 1.56132337e-26,
+ 1.09999873e-04,
+ 4.59999476e-04,
+ 5.69999349e-04,
+ 6.79999222e-04,
+};
pylith::materials::GenMaxwellIsotropic3DTimeDepData::GenMaxwellIsotropic3DTimeDepData(void)
{ // constructor
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh 2011-05-31 05:27:02 UTC (rev 18497)
@@ -101,7 +101,7 @@
static const double _initialStrain[];
- static const double* _stateVarsUpdated;
+ static const double _stateVarsUpdated[];
};
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py 2011-05-31 05:27:02 UTC (rev 18497)
@@ -46,8 +46,8 @@
"""
ElasticMaterialApp.__init__(self, name)
- import pdb
- pdb.set_trace()
+ # import pdb
+ # pdb.set_trace()
numLocs = 2
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py 2011-05-31 05:27:02 UTC (rev 18497)
@@ -48,8 +48,8 @@
"""
ElasticMaterialApp.__init__(self, name)
- import pdb
- pdb.set_trace()
+ # import pdb
+ # pdb.set_trace()
numLocs = 2
self.dt = 2.0e5
@@ -113,8 +113,8 @@
strainTA = [strainA[0], strainA[1], 0.0, strainA[2]]
strainTB = [strainB[0], strainB[1], 0.0, strainB[2]]
- maxwellTimeA = [0.0, 0.0, 0.0]
- maxwellTimeB = [0.0, 0.0, 0.0]
+ maxwellTimeA = [1.0e30, 1.0e30, 1.0e30]
+ maxwellTimeB = [1.0e30, 1.0e30, 1.0e30]
visStrainA = numpy.zeros( (numMaxwellModels, 4), dtype=numpy.float64)
visStrainB = numpy.zeros( (numMaxwellModels, 4), dtype=numpy.float64)
for imodel in xrange(numMaxwellModels):
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py 2011-05-31 05:27:02 UTC (rev 18497)
@@ -45,8 +45,8 @@
"""
ElasticMaterialApp.__init__(self, name)
- import pdb
- pdb.set_trace()
+ # import pdb
+ # pdb.set_trace()
numLocs = 2
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py 2011-05-31 01:38:18 UTC (rev 18496)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py 2011-05-31 05:27:02 UTC (rev 18497)
@@ -46,8 +46,8 @@
"""
ElasticMaterialApp.__init__(self, name)
- import pdb
- pdb.set_trace()
+ # import pdb
+ # pdb.set_trace()
numLocs = 2
More information about the CIG-COMMITS
mailing list