[cig-commits] r15429 - in short/3D/PyLith/trunk: libsrc/materials unittests/libtests/materials unittests/libtests/materials/data
willic3 at geodynamics.org
willic3 at geodynamics.org
Sun Jul 5 19:30:19 PDT 2009
Author: willic3
Date: 2009-07-05 19:30:18 -0700 (Sun, 05 Jul 2009)
New Revision: 15429
Modified:
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc
Log:
Modifications related primarily to initial stress and strain.
Unit tests using initial strain now all pass, but full-scale testing is
needed.
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2009-07-06 02:30:18 UTC (rev 15429)
@@ -355,9 +355,9 @@
const double s123 = lambda * (e11 + e22 + e33);
- stress[0] = s123 + mu2*e11 + initialStress[0];
- stress[1] = s123 + mu2*e22 + initialStress[1];
- stress[2] = s123 + mu2*e33 + initialStress[2];
+ stress[0] = s123 + mu2 * e11 + initialStress[0];
+ stress[1] = s123 + mu2 * e22 + initialStress[1];
+ stress[2] = s123 + mu2 * e33 + initialStress[2];
stress[3] = mu2 * e12 + initialStress[3];
stress[4] = mu2 * e23 + initialStress[4];
stress[5] = mu2 * e13 + initialStress[5];
@@ -406,15 +406,34 @@
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] +
+ 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 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];
+ 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 };
@@ -433,13 +452,13 @@
double devStressTpdt = 0.0;
for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStressTpdt = mu2 * _viscousStrain[iComp];
+ devStressTpdt = mu2 * (_viscousStrain[iComp] - devStrainInitial[iComp]) +
+ devStressInitial[iComp];
- stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
- initialStress[iComp];
+ stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt;
} // for
- PetscLogFlops(7 + 4 * tensorSize);
+ PetscLogFlops(22 + 5 * tensorSize);
} // _calcStressViscoelastic
// ----------------------------------------------------------------------
@@ -544,8 +563,8 @@
double dq = ViscoelasticMaxwell::viscousStrainParam(_dt, maxwellTime);
const double visFac = mu * dq / 3.0;
- elasticConsts[ 0] = bulkModulus + 4.0*visFac; // C1111
- elasticConsts[ 1] = bulkModulus - 2.0*visFac; // C1122
+ elasticConsts[ 0] = bulkModulus + 4.0 * visFac; // C1111
+ elasticConsts[ 1] = bulkModulus - 2.0 * visFac; // C1122
elasticConsts[ 2] = elasticConsts[1]; // C1133
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
@@ -570,7 +589,7 @@
} // _calcElasticConstsViscoelastic
// ----------------------------------------------------------------------
-// Update state variables.
+// Update state variables as an elastic material.
void
pylith::materials::MaxwellIsotropic3D::_updateStateVarsElastic(
double* const stateVars,
@@ -620,7 +639,7 @@
} // _updateStateVarsElastic
// ----------------------------------------------------------------------
-// Update state variables.
+// Update state variables as a viscoelastic material.
void
pylith::materials::MaxwellIsotropic3D::_updateStateVarsViscoelastic(
double* const stateVars,
@@ -683,7 +702,6 @@
// ----------------------------------------------------------------------
// Compute viscous strain for current time step.
-// material.
void
pylith::materials::MaxwellIsotropic3D::_computeStateVars(
const double* stateVars,
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2009-07-06 02:30:18 UTC (rev 15429)
@@ -15,6 +15,7 @@
#include "MaxwellPlaneStrain.hh" // implementation of object methods
#include "ViscoelasticMaxwell.hh" // USES computeVisStrain
+#include "Metadata.hh" // USES Metadata
#include "pylith/utils/array.hh" // USES double_array
@@ -32,67 +33,113 @@
namespace materials {
namespace _MaxwellPlaneStrain{
+ // Dimension of material
+ const int dimension = 3;
+
/// Number of entries in stress/strain tensors.
- const int tensorSize = 4;
+ const int tensorSize = 3;
/// Number of entries in derivative of elasticity matrix.
const int numElasticConsts = 6;
/// Number of physical properties.
- const int numProperties = 6;
+ const int numProperties = 4;
/// Physical properties.
const Material::PropMetaData properties[] = {
- { "density", 1, OTHER_FIELD },
- { "mu", 1, OTHER_FIELD },
- { "lambda", 1, OTHER_FIELD },
- { "maxwell_time", 1, OTHER_FIELD },
- { "total_strain", 4, OTHER_FIELD },
- { "viscous_strain", 4, OTHER_FIELD },
+ { "density", 1, pylith::topology::FieldBase::SCALAR },
+ { "mu", 1, pylith::topology::FieldBase::SCALAR },
+ { "lambda", 1, pylith::topology::FieldBase::SCALAR },
+ { "maxwell_time", 1, pylith::topology::FieldBase::SCALAR },
};
- /// Indices (order) of properties.
- const int pidDensity = 0;
- const int pidMu = pidDensity + 1;
- const int pidLambda = pidMu + 1;
- const int pidMaxwellTime = pidLambda + 1;
- const int pidStrainT = pidMaxwellTime + 1;
- const int pidVisStrain = pidStrainT + tensorSize;
- /// Values expected in spatial database
- const int numDBValues = 4;
- const char* namesDBValues[] = {"density", "vs", "vp" , "viscosity"};
+ /// Values expected in properties spatial database
+ const int numDBProperties = 4;
+ const char* dbProperties[] = {"density", "vs", "vp" , "viscosity"};
- /// Indices (order) of database values
- const int didDensity = 0;
- const int didVs = 1;
- const int didVp = 2;
- const int didViscosity = 3;
+ /// Number of state variables.
+ const int numStateVars = 2;
- /// Initial state values expected in spatial database
- const int numInitialStateDBValues = tensorSize;
- const char* namesInitialStateDBValues[] = { "stress_xx", "stress_yy",
- "stress_zz", "stress_xy" };
+ /// Size of viscous_strain state variable.
+ const int viscousStrainSize = 4;
+ /// State variables.
+ const Metadata::ParamDescription stateVars[] = {
+ { "total_strain", tensorSize, pylith::topology::FieldBase::TENSOR },
+ { "viscous_strain", viscousStrainSize,
+ pylith::topology::FieldBase::TENSOR },
+ };
+
+ /// Values expected in state variables spatial database
+ const int numDBStateVars = tensorSize + viscousStrainSize;
+ const char* dbStateVars[] = { "total-strain-xx",
+ "total-strain-yy",
+ "total-strain-xy",
+ "viscous-strain-xx",
+ "viscous-strain-yy",
+ "viscous-strain-zz",
+ "viscous-strain-xy" };
+
} // _MaxwellPlaneStrain
} // materials
} // pylith
+// Indices of physical properties
+const int pylith::materials::MaxwellPlaneStrain::p_density = 0;
+
+const int pylith::materials::MaxwellPlaneStrain::p_mu =
+ pylith::materials::MaxwellPlaneStrain::p_density + 1;
+
+const int pylith::materials::MaxwellPlaneStrain::p_lambda =
+ pylith::materials::MaxwellPlaneStrain::p_mu + 1;
+
+const int pylith::materials::MaxwellPlaneStrain::p_maxwellTime =
+ pylith::materials::MaxwellPlaneStrain::p_lambda + 1;
+
+// Indices of database values (order must match dbProperties)
+const int pylith::materials::MaxwellPlaneStrain::db_density = 0;
+
+const int pylith::materials::MaxwellPlaneStrain::db_vs =
+ pylith::materials::MaxwellPlaneStrain::db_density + 1;
+
+const int pylith::materials::MaxwellPlaneStrain::db_vp =
+ pylith::materials::MaxwellPlaneStrain::db_vs + 1;
+
+const int pylith::materials::MaxwellPlaneStrain::db_viscosity =
+ pylith::materials::MaxwellPlaneStrain::db_vp + 1;
+
+// Indices of state variables
+const int pylith::materials::MaxwellPlaneStrain::s_totalStrain = 0;
+
+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_viscousStrain =
+ pylith::materials::MaxwellPlaneStrain::db_totalStrain + 3;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::MaxwellPlaneStrain::MaxwellPlaneStrain(void) :
- ElasticMaterial(_MaxwellPlaneStrain::tensorSize,
+ ElasticMaterial(_MaxwellPlaneStrain::dimension,
+ _MaxwellPlaneStrain::tensorSize,
_MaxwellPlaneStrain::numElasticConsts,
- _MaxwellPlaneStrain::namesDBValues,
- _MaxwellPlaneStrain::namesInitialStateDBValues,
- _MaxwellPlaneStrain::numDBValues,
- _MaxwellPlaneStrain::properties,
- _MaxwellPlaneStrain::numProperties),
- _calcElasticConstsFn(&pylith::materials::MaxwellPlaneStrain::_calcElasticConstsElastic),
- _calcStressFn(&pylith::materials::MaxwellPlaneStrain::_calcStressElastic),
- _updatePropertiesFn(&pylith::materials::MaxwellPlaneStrain::_updatePropertiesElastic)
+ Metadata(_MaxwellPlaneStrain::properties,
+ _MaxwellPlaneStrain::numProperties,
+ _MaxwellPlaneStrain::dbProperties,
+ _MaxwellPlaneStrain::numDBProperties,
+ _MaxwellPlaneStrain::stateVars,
+ _MaxwellPlaneStrain::numStateVars,
+ _MaxwellPlaneStrain::dbStateVars,
+ _MaxwellPlaneStrain::numDBStateVars)),
+ _calcElasticConstsFn(0),
+ _calcStressFn(0),
+ _updateStateVarsFn(0)
{ // constructor
- _dimension = 2;
- _visStrain.resize(_MaxwellPlaneStrain::tensorSize);
+ useElasticBehavior(true);
+ _viscousStrain.resize(tensorSize);
} // constructor
// ----------------------------------------------------------------------
@@ -102,6 +149,29 @@
} // destructor
// ----------------------------------------------------------------------
+// Set whether elastic or inelastic constitutive relations are used.
+void
+pylith::materials::MaxwellPlaneStrain::useElasticBehavior(const bool flag)
+{ // useElasticBehavior
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcElasticConstsElastic;
+ _updateStateVarsFn =
+ &pylith::materials::MaxwellPlaneStrain::_updateStateVarsElastic;
+
+ } else {
+ _calcStressFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcStressViscoelastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcElasticConstsViscoelastic;
+ _updateStateVarsFn =
+ &pylith::materials::MaxwellPlaneStrain::_updateStateVarsViscoelastic;
+ } // if/else
+} // useElasticBehavior
+
+// ----------------------------------------------------------------------
// Compute properties from values in spatial database.
void
pylith::materials::MaxwellPlaneStrain::_dbToProperties(
@@ -110,12 +180,12 @@
{ // _dbToProperties
assert(0 != propValues);
const int numDBValues = dbValues.size();
- assert(_MaxwellPlaneStrain::numDBValues == numDBValues);
+ assert(_MaxwellPlaneStrain::numDBProperties == numDBValues);
- const double density = dbValues[_MaxwellPlaneStrain::didDensity];
- const double vs = dbValues[_MaxwellPlaneStrain::didVs];
- const double vp = dbValues[_MaxwellPlaneStrain::didVp];
- const double viscosity = dbValues[_MaxwellPlaneStrain::didViscosity];
+ const double density = dbValues[db_density];
+ const double vs = dbValues[db_vs];
+ const double vp = dbValues[db_vp];
+ const double viscosity = dbValues[db_viscosity];
if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || viscosity <= 0.0) {
std::ostringstream msg;
@@ -141,13 +211,13 @@
} // if
assert(mu > 0);
- const double maxwelltime = viscosity / mu;
- assert(maxwelltime > 0.0);
+ const double maxwellTime = viscosity / mu;
+ assert(maxwellTime > 0.0);
- propValues[_MaxwellPlaneStrain::pidDensity] = density;
- propValues[_MaxwellPlaneStrain::pidMu] = mu;
- propValues[_MaxwellPlaneStrain::pidLambda] = lambda;
- propValues[_MaxwellPlaneStrain::pidMaxwellTime] = maxwelltime;
+ propValues[p_density] = density;
+ propValues[p_mu] = mu;
+ propValues[p_lambda] = lambda;
+ propValues[p_maxwellTime] = maxwellTime;
PetscLogFlops(7);
} // _dbToProperties
@@ -160,23 +230,19 @@
{ // _nondimProperties
assert(0 != _normalizer);
assert(0 != values);
- assert(nvalues == _totalPropsQuadPt);
+ assert(nvalues == _numPropsQuadPt);
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
const double timeScale = _normalizer->timeScale();
- values[_MaxwellPlaneStrain::pidDensity] =
- _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidDensity],
- densityScale);
- values[_MaxwellPlaneStrain::pidMu] =
- _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidMu],
- pressureScale);
- values[_MaxwellPlaneStrain::pidLambda] =
- _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidLambda],
- pressureScale);
- values[_MaxwellPlaneStrain::pidMaxwellTime] =
- _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidMaxwellTime],
- timeScale);
+ values[p_density] =
+ _normalizer->nondimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->nondimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->nondimensionalize(values[p_lambda], pressureScale);
+ values[p_maxwellTime] =
+ _normalizer->nondimensionalize(values[p_maxwellTime], timeScale);
PetscLogFlops(4);
} // _nondimProperties
@@ -189,128 +255,59 @@
{ // _dimProperties
assert(0 != _normalizer);
assert(0 != values);
- assert(nvalues == _totalPropsQuadPt);
+ assert(nvalues == _numPropsQuadPt);
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
const double timeScale = _normalizer->timeScale();
- values[_MaxwellPlaneStrain::pidDensity] =
- _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidDensity],
- densityScale);
- values[_MaxwellPlaneStrain::pidMu] =
- _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidMu],
- pressureScale);
- values[_MaxwellPlaneStrain::pidLambda] =
- _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidLambda],
- pressureScale);
- values[_MaxwellPlaneStrain::pidMaxwellTime] =
- _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidMaxwellTime],
- timeScale);
+ values[p_density] =
+ _normalizer->dimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->dimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->dimensionalize(values[p_lambda], pressureScale);
+ values[p_maxwellTime] =
+ _normalizer->dimensionalize(values[p_maxwellTime], timeScale);
PetscLogFlops(4);
} // _dimProperties
// ----------------------------------------------------------------------
-// Nondimensionalize initial state.
+// Compute initial state variables from values in spatial database.
void
-pylith::materials::MaxwellPlaneStrain::_nondimInitState(double* const values,
- const int nvalues) const
-{ // _nondimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _MaxwellPlaneStrain::numInitialStateDBValues);
+pylith::materials::MaxwellPlaneStrain::_dbToStateVars(
+ double* const stateValues,
+ const double_array& dbValues) const
+{ // _dbToStateVars
+ assert(0 != stateValues);
+ const int numDBValues = dbValues.size();
+ assert(_MaxwellPlaneStrain::numDBStateVars == numDBValues);
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->nondimensionalize(values, nvalues, pressureScale);
+ const int totalSize = _tensorSize + 4;
+ assert(totalSize == _numVarsQuadPt);
+ assert(totalSize == numDBValues);
+ memcpy(stateValues, &dbValues[0], totalSize*sizeof(double));
- PetscLogFlops(nvalues);
-} // _nondimInitState
+ PetscLogFlops(0);
+} // _dbToStateVars
// ----------------------------------------------------------------------
-// Dimensionalize initial state.
-void
-pylith::materials::MaxwellPlaneStrain::_dimInitState(double* const values,
- const int nvalues) const
-{ // _dimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _MaxwellPlaneStrain::numInitialStateDBValues);
-
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->dimensionalize(values, nvalues, pressureScale);
-
- PetscLogFlops(nvalues);
-} // _dimInitState
-
-// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::MaxwellPlaneStrain::_calcDensity(double* const density,
const double* properties,
- const int numProperties)
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars)
{ // _calcDensity
assert(0 != density);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
- density[0] = properties[_MaxwellPlaneStrain::pidDensity];
+ density[0] = properties[p_density];
} // _calcDensity
// ----------------------------------------------------------------------
-// Compute viscous strain for current time step.
-// material.
-void
-pylith::materials::MaxwellPlaneStrain::_computeStateVars(
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _computeStateVars
- assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
- assert(0 != totalStrain);
- assert(_MaxwellPlaneStrain::tensorSize == strainSize);
- assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
-
- const int tensorSize = _MaxwellPlaneStrain::tensorSize;
- const double maxwelltime = properties[_MaxwellPlaneStrain::pidMaxwellTime];
-
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
- const double e12 = totalStrain[3];
-
- const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
-
- const double diag[] = { 1.0, 1.0, 1.0, 0.0 };
-
- const double meanStrainT =
- (properties[_MaxwellPlaneStrain::pidStrainT+0] +
- properties[_MaxwellPlaneStrain::pidStrainT+1] +
- properties[_MaxwellPlaneStrain::pidStrainT+2])/3.0;
-
- // Time integration.
- double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
- const double expFac = exp(-_dt/maxwelltime);
-
- double devStrainTpdt = 0.0;
- double devStrainT = 0.0;
-
- for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
- devStrainT = properties[_MaxwellPlaneStrain::pidStrainT+iComp] -
- diag[iComp] * meanStrainT;
- _visStrain[iComp] = expFac *
- properties[_MaxwellPlaneStrain::pidVisStrain + iComp] +
- dq * (devStrainTpdt - devStrainT);
- } // for
-
- PetscLogFlops(8 + 7 * tensorSize);
-} // _computeStateVars
-
-// ----------------------------------------------------------------------
// Compute stress tensor at location from properties as an elastic
// material.
void
@@ -319,38 +316,44 @@
const int stressSize,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars)
{ // _calcStressElastic
assert(0 != stress);
assert(_MaxwellPlaneStrain::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
- assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStrainSize);
- const double mu = properties[_MaxwellPlaneStrain::pidMu];
- const double lambda = properties[_MaxwellPlaneStrain::pidLambda];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0 * mu;
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
- const double e12 = totalStrain[3];
+ const double e11 = totalStrain[0] - initialStrain[0];
+ const double e22 = totalStrain[1] - initialStrain[1];
+ const double e12 = totalStrain[2] - initialStrain[2];
- const double traceStrainTpdt = e11 + e22 + e33;
- const double s123 = lambda * traceStrainTpdt;
+ const double s12 = lambda * (e11 + e22);
- stress[0] = s123 + mu2*e11 + initialState[0];
- stress[1] = s123 + mu2*e22 + initialState[1];
- stress[2] = s123 + initialState[2];
- stress[3] = mu2 * e12 + initialState[3];
+ stress[0] = s12 + mu2 * e11 + initialStress[0];
+ stress[1] = s12 + mu2 * e22 + initialStress[1];
+ stress[2] = mu2 * e12 + initialStress[2];
- PetscLogFlops(13);
+ PetscLogFlops(14);
} // _calcStressElastic
// ----------------------------------------------------------------------
@@ -362,32 +365,42 @@
const int stressSize,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars)
{ // _calcStressViscoelastic
assert(0 != stress);
assert(_MaxwellPlaneStrain::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numPropsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
- assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStrainSize);
const int tensorSize = _MaxwellPlaneStrain::tensorSize;
- const double mu = properties[_MaxwellPlaneStrain::pidMu];
- const double lambda = properties[_MaxwellPlaneStrain::pidLambda];
- const double maxwelltime = properties[_MaxwellPlaneStrain::pidMaxwellTime];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+ const double maxwellTime = properties[p_maxwellTime];
const double mu2 = 2.0 * mu;
- const double bulkModulus = lambda + mu2/3.0;
+ const double bulkModulus = lambda + mu2 / 3.0;
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
+ const double e11 = totalStrain[0] - initialStrain[0];
+ const double e22 = totalStrain[1] - initialStrain[1];
+ const double e12 = totalStrain[2] - initialStrain[2];
+ const double e33 = 0.0;
const double traceStrainTpdt = e11 + e22 + e33;
const double meanStrainTpdt = traceStrainTpdt/3.0;
@@ -397,28 +410,27 @@
// Get viscous strains
if (computeStateVars) {
- pylith::materials::MaxwellPlaneStrain::_computeStateVars(properties,
- numProperties,
- totalStrain,
- strainSize,
- initialState,
- initialStateSize);
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
} else {
- memcpy(&_visStrain[0], &properties[_MaxwellPlaneStrain::pidVisStrain],
- tensorSize * sizeof(double));
+ memcpy(&_viscousStrain[0], &stateVars[s_viscousStrain],
+ 4 * sizeof(double));
} // else
// Compute new stresses
double devStressTpdt = 0.0;
for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStressTpdt = mu2 * _visStrain[iComp];
+ devStressTpdt = mu2 * _viscousStrain[iComp];
stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
initialState[iComp];
} // for
- PetscLogFlops(7 + 4 * tensorSize);
+ PetscLogFlops(10 + 4 * tensorSize);
} // _calcStressViscoelastic
// ----------------------------------------------------------------------
@@ -429,174 +441,270 @@
const int numElasticConsts,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConstsElastic
assert(0 != elasticConsts);
assert(_MaxwellPlaneStrain::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numPropsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
- assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStrainSize);
- const double mu = properties[_MaxwellPlaneStrain::pidMu];
- const double lambda = properties[_MaxwellPlaneStrain::pidLambda];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
elasticConsts[ 0] = lambda2mu; // C1111
elasticConsts[ 1] = lambda; // C1122
- elasticConsts[ 2] = lambda; // C1133
- elasticConsts[ 3] = 0; // C1112
- elasticConsts[ 4] = lambda2mu; // C2222
- elasticConsts[ 5] = lambda; // C2233
- elasticConsts[ 6] = 0; // C2212
- elasticConsts[ 7] = lambda2mu; // C3333
- elasticConsts[ 8] = 0; // C3312
- elasticConsts[ 9] = mu2; // C1212
+ elasticConsts[ 2] = 0; // C1112
+ elasticConsts[ 3] = lambda2mu; // C2222
+ elasticConsts[ 4] = 0; // C2212
+ elasticConsts[ 5] = mu2; // C1212
PetscLogFlops(2);
} // _calcElasticConstsElastic
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from properties
-// as an elastic material.
+// as a viscoelastic material.
void
pylith::materials::MaxwellPlaneStrain::_calcElasticConstsViscoelastic(
double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConstsViscoelastic
assert(0 != elasticConsts);
assert(_MaxwellPlaneStrain::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numPropsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
- assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStrainSize);
- const double mu = properties[_MaxwellPlaneStrain::pidMu];
- const double lambda = properties[_MaxwellPlaneStrain::pidLambda];
- const double maxwelltime = properties[_MaxwellPlaneStrain::pidMaxwellTime];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+ const double maxwellTime = properties[p_maxwellTime];
const double mu2 = 2.0 * mu;
- const double bulkModulus = lambda + mu2/3.0;
+ const double bulkModulus = lambda + mu2 / 3.0;
- double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
+ double dq = ViscoelasticMaxwell::viscousStrainParam(_dt, maxwellTime);
- const double visFac = mu*dq/3.0;
- elasticConsts[ 0] = bulkModulus + 4.0*visFac; // C1111
- elasticConsts[ 1] = bulkModulus - 2.0*visFac; // C1122
- elasticConsts[ 2] = elasticConsts[1]; // C1133
- elasticConsts[ 3] = 0; // C1112
- elasticConsts[ 4] = elasticConsts[0]; // C2222
- elasticConsts[ 5] = elasticConsts[1]; // C2233
- elasticConsts[ 6] = 0; // C2212
- elasticConsts[ 7] = elasticConsts[0]; // C3333
- elasticConsts[ 8] = 0; // C3312
- elasticConsts[ 9] = 6.0 * visFac; // C1212
+ const double visFac = mu * dq / 3.0;
+ elasticConsts[ 0] = bulkModulus + 4.0 * visFac; // C1111
+ elasticConsts[ 1] = bulkModulus - 2.0 * visFac; // C1122
+ elasticConsts[ 2] = 0; // C1112
+ elasticConsts[ 3] = elasticConsts[0]; // C2222
+ elasticConsts[ 4] = 0; // C2212
+ elasticConsts[ 5] = 6.0 * visFac; // C1212
PetscLogFlops(10);
} // _calcElasticConstsViscoelastic
// ----------------------------------------------------------------------
-// Get stable time step for implicit time integration.
-double
-pylith::materials::MaxwellPlaneStrain::_stableTimeStepImplicit(const double* properties,
- const int numProperties) const
-{ // _stableTimeStepImplicit
- assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
-
- const double maxwellTime =
- properties[_MaxwellPlaneStrain::pidMaxwellTime];
- const double dtStable = 0.1*maxwellTime;
-
- return dtStable;
-} // _stableTimeStepImplicit
-
-// ----------------------------------------------------------------------
-// Update state variables.
+// Update state variables as an elastic material.
void
-pylith::materials::MaxwellPlaneStrain::_updatePropertiesElastic(
- double* const properties,
+pylith::materials::MaxwellPlaneStrain::_updateStateVarsElastic(
+ double* const stateVars,
+ const int numStateVars,
+ const double* properties,
const int numProperties,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _updatePropertiesElastic
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVarsElastic
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStrainSize);
- const double maxwelltime = properties[_MaxwellPlaneStrain::pidMaxwellTime];
+ const int tensorSize = _tensorSize;
+ const double maxwellTime = properties[p_maxwellTime];
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
+ const double strainTpdt[] = {totalStrain[0],
+ totalStrain[1],
+ 0.0,
+ totalStrain[2]};
- const double traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
+ const double traceStrainTpdt = strainTpdt[0] + strainTpdt[1] + strainTpdt[2];
+ const double meanStrainTpdt = traceStrainTpdt / 3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0 };
- for (int iComp=0; iComp < _MaxwellPlaneStrain::tensorSize; ++iComp) {
- properties[_MaxwellPlaneStrain::pidStrainT+iComp] = totalStrain[iComp];
- properties[_MaxwellPlaneStrain::pidVisStrain+iComp] =
- totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ stateVars[s_totalStrain] = totalStrain[0];
+ stateVars[s_totalStrain + 1] = totalStrain[1];
+ stateVars[s_totalStrain + 2] = totalStrain[2];
+
+ for (int iComp=0; iComp < 4; ++iComp) {
+ stateVars[s_viscousStrain + iComp] =
+ strainTpdt[iComp] - diag[iComp] * meanStrainTpdt;
} // for
- PetscLogFlops(3 + 2 * _MaxwellPlaneStrain::tensorSize);
+ PetscLogFlops(11);
_needNewJacobian = true;
-} // _updatePropertiesElastic
+} // _updateStateVarsElastic
// ----------------------------------------------------------------------
-// Update state variables.
+// Update state variables as a viscoelastic material.
void
-pylith::materials::MaxwellPlaneStrain::_updatePropertiesViscoelastic(
- double* const properties,
+pylith::materials::MaxwellPlaneStrain::_updateStateVarsViscoelastic(
+ double* const stateVars,
+ const int numStateVars,
+ const double* properties,
const int numProperties,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _updatePropertiesViscoelastic
+ const double* initialStress,
+ const double* initialStress,
+ const int initialStrainSize,
+ const int initialStrainSize)
+{ // _updateStateVarsViscoelastic
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
- assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStrainSize);
- const int tensorSize = _MaxwellPlaneStrain::tensorSize;
+ const int tensorSize = _tensorSize;
- pylith::materials::MaxwellPlaneStrain::_computeStateVars(properties,
- numProperties,
- totalStrain,
- strainSize,
- initialState,
- initialStateSize);
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
- memcpy(&properties[_MaxwellPlaneStrain::pidVisStrain],
- &_visStrain[0],
- tensorSize * sizeof(double));
- memcpy(&properties[_MaxwellPlaneStrain::pidStrainT],
- &totalStrain[0],
- tensorSize * sizeof(double));
+ memcpy(&stateVars[s_totalStrain], totalStrain, tensorSize * sizeof(double));
+ memcpy(&stateVars[s_viscousStrain], &_viscousStrain[0], 4 * sizeof(double));
+
_needNewJacobian = false;
-} // _updatePropertiesViscoelastic
+} // _updateStateVarsViscoelastic
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+pylith::materials::MaxwellPlaneStrain::_stableTimeStepImplicit(
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const
+{ // _stableTimeStepImplicit
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ const double maxwellTime = properties[p_maxwellTime];
+ const double dtStable = 0.1*maxwellTime;
+
+ return dtStable;
+} // _stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Compute viscous strain for current time step.
+void
+pylith::materials::MaxwellPlaneStrain::_computeStateVars(
+ const double* stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _computeStateVars
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ const int tensorSize = _tensorSize;
+ const double maxwellTime = properties[p_maxwellTime];
+
+ const double e11 = totalStrain[0];
+ const double e22 = totalStrain[1];
+ const double e33 = totalStrain[2];
+ const double e12 = totalStrain[3];
+
+ const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0 };
+
+ const double meanStrainT =
+ (properties[_MaxwellPlaneStrain::pidStrainT+0] +
+ properties[_MaxwellPlaneStrain::pidStrainT+1] +
+ properties[_MaxwellPlaneStrain::pidStrainT+2])/3.0;
+
+ // Time integration.
+ double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
+ const double expFac = exp(-_dt/maxwelltime);
+
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ devStrainT = properties[_MaxwellPlaneStrain::pidStrainT+iComp] -
+ diag[iComp] * meanStrainT;
+ _visStrain[iComp] = expFac *
+ properties[_MaxwellPlaneStrain::pidVisStrain + iComp] +
+ dq * (devStrainTpdt - devStrainT);
+ } // for
+
+ PetscLogFlops(8 + 7 * tensorSize);
+} // _computeStateVars
+
+
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2009-07-06 02:30:18 UTC (rev 15429)
@@ -122,66 +122,10 @@
void
pylith::materials::TestMaxwellIsotropic3D::test_updateStateVarsElastic(void)
{ // test_updateStateVarsElastic
- // :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);
- MaxwellIsotropic3D 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 };
- for (int i=0; i < tensorSize; ++i)
- stateVarsE[s_viscousStrain+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();
} // test_updateStateVarsElastic
// ----------------------------------------------------------------------
@@ -219,90 +163,15 @@
void
pylith::materials::TestMaxwellIsotropic3D::test_updateStateVarsTimeDep(void)
{ // test_updateStateVarsTimeDep
- // :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);
- MaxwellIsotropic3D material;
- material.useElasticBehavior(false);
-
delete _dataElastic; _dataElastic = new MaxwellIsotropic3DTimeDepData();
- 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));
-
- 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 meanStrainTpdt =
- (strain[0] + strain[1] + strain[2]) / 3.0;
- const double meanStrainT =
- (stateVars[s_totalStrain+0] +
- stateVars[s_totalStrain+1] +
- stateVars[s_totalStrain+2]) / 3.0;
-
- const int p_maxwellTime = 3;
- const double maxwellTime = properties[p_maxwellTime];
- const double dq = maxwellTime*(1.0-exp(-dt/maxwellTime))/dt;
- const double expFac = exp(-dt/maxwellTime);
- double devStrainTpdt = 0.0;
- double devStrainT = 0.0;
-
- const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
- for (int i=0; i < tensorSize; ++i) {
- devStrainTpdt = strain[i] - diag[i]*meanStrainTpdt;
- devStrainT = stateVars[s_totalStrain+i] - diag[i]*meanStrainT;
- stateVarsE[s_viscousStrain+i] =
- expFac * stateVars[s_viscousStrain+i] +
- dq * (devStrainTpdt - devStrainT);
- } //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)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(stateVarsE[i], stateVars[i], tolerance);
- } // for
} // test_updateStateVarsTimeDep
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc 2009-07-06 02:30:18 UTC (rev 15429)
@@ -185,5 +185,32 @@
TestElasticMaterial::test_stableTimeStepImplicit();
} // test_stableTimeStepImplicit
+// ----------------------------------------------------------------------
+// Test hasProperty().
+void
+pylith::materials::TestPowerLaw3D::testHasProperty(void)
+{ // testHasProperty
+ PowerLaw3D material;
+ CPPUNIT_ASSERT(material.hasProperty("mu"));
+ CPPUNIT_ASSERT(material.hasProperty("lambda"));
+ CPPUNIT_ASSERT(material.hasProperty("density"));
+ CPPUNIT_ASSERT(material.hasProperty("viscosity_coeff"));
+ CPPUNIT_ASSERT(material.hasProperty("power_law_exponent"));
+ CPPUNIT_ASSERT(!material.hasProperty("aaa"));
+} // testHasProperty
+
+// ----------------------------------------------------------------------
+// Test hasStateVar().
+void
+pylith::materials::TestPowerLaw3D::testHasStateVar(void)
+{ // testHasStateVar
+ PowerLaw3D material;
+
+ CPPUNIT_ASSERT(material.hasStateVar("viscous_strain"));
+ CPPUNIT_ASSERT(material.hasStateVar("stress"));
+ CPPUNIT_ASSERT(!material.hasStateVar("total_strain"));
+} // testHasStateVar
+
+
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.hh 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.hh 2009-07-06 02:30:18 UTC (rev 15429)
@@ -62,6 +62,9 @@
CPPUNIT_TEST( test_updateStateVarsElastic );
CPPUNIT_TEST( test_updateStateVarsTimeDep );
+ CPPUNIT_TEST( testHasProperty );
+ CPPUNIT_TEST( testHasStateVar );
+
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
@@ -100,6 +103,12 @@
/// Test _stableTimeStepImplicit()
void test_stableTimeStepImplicit(void);
+ /// Test hasProperty()
+ void testHasProperty(void);
+
+ /// Test hasStateVar()
+ void testHasStateVar(void);
+
}; // class TestPowerLaw3D
#endif // pylith_materials_testpowerlaw3d_hh
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py 2009-07-06 02:30:18 UTC (rev 15429)
@@ -68,7 +68,7 @@
viscosityA = 1.0e18
strainA = [1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4]
initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
- initialStrainA = [3.1e-4, 3.2e-4, 3.3e-4, 3.4e-4, 3.5e-4, 3.6e-4]
+ 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
maxwellTimeA = viscosityA / muA
@@ -79,7 +79,7 @@
viscosityB = 1.0e18
strainB = [4.1e-4, 4.2e-4, 4.3e-4, 4.4e-4, 4.5e-4, 4.6e-4]
initialStressB = [5.1e4, 5.2e4, 5.3e4, 5.4e4, 5.5e4, 5.6e4]
- initialStrainB = [6.1e-4, 6.2e-4, 6.3e-4, 6.4e-4, 6.5e-4, 6.6e-4]
+ 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
maxwellTimeB = viscosityB / muB
@@ -128,13 +128,15 @@
dtype=numpy.float64)
self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
- self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts), \
- dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts), \
+ dtype=numpy.float64)
+ self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + tensorSize), \
+ dtype=numpy.float64)
- (self.elasticConsts[0,:], self.stress[0,:]) = \
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
self._calcStress(strainA, muA, lambdaA, \
initialStressA, initialStrainA)
- (self.elasticConsts[1,:], self.stress[1,:]) = \
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
self._calcStress(strainB, muB, lambdaB, \
initialStressB, initialStrainB)
self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
@@ -184,7 +186,16 @@
[C1113, C2213, C3313, C1213, C2313, C1313] ],
dtype=numpy.float64)
stress = numpy.dot(elastic, strain-initialStrain) + initialStress
- return (elasticConsts, numpy.ravel(stress))
+ meanStrain = (strain[0] + strain[1] + strain[2])/3.0
+ viscousStrain = [strain[0] - meanStrain,
+ strain[1] - meanStrain,
+ strain[2] - meanStrain,
+ strain[3],
+ strain[4],
+ strain[5]]
+ stateVarsUpdated = numpy.array( [strain, viscousStrain],
+ dtype=numpy.float64)
+ return (elasticConsts, numpy.ravel(stress), numpy.ravel(stateVarsUpdated))
# MAIN /////////////////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2009-07-06 02:30:18 UTC (rev 15429)
@@ -210,18 +210,18 @@
};
const double pylith::materials::MaxwellIsotropic3DElasticData::_stress[] = {
- -2.24790000e+07,
- -2.24780000e+07,
- -2.24770000e+07,
- -8.97600000e+06,
- -8.97500000e+06,
- -8.97400000e+06,
- -2.82900000e+06,
- -2.82800000e+06,
- -2.82700000e+06,
- -1.09800000e+06,
- -1.09700000e+06,
- -1.09600000e+06,
+ 9.51600000e+06,
+ 9.92200000e+06,
+ 1.03280000e+07,
+ 4.79400000e+06,
+ 5.20000000e+06,
+ 5.60600000e+06,
+ 5.15436000e+06,
+ 5.20720000e+06,
+ 5.26004000e+06,
+ 2.21976000e+06,
+ 2.27260000e+06,
+ 2.32544000e+06,
};
const double pylith::materials::MaxwellIsotropic3DElasticData::_elasticConsts[] = {
@@ -285,21 +285,46 @@
};
const double pylith::materials::MaxwellIsotropic3DElasticData::_initialStrain[] = {
- 3.10000000e-04,
- 3.20000000e-04,
- 3.30000000e-04,
- 3.40000000e-04,
- 3.50000000e-04,
- 3.60000000e-04,
- 6.10000000e-04,
- 6.20000000e-04,
- 6.30000000e-04,
- 6.40000000e-04,
- 6.50000000e-04,
- 6.60000000e-04,
+ 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::MaxwellIsotropic3DElasticData::_stateVarsUpdated = 0;
+const double pylith::materials::MaxwellIsotropic3DElasticData::_stateVarsUpdated[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.30000000e-04,
+ 1.40000000e-04,
+ 1.50000000e-04,
+ 1.60000000e-04,
+ -1.00000000e-05,
+ 1.35525272e-20,
+ 1.00000000e-05,
+ 1.40000000e-04,
+ 1.50000000e-04,
+ 1.60000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.30000000e-04,
+ 4.40000000e-04,
+ 4.50000000e-04,
+ 4.60000000e-04,
+ -1.00000000e-05,
+ 0.00000000e+00,
+ 1.00000000e-05,
+ 4.40000000e-04,
+ 4.50000000e-04,
+ 4.60000000e-04,
+};
pylith::materials::MaxwellIsotropic3DElasticData::MaxwellIsotropic3DElasticData(void)
{ // constructor
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2009-07-06 02:30:18 UTC (rev 15429)
@@ -95,7 +95,7 @@
static const double _initialStrain[];
- static const double* _stateVarsUpdated;
+ static const double _stateVarsUpdated[];
};
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py 2009-07-06 02:30:18 UTC (rev 15429)
@@ -39,6 +39,9 @@
"""
ElasticMaterialApp.__init__(self, name)
+ # import pdb
+ # pdb.set_trace()
+
numLocs = 2
self.dimension = dimension
@@ -72,7 +75,7 @@
viscosityA = 1.0e18
strainA = [1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4]
initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
- #initialStrainA = [3.6e-4, 3.5e-4, 3.4e-4, 3.3e-4, 3.2e-4, 3.1e-4]
+ initialStrainA = [3.6e-5, 3.5e-5, 3.4e-5, 3.3e-5, 3.2e-5, 3.1e-5]
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
maxwellTimeA = viscosityA / muA
@@ -85,15 +88,15 @@
strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
strainB = [4.1e-4, 4.2e-4, 4.3e-4, 4.4e-4, 4.5e-4, 4.6e-4]
initialStressB = [5.1e4, 5.2e4, 5.3e4, 5.4e4, 5.5e4, 5.6e4]
- #initialStrainB = [6.1e-4, 6.2e-4, 6.3e-4, 6.6e-4, 6.5e-4, 6.4e-4]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.3e-5, 6.6e-5, 6.5e-5, 6.4e-5]
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
maxwellTimeB = viscosityB / 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]
+ # 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]
diag = numpy.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0],
dtype=numpy.float64)
@@ -118,8 +121,10 @@
density0 = self.densityScale
time0 = self.timeScale
self.propertiesNondim = \
- numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, maxwellTimeA/time0],
- [densityB/density0, muB/mu0, lambdaB/mu0, maxwellTimeB/time0] ],
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, \
+ maxwellTimeA/time0],
+ [densityB/density0, muB/mu0, lambdaB/mu0, \
+ maxwellTimeB/time0] ],
dtype=numpy.float64)
self.initialStress = numpy.array([initialStressA,
@@ -139,8 +144,8 @@
# (viscous_strain) are defined by the assigned elastic strain.
totalStrainA = strainA
totalStrainB = strainB
- viscousStrainA = numpy.array(strainA) - diag*meanStrainA
- viscousStrainB = numpy.array(strainB) - diag*meanStrainB
+ viscousStrainA = numpy.array(strainA) - diag * meanStrainA
+ viscousStrainB = numpy.array(strainB) - diag * meanStrainB
self.stateVars = numpy.array([ [totalStrainA, viscousStrainA],
[totalStrainB, viscousStrainB] ],
dtype=numpy.float64)
@@ -149,15 +154,17 @@
self.strain = numpy.array([strainA, strainB],
dtype=numpy.float64)
self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
- self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
+ self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + tensorSize),
+ dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts),
dtype=numpy.float64)
- (self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA,
- muA, lambdaA, maxwellTimeA,
- totalStrainA, viscousStrainA,
- initialStressA, initialStrainA)
- (self.elasticConsts[1,:], self.stress[1,:]) = \
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA,
+ muA, lambdaA, maxwellTimeA,
+ totalStrainA, viscousStrainA,
+ initialStressA, initialStrainA)
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
self._calcStress(strainB,
muB, lambdaB, maxwellTimeB,
totalStrainB, viscousStrainB,
@@ -167,31 +174,83 @@
return
- def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, totalStrainV,
- viscousStrainV, initialStressV, initialStrainV):
+ def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
+ muV, lambdaV, maxwellTimeV, dqV, totalStrainT,
+ viscousStrainT, initialStress, initialStrain):
"""
- Compute stress and derivative of elasticity matrix.
- This assumes behavior is always viscoelastic.
+ 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,
+ maxwellTimeV,
+ dqV,
+ totalStrainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain)
+ return stressTpdt[stressComp]
+
+
+ def _computeStress(self, strainTpdt, muV, lambdaV, maxwellTimeV, dqV,
+ strainT, viscousStrainT,
+ initialStress, initialStrain):
+ """
+ Function to compute stresses and viscous strains for a given strain.
+ """
import math
- bulkModulus = lambdaV + 2.0 * muV/3.0
+ bulkModulus = lambdaV + 2.0 * muV / 3.0
diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
- totalStrainR = numpy.array(totalStrainV) - numpy.array(initialStrainV)
- print strainV
- print initialStrainV
- print totalStrainR
+ # Initial stresses and strains
+ meanStrainInitial = \
+ (initialStrain[0] + initialStrain[1] + initialStrain[2]) / 3.0
+ meanStressInitial = \
+ (initialStress[0] + initialStress[1] + initialStress[2]) / 3.0
- traceStrainT = totalStrainR[0] + totalStrainR[1] + totalStrainR[2]
- traceStrainTpdt = strainV[0] + strainV[1] + strainV[2]
- meanStrainT = traceStrainT / 3.0
- meanStrainTpdt = traceStrainTpdt / 3.0
- meanStressTpdt = bulkModulus * traceStrainTpdt
+ devStrainInitial = initialStrain - numpy.array(diag) * meanStrainInitial
+ devStressInitial = initialStress - numpy.array(diag) * meanStressInitial
+
+ meanStrainT = (strainT[0] + strainT[1] + strainT[2]) / 3.0
+ meanStrainTpdt = (strainTpdt[0] + strainTpdt[1] + strainTpdt[2]) / 3.0
+ meanStressTpdt = 3.0 * bulkModulus * \
+ (meanStrainTpdt - meanStrainInitial) + meanStressInitial
+
+ stressTpdt = numpy.zeros( (tensorSize), dtype=numpy.float64)
+ viscousStrainTpdt = numpy.zeros( (tensorSize), dtype=numpy.float64)
+
+ expFac = math.exp(-self.dt/maxwellTimeV)
+ elasFac = 2.0 * muV
+ devStrainTpdt = 0.0
+ devStrainT = 0.0
+ devStressTpdt = 0.0
+ for iComp in range(tensorSize):
+ devStrainTpdt = strainTpdt[iComp] - diag[iComp] * meanStrainTpdt
+ devStrainT = strainT[iComp] - diag[iComp] * meanStrainT
+ viscousStrainTpdt[iComp] = expFac * viscousStrainT[iComp] + \
+ dqV * (devStrainTpdt - devStrainT)
+ devStressTpdt = elasFac * \
+ (viscousStrainTpdt[iComp] - devStrainInitial[iComp]) + \
+ devStressInitial[iComp]
+ stressTpdt[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt
+
+ return stressTpdt, viscousStrainTpdt
+
+
+ def _computeViscousFactor(self, maxwellTime):
+ """
+ Compute viscous strain factor for a given Maxwell time.
+ """
+ import math
+
timeFrac = 1.0e-5
numTerms = 5
dq = 0.0
- if maxwellTimeV < timeFrac*self.dt:
+ if maxwellTime < timeFrac*self.dt:
fSign = 1.0
factorial = 1.0
fraction = 1.0
@@ -199,61 +258,75 @@
for iTerm in range(2, numTerms + 1):
factorial *= iTerm
fSign *= -1.0
- fraction *= self.dt/maxwellTimeV
+ fraction *= self.dt/maxwellTime
dq += fSign*fraction/factorial
else:
- dq = maxwellTimeV*(1.0-math.exp(-self.dt/maxwellTimeV))/self.dt
+ dq = maxwellTime*(1.0-math.exp(-self.dt/maxwellTime))/self.dt
- visFac = muV*dq/3.0
+ return dq
- C1111 = bulkModulus + 4.0*visFac
- C1122 = bulkModulus - 2.0*visFac
- C1133 = C1122
- C1112 = 0.0
- C1123 = 0.0
- C1113 = 0.0
- C2222 = C1111
- C2233 = C1122
- C2212 = 0.0
- C2223 = 0.0
- C2213 = 0.0
- C3333 = C1111
- C3312 = 0.0
- C3323 = 0.0
- C3313 = 0.0
- C1212 = 6.0*visFac
- C1223 = 0.0
- C1213 = 0.0
- C2323 = C1212
- C2313 = 0.0
- C1313 = C1212
- elasticConsts = numpy.array([C1111, C1122, C1133, C1112, C1123, C1113,
- C2222, C2233, C2212, C2223, C2213,
- C3333, C3312, C3323, C3313,
- C1212, C1223, C1213,
- C2323, C2313,
- C1313], dtype=numpy.float64)
+
+ def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, totalStrainV,
+ viscousStrainV, initialStressV, initialStrainV):
+ """
+ Compute stress, derivative of elasticity matrix, and updated state
+ variables. This assumes behavior is always viscoelastic.
+ """
+ import scipy.misc
- stressV = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ # Define some numpy arrays
+ strainTpdt = numpy.array(strainV, dtype=numpy.float64)
+ strainT = numpy.array(totalStrainV, dtype=numpy.float64)
+ viscousStrainT = numpy.array(viscousStrainV, dtype=numpy.float64)
+ initialStress = numpy.array(initialStressV, dtype=numpy.float64)
+ initialStrain = numpy.array(initialStrainV, dtype=numpy.float64)
+
+ # Compute current stress and viscous strain.
+ dqV = self._computeViscousFactor(maxwellTimeV)
+ stressTpdt, viscousStrainTpdt = self._computeStress(strainTpdt,
+ muV,
+ lambdaV,
+ maxwellTimeV,
+ dqV,
+ strainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain)
+
+ # Form updated state variables
+ stateVarsUpdated = numpy.array( [strainTpdt, viscousStrainTpdt],
+ dtype=numpy.float64)
- expFac = math.exp(-self.dt/maxwellTimeV)
- print "expFac:",expFac
- print "viscousStrain",viscousStrainV
- elasFac = 2.0*muV
- devStrainTpdt = 0.0
- devStrainT = 0.0
- devStressTpdt = 0.0
- viscousStrain = 0.0
- for iComp in range(tensorSize):
- devStrainTpdt = strainV[iComp] - diag[iComp]*meanStrainTpdt
- devStrainT = totalStrainR[iComp] - diag[iComp]*meanStrainT
- viscousStrain = expFac*viscousStrainV[iComp] + dq*(devStrainTpdt - devStrainT)
- devStressTpdt = elasFac*viscousStrain
- stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt + \
- initialStressV[iComp]
-
- stress = numpy.reshape(stressV, (tensorSize,1))
- return (elasticConsts, numpy.ravel(stress))
+ # Compute components of tangent constitutive matrix using numerical
+ # derivatives.
+ derivDx = 1.0e-12
+ derivOrder = 3
+ elasticConstsList = []
+
+ for stressComp in range(tensorSize):
+ for strainComp in range(stressComp, tensorSize):
+ dStressDStrain = stressComp + strainComp
+ dStressDStrain = scipy.misc.derivative(self._calcStressComponent,
+ strainTpdt[strainComp],
+ dx=derivDx,
+ args=(strainComp,
+ stressComp,
+ strainTpdt,
+ muV,
+ lambdaV,
+ maxwellTimeV,
+ dqV,
+ strainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain),
+ order=derivOrder)
+ elasticConstsList.append(dStressDStrain)
+
+ elasticConsts = numpy.array(elasticConstsList, dtype=numpy.float64)
+
+ return (elasticConsts, numpy.ravel(stressTpdt),
+ numpy.ravel(stateVarsUpdated))
# MAIN /////////////////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2009-07-06 02:30:18 UTC (rev 15429)
@@ -198,63 +198,63 @@
};
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_stress[] = {
- 1.30730205e+07,
- 1.35220000e+07,
- 1.39709795e+07,
- 6.29571369e+06,
- 6.74469324e+06,
- 7.19367279e+06,
- 6.04140332e+06,
- 6.10000000e+06,
- 6.15859668e+06,
- 2.58825402e+06,
- 2.64685071e+06,
- 2.70544739e+06,
+ 9.09052045e+06,
+ 9.58450000e+06,
+ 1.00784795e+07,
+ 4.81071369e+06,
+ 5.30469324e+06,
+ 5.79867279e+06,
+ 5.15436332e+06,
+ 5.20720000e+06,
+ 5.26003668e+06,
+ 2.20809402e+06,
+ 2.27245071e+06,
+ 2.33680739e+06,
};
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_elasticConsts[] = {
- 6.74326011e+10,
- 2.25336994e+10,
- 2.25336994e+10,
+ 6.74326010e+10,
+ 2.25336999e+10,
+ 2.25336999e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 6.74326011e+10,
- 2.25336994e+10,
+ 6.74326010e+10,
+ 2.25336999e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 6.74326011e+10,
+ 6.74326010e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 4.48989017e+10,
+ 4.48989016e+10,
0.00000000e+00,
0.00000000e+00,
- 4.48989017e+10,
+ 4.48989016e+10,
0.00000000e+00,
- 4.48989017e+10,
- 8.63988941e+09,
- 2.88005529e+09,
- 2.88005529e+09,
+ 4.48989011e+10,
+ 8.63988977e+09,
+ 2.88005546e+09,
+ 2.88005546e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 8.63988941e+09,
- 2.88005529e+09,
+ 8.63988977e+09,
+ 2.88005546e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 8.63988941e+09,
+ 8.63988930e+09,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 5.75983412e+09,
+ 5.75983408e+09,
0.00000000e+00,
0.00000000e+00,
- 5.75983412e+09,
+ 5.75983408e+09,
0.00000000e+00,
- 5.75983412e+09,
+ 5.75983408e+09,
};
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_initialStress[] = {
@@ -273,22 +273,47 @@
};
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_initialStrain[] = {
+ 3.60000000e-05,
+ 3.50000000e-05,
+ 3.40000000e-05,
+ 3.30000000e-05,
+ 3.20000000e-05,
+ 3.10000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.30000000e-05,
+ 6.60000000e-05,
+ 6.50000000e-05,
+ 6.40000000e-05,
+};
+
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_stateVarsUpdated[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.30000000e-04,
+ 1.40000000e-04,
+ 1.50000000e-04,
+ 1.60000000e-04,
+ -9.95510110e-06,
+ 1.34916778e-20,
+ 9.95510110e-06,
+ 1.39371415e-04,
+ 1.49326516e-04,
+ 1.59281618e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.30000000e-04,
+ 4.40000000e-04,
+ 4.50000000e-04,
+ 4.60000000e-04,
+ -9.99942402e-06,
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,
+ 9.99942402e-06,
+ 4.39974657e-04,
+ 4.49974081e-04,
+ 4.59973505e-04,
};
-const double* pylith::materials::MaxwellIsotropic3DTimeDepData::_stateVarsUpdated = 0;
-
pylith::materials::MaxwellIsotropic3DTimeDepData::MaxwellIsotropic3DTimeDepData(void)
{ // constructor
dimension = _dimension;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2009-07-06 02:30:18 UTC (rev 15429)
@@ -95,7 +95,7 @@
static const double _initialStrain[];
- static const double* _stateVarsUpdated;
+ static const double _stateVarsUpdated[];
};
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py 2009-07-06 02:30:18 UTC (rev 15429)
@@ -79,7 +79,7 @@
powerLawExpA = 1.0
strainA = [1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4]
initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
- initialStrainA = [3.6e-4, 3.5e-4, 3.4e-4, 3.3e-4, 3.2e-4, 3.1e-4]
+ initialStrainA = [3.6e-5, 3.5e-5, 3.4e-5, 3.3e-5, 3.2e-5, 3.1e-5]
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
@@ -90,13 +90,13 @@
powerLawExpB = 3.0
strainB = [4.1e-4, 4.2e-4, 4.3e-4, 4.4e-4, 4.5e-4, 4.6e-4]
initialStressB = [5.1e4, 5.2e4, 5.3e4, 5.4e4, 5.5e4, 5.6e4]
- initialStrainB = [6.1e-4, 6.2e-4, 6.3e-4, 6.6e-4, 6.5e-4, 6.4e-4]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.3e-5, 6.6e-5, 6.5e-5, 6.4e-5]
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
# 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]
+ # 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]
self.lengthScale = 1.0e+3
self.pressureScale = muA
@@ -351,7 +351,8 @@
initialStrain[2])/3.0
# Values for current time step
- meanStrainTpdt = (strainTpdt[0] + strainTpdt[1] + strainTpdt[2])/3.0
+ meanStrainTpdt = (strainTpdt[0] + strainTpdt[1] + strainTpdt[2])/3.0 - \
+ meanStrainInitial
meanStressTpdt = 3.0 * bulkModulus * meanStrainTpdt
strainPPTpdt = strainTpdt - meanStrainTpdt * diag - \
visStrainT - initialStrain
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc 2009-07-06 02:12:35 UTC (rev 15428)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc 2009-07-06 02:30:18 UTC (rev 15429)
@@ -206,33 +206,33 @@
};
const double pylith::materials::PowerLaw3DTimeDepData::_stress[] = {
- 1.12311566e+07,
- 1.16362430e+07,
- 1.20413293e+07,
- 4.33417161e+06,
- 4.73925792e+06,
- 5.14434423e+06,
- 5.97804000e+06,
- 6.03088000e+06,
- 6.08372000e+06,
- 2.50776000e+06,
- 2.56060000e+06,
- 2.61344000e+06,
+ 7.24875767e+06,
+ 7.69874295e+06,
+ 8.14872824e+06,
+ 2.85250536e+06,
+ 3.30249065e+06,
+ 3.75247593e+06,
+ 5.09100000e+06,
+ 5.13808000e+06,
+ 5.18516000e+06,
+ 2.12760000e+06,
+ 2.18620000e+06,
+ 2.24480000e+06,
};
const double pylith::materials::PowerLaw3DTimeDepData::_elasticConsts[] = {
- 6.74326513e+10,
- 2.25336738e+10,
- 2.25336738e+10,
+ 6.74326518e+10,
+ 2.25336747e+10,
+ 2.25336747e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 6.74326513e+10,
- 2.25336738e+10,
+ 6.74326518e+10,
+ 2.25336747e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 6.74326522e+10,
+ 6.74326518e+10,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
@@ -281,45 +281,45 @@
};
const double pylith::materials::PowerLaw3DTimeDepData::_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.60000000e-05,
+ 3.50000000e-05,
+ 3.40000000e-05,
+ 3.30000000e-05,
+ 3.20000000e-05,
+ 3.10000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.30000000e-05,
+ 6.60000000e-05,
+ 6.50000000e-05,
+ 6.40000000e-05,
};
const double pylith::materials::PowerLaw3DTimeDepData::_stateVarsUpdated[] = {
- 4.08854078e-05,
+ 4.08831629e-05,
4.19057121e-05,
- 4.29260165e-05,
- 4.42184086e-05,
- 4.52387129e-05,
- 4.62590172e-05,
- 1.12311566e+07,
- 1.16362430e+07,
- 1.20413293e+07,
- 4.33417161e+06,
- 4.73925792e+06,
- 5.14434423e+06,
+ 4.29282614e-05,
+ 4.41443253e-05,
+ 4.51668745e-05,
+ 4.61894238e-05,
+ 7.24875767e+06,
+ 7.69874295e+06,
+ 8.14872824e+06,
+ 2.85250536e+06,
+ 3.30249065e+06,
+ 3.75247593e+06,
1.10000000e-05,
1.20000000e-05,
1.30000000e-05,
- 1.40000007e-05,
- 1.50000007e-05,
- 1.60000007e-05,
- 5.97804000e+06,
- 6.03088000e+06,
- 6.08372000e+06,
- 2.50776000e+06,
- 2.56060000e+06,
- 2.61344000e+06,
+ 1.40000004e-05,
+ 1.50000004e-05,
+ 1.60000004e-05,
+ 5.09100000e+06,
+ 5.13808000e+06,
+ 5.18516000e+06,
+ 2.12760000e+06,
+ 2.18620000e+06,
+ 2.24480000e+06,
};
pylith::materials::PowerLaw3DTimeDepData::PowerLaw3DTimeDepData(void)
More information about the CIG-COMMITS
mailing list