[cig-commits] r15221 - in short/3D/PyLith/trunk: . libsrc libsrc/materials modulesrc/materials pylith/materials unittests/pytests/materials
brad at geodynamics.org
brad at geodynamics.org
Fri Jun 12 16:35:47 PDT 2009
Author: brad
Date: 2009-06-12 16:35:46 -0700 (Fri, 12 Jun 2009)
New Revision: 15221
Added:
short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellIsotropic3D.i
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
short/3D/PyLith/trunk/modulesrc/materials/Makefile.am
short/3D/PyLith/trunk/modulesrc/materials/materials.i
short/3D/PyLith/trunk/pylith/materials/GenMaxwellIsotropic3D.py
short/3D/PyLith/trunk/unittests/pytests/materials/TestGenMaxwellIsotropic3D.py
short/3D/PyLith/trunk/unittests/pytests/materials/testmaterials.py
Log:
Updated generalized Maxwell model (had not been updated to match new interface for materials).
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/TODO 2009-06-12 23:35:46 UTC (rev 15221)
@@ -2,6 +2,11 @@
CURRENT ISSUES/PRIORITIES
======================================================================
+CRITICAL
+ petsc-dev doesn't build on the cygwin buildbot (`_vsnprintf' undeclared)
+ http://www.geodynamics.org:8009/?project=PETSc
+ Generalized maxwell model
+
BUGS
malloc during assembly when running in parallel
@@ -43,9 +48,6 @@
Charles
3-D Power-law rheology
- libtests
- Python
- pytests
full-scale test (axial compression/extension?)
Generalized Maxwell materials (initial stress/strain)
CURRENTLY DISABLED!!!!
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-12 23:35:46 UTC (rev 15221)
@@ -80,6 +80,7 @@
materials/ElasticPlaneStress.cc \
materials/ElasticIsotropic3D.cc \
materials/ViscoelasticMaxwell.cc \
+ materials/GenMaxwellIsotropic3D.cc \
materials/MaxwellIsotropic3D.cc \
materials/PowerLaw3D.cc \
meshio/BinaryIO.cc \
@@ -109,7 +110,6 @@
utils/EventLogger.cc \
utils/TestArray.cc
-# materials/GenMaxwellIsotropic3D.cc \
# topology/MeshRefiner.cc \
# topology/RefineUniform.cc
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2009-06-12 23:35:46 UTC (rev 15221)
@@ -15,6 +15,7 @@
#include "GenMaxwellIsotropic3D.hh" // implementation of object methods
#include "ViscoelasticMaxwell.hh" // USES computeVisStrain
+#include "Metadata.hh" // USES Metadata
#include "pylith/utils/array.hh" // USES double_array
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
@@ -33,88 +34,168 @@
namespace materials {
namespace _GenMaxwellIsotropic3D{
+ // Dimension of material.
+ const int dimension = 3;
+
/// Number of entries in stress/strain tensors.
const int tensorSize = 6;
/// Number of entries in derivative of elasticity matrix.
const int numElasticConsts = 21;
+ /// Number of physical properties.
+ const int numProperties = 7;
+
/// Number of Maxwell models in parallel.
const int numMaxwellModels = 3;
- /// Number of physical properties.
- const int numProperties = 7;
-
/// Physical properties.
- const Material::PropMetaData properties[] = {
- { "density", 1, OTHER_FIELD },
- { "mu", 1, OTHER_FIELD },
- { "lambda", 1, OTHER_FIELD },
- { "shear_ratio", numMaxwellModels, OTHER_FIELD },
- { "maxwell_time", numMaxwellModels, OTHER_FIELD },
- { "total_strain", tensorSize, OTHER_FIELD },
- { "viscous_strain", numMaxwellModels*tensorSize, OTHER_FIELD },
+ const Metadata::ParamDescription properties[] = {
+ { "density", 1, pylith::topology::FieldBase::SCALAR },
+ { "mu", 1, pylith::topology::FieldBase::SCALAR },
+ { "lambda", 1, pylith::topology::FieldBase::SCALAR },
+ { "shear_ratio", numMaxwellModels, pylith::topology::FieldBase::OTHER },
+ { "maxwell_time", numMaxwellModels, pylith::topology::FieldBase::OTHER },
};
- /// Indices (order) of properties.
- const int pidDensity = 0;
- const int pidMuTot = pidDensity + 1;
- const int pidLambdaTot = pidMuTot + 1;
- const int pidShearRatio = pidLambdaTot + 1;
- const int pidMaxwellTime = pidShearRatio + numMaxwellModels;
- const int pidStrainT = pidMaxwellTime + numMaxwellModels;
- const int pidVisStrain = pidStrainT + tensorSize;
+ // Values expected in properties spatial database. :KLUDGE: Not
+ // generalized over number of models.
+ const int numDBProperties = 3 + 2*numMaxwellModels;
+ const char* dbProperties[] = {
+ "density", "vs", "vp" ,
+ "shear_ratio_1",
+ "shear_ratio_2",
+ "shear_ratio_3",
+ "viscosity_1",
+ "viscosity_2",
+ "viscosity_3",
+ };
+
+ /// Number of state variables.
+ const int numStateVars = 1+numMaxwellModels;
+
+ /// State variables. :KLUDGE: Not generalized over number of models.
+ const Metadata::ParamDescription stateVars[] = {
+ { "total_strain", tensorSize, pylith::topology::FieldBase::TENSOR },
+ { "viscous_strain_1", tensorSize, pylith::topology::FieldBase::TENSOR },
+ { "viscous_strain_2", tensorSize, pylith::topology::FieldBase::TENSOR },
+ { "viscous_strain_3", tensorSize, pylith::topology::FieldBase::TENSOR },
+ };
- /// Values expected in spatial database
- const int numDBValues = 3 + 2*numMaxwellModels;
+ // Values expected in state variables spatial database
+ const int numDBStateVars = tensorSize + numMaxwellModels*tensorSize;
+ const char* dbStateVars[] = {"total-strain-xx",
+ "total-strain-yy",
+ "total-strain-zz",
+ "total-strain-xy",
+ "total-strain-yz",
+ "total-strain-xz",
+ "viscous-strain1-xx",
+ "viscous-strain1-yy",
+ "viscous-strain1-zz",
+ "viscous-strain1-xy",
+ "viscous-strain1-yz",
+ "viscous-strain1-xz",
+ "viscous-strain2-xx",
+ "viscous-strain2-yy",
+ "viscous-strain2-zz",
+ "viscous-strain2-xy",
+ "viscous-strain2-yz",
+ "viscous-strain2-xz",
+ "viscous-strain3-xx",
+ "viscous-strain3-yy",
+ "viscous-strain3-zz",
+ "viscous-strain3-xy",
+ "viscous-strain3-yz",
+ "viscous-strain3-xz",
+ };
- /// NOTE: I haven't generalized the database names to the number
- /// of Maxwell models like I have for everything else.
- const char* namesDBValues[] = {"density", "vs", "vp" ,
- "shear_ratio_1",
- "shear_ratio_2",
- "shear_ratio_3",
- "viscosity_1",
- "viscosity_2",
- "viscosity_3"};
-
- /// Indices (order) of database values.
- const int didDensity = 0;
- const int didVs = 1;
- const int didVp = 2;
- const int didShearRatio1 = 3;
- const int didShearRatio2 = 4;
- const int didShearRatio3 = 5;
- const int didViscosity1 = 6;
- const int didViscosity2 = 7;
- const int didViscosity3 = 8;
-
- /// Initial state values expected in spatial database
- const int numInitialStateDBValues = tensorSize;
- const char* namesInitialStateDBValues[] = { "stress_xx", "stress_yy",
- "stress_zz", "stress_xy",
- "stress_yz", "stress_xz" };
-
} // _GenMaxwellIsotropic3D
} // materials
} // pylith
+// Indices of physical properties
+const int pylith::materials::GenMaxwellIsotropic3D::p_density = 0;
+
+const int pylith::materials::GenMaxwellIsotropic3D::p_muEff =
+ pylith::materials::GenMaxwellIsotropic3D::p_density + 1;
+
+const int pylith::materials::GenMaxwellIsotropic3D::p_lambdaEff =
+ pylith::materials::GenMaxwellIsotropic3D::p_muEff + 1;
+
+const int pylith::materials::GenMaxwellIsotropic3D::p_shearRatio =
+ pylith::materials::GenMaxwellIsotropic3D::p_lambdaEff + 1;
+
+const int pylith::materials::GenMaxwellIsotropic3D::p_maxwellTime =
+ pylith::materials::GenMaxwellIsotropic3D::p_shearRatio +
+ pylith::materials::_GenMaxwellIsotropic3D::numMaxwellModels;
+
+// Indices of database values (order must match dbProperties)
+const int pylith::materials::GenMaxwellIsotropic3D::db_density = 0;
+
+const int pylith::materials::GenMaxwellIsotropic3D::db_vs =
+ pylith::materials::GenMaxwellIsotropic3D::db_density + 1;
+
+const int pylith::materials::GenMaxwellIsotropic3D::db_vp =
+ pylith::materials::GenMaxwellIsotropic3D::db_vs + 1;
+
+const int pylith::materials::GenMaxwellIsotropic3D::db_shearRatio =
+ pylith::materials::GenMaxwellIsotropic3D::db_vp +
+ pylith::materials::_GenMaxwellIsotropic3D::numMaxwellModels;
+
+const int pylith::materials::GenMaxwellIsotropic3D::db_viscosity =
+ pylith::materials::GenMaxwellIsotropic3D::db_shearRatio +
+ pylith::materials::_GenMaxwellIsotropic3D::numMaxwellModels;
+
+// Indices of state variables
+const int pylith::materials::GenMaxwellIsotropic3D::s_totalStrain = 0;
+
+const int pylith::materials::GenMaxwellIsotropic3D::s_viscousStrain1 =
+ pylith::materials::GenMaxwellIsotropic3D::s_totalStrain +
+ pylith::materials::_GenMaxwellIsotropic3D::tensorSize;
+
+const int pylith::materials::GenMaxwellIsotropic3D::s_viscousStrain2 =
+ pylith::materials::GenMaxwellIsotropic3D::s_viscousStrain1 +
+ pylith::materials::_GenMaxwellIsotropic3D::tensorSize;
+
+const int pylith::materials::GenMaxwellIsotropic3D::s_viscousStrain3 =
+ pylith::materials::GenMaxwellIsotropic3D::s_viscousStrain2 +
+ pylith::materials::_GenMaxwellIsotropic3D::tensorSize;
+
+// Indices of database values (order must match dbStateVars)
+const int pylith::materials::GenMaxwellIsotropic3D::db_totalStrain = 0;
+
+const int pylith::materials::GenMaxwellIsotropic3D::db_viscousStrain1 =
+ pylith::materials::GenMaxwellIsotropic3D::db_totalStrain +
+ pylith::materials::_GenMaxwellIsotropic3D::tensorSize;
+
+const int pylith::materials::GenMaxwellIsotropic3D::db_viscousStrain2 =
+ pylith::materials::GenMaxwellIsotropic3D::db_viscousStrain1 +
+ pylith::materials::_GenMaxwellIsotropic3D::tensorSize;
+
+const int pylith::materials::GenMaxwellIsotropic3D::db_viscousStrain3 =
+ pylith::materials::GenMaxwellIsotropic3D::db_viscousStrain2 +
+ pylith::materials::_GenMaxwellIsotropic3D::tensorSize;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::GenMaxwellIsotropic3D::GenMaxwellIsotropic3D(void) :
- ElasticMaterial(_GenMaxwellIsotropic3D::tensorSize,
+ ElasticMaterial(_GenMaxwellIsotropic3D::dimension,
+ _GenMaxwellIsotropic3D::tensorSize,
_GenMaxwellIsotropic3D::numElasticConsts,
- _GenMaxwellIsotropic3D::namesDBValues,
- _GenMaxwellIsotropic3D::namesInitialStateDBValues,
- _GenMaxwellIsotropic3D::numDBValues,
- _GenMaxwellIsotropic3D::properties,
- _GenMaxwellIsotropic3D::numProperties),
- _calcElasticConstsFn(&pylith::materials::GenMaxwellIsotropic3D::_calcElasticConstsElastic),
- _calcStressFn(&pylith::materials::GenMaxwellIsotropic3D::_calcStressElastic),
- _updatePropertiesFn(&pylith::materials::GenMaxwellIsotropic3D::_updatePropertiesElastic)
+ Metadata(_GenMaxwellIsotropic3D::properties,
+ _GenMaxwellIsotropic3D::numProperties,
+ _GenMaxwellIsotropic3D::dbProperties,
+ _GenMaxwellIsotropic3D::numDBProperties,
+ _GenMaxwellIsotropic3D::stateVars,
+ _GenMaxwellIsotropic3D::numStateVars,
+ _GenMaxwellIsotropic3D::dbStateVars,
+ _GenMaxwellIsotropic3D::numDBStateVars)),
+ _calcElasticConstsFn(0),
+ _calcStressFn(0),
+ _updateStateVarsFn(0)
{ // constructor
- _dimension = 3;
- _visStrain.resize(_GenMaxwellIsotropic3D::numMaxwellModels *
- _GenMaxwellIsotropic3D::tensorSize);
+ useElasticBehavior(true);
+ _viscousStrain.resize(_GenMaxwellIsotropic3D::numMaxwellModels*_tensorSize);
} // constructor
// ----------------------------------------------------------------------
@@ -124,6 +205,29 @@
} // destructor
// ----------------------------------------------------------------------
+// Set whether elastic or inelastic constitutive relations are used.
+void
+pylith::materials::GenMaxwellIsotropic3D::useElasticBehavior(const bool flag)
+{ // useElasticBehavior
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::GenMaxwellIsotropic3D::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::GenMaxwellIsotropic3D::_calcElasticConstsElastic;
+ _updateStateVarsFn =
+ &pylith::materials::GenMaxwellIsotropic3D::_updateStateVarsElastic;
+
+ } else {
+ _calcStressFn =
+ &pylith::materials::GenMaxwellIsotropic3D::_calcStressViscoelastic;
+ _calcElasticConstsFn =
+ &pylith::materials::GenMaxwellIsotropic3D::_calcElasticConstsViscoelastic;
+ _updateStateVarsFn =
+ &pylith::materials::GenMaxwellIsotropic3D::_updateStateVarsViscoelastic;
+ } // if/else
+} // useElasticBehavior
+
+// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
pylith::materials::GenMaxwellIsotropic3D::_dbToProperties(
@@ -132,11 +236,11 @@
{ // _dbToProperties
assert(0 != propValues);
const int numDBValues = dbValues.size();
- assert(_GenMaxwellIsotropic3D::numDBValues == numDBValues);
+ assert(_GenMaxwellIsotropic3D::numDBProperties == numDBValues);
- const double density = dbValues[_GenMaxwellIsotropic3D::didDensity];
- const double vs = dbValues[_GenMaxwellIsotropic3D::didVs];
- const double vp = dbValues[_GenMaxwellIsotropic3D::didVp];
+ const double density = dbValues[db_density];
+ const double vs = dbValues[db_vs];
+ const double vp = dbValues[db_vp];
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
@@ -160,18 +264,30 @@
<< "vs: " << vs << "\n";
throw std::runtime_error(msg.str());
} // if
-
assert(mu > 0);
- propValues[_GenMaxwellIsotropic3D::pidDensity] = density;
- propValues[_GenMaxwellIsotropic3D::pidMuTot] = mu;
- propValues[_GenMaxwellIsotropic3D::pidLambdaTot] = lambda;
+ propValues[p_density] = density;
+ propValues[p_muEff] = mu;
+ propValues[p_lambdaEff] = lambda;
+ double visFrac = 0.0;
+ for (int imodel = 0; imodel < numMaxwellModels; ++imodel)
+ visFrac += propValues[db_shearRatio + imodel];
+ if (visFrac > 1.0) {
+ std::ostringstream msg;
+ msg << "Shear modulus ratios sum to a value greater than 1.0 for\n"
+ << "Generalized Maxwell model.\n"
+ << "Ratio 1: " << propValues[db_shearRatio ] << "\n"
+ << "Ratio 2: " << propValues[db_shearRatio+1] << "\n"
+ << "Ratio 3: " << propValues[db_shearRatio+2] << "\n"
+ << "Total: " << visFrac << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
// Loop over number of Maxwell models.
- for (int model =0;
- model < numMaxwellModels; ++model) {
- double muRatio = dbValues[_GenMaxwellIsotropic3D::didShearRatio1 + model];
- double viscosity = dbValues[_GenMaxwellIsotropic3D::didViscosity1 + model];
+ for (int imodel =0; imodel < numMaxwellModels; ++imodel) {
+ double muRatio = dbValues[db_shearRatio + imodel];
+ double viscosity = dbValues[db_viscosity + imodel];
double muFac = muRatio*mu;
double maxwellTime = 1.0e30;
if (muFac > 0.0)
@@ -185,8 +301,8 @@
<< "maxwellTime: " << maxwellTime << "\n";
throw std::runtime_error(msg.str());
} // if
- propValues[_GenMaxwellIsotropic3D::pidShearRatio + model] = muRatio;
- propValues[_GenMaxwellIsotropic3D::pidMaxwellTime + model] = maxwellTime;
+ propValues[p_shearRatio + imodel] = muRatio;
+ propValues[p_maxwellTime + imodel] = maxwellTime;
} // for
PetscLogFlops(6+2*numMaxwellModels);
@@ -200,24 +316,21 @@
{ // _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[_GenMaxwellIsotropic3D::pidDensity] =
- _normalizer->nondimensionalize(values[_GenMaxwellIsotropic3D::pidDensity],
- densityScale);
- values[_GenMaxwellIsotropic3D::pidMuTot] =
- _normalizer->nondimensionalize(values[_GenMaxwellIsotropic3D::pidMuTot],
- pressureScale);
- values[_GenMaxwellIsotropic3D::pidLambdaTot] =
- _normalizer->nondimensionalize(values[_GenMaxwellIsotropic3D::pidLambdaTot],
- pressureScale);
+ values[p_density] =
+ _normalizer->nondimensionalize(values[p_density], densityScale);
+ values[p_muEff] =
+ _normalizer->nondimensionalize(values[p_muEff], pressureScale);
+ values[p_lambdaEff] =
+ _normalizer->nondimensionalize(values[p_lambdaEff], pressureScale);
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
- _normalizer->nondimensionalize(&values[_GenMaxwellIsotropic3D::pidMaxwellTime],
+ _normalizer->nondimensionalize(&values[p_maxwellTime],
numMaxwellModels, timeScale);
-
+
PetscLogFlops(3+1*numMaxwellModels);
} // _nondimProperties
@@ -229,205 +342,115 @@
{ // _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[_GenMaxwellIsotropic3D::pidDensity] =
- _normalizer->dimensionalize(values[_GenMaxwellIsotropic3D::pidDensity],
- densityScale);
- values[_GenMaxwellIsotropic3D::pidMuTot] =
- _normalizer->dimensionalize(values[_GenMaxwellIsotropic3D::pidMuTot],
- pressureScale);
- values[_GenMaxwellIsotropic3D::pidLambdaTot] =
- _normalizer->dimensionalize(values[_GenMaxwellIsotropic3D::pidLambdaTot],
- pressureScale);
+ values[p_density] =
+ _normalizer->dimensionalize(values[p_density], densityScale);
+ values[p_muEff] =
+ _normalizer->dimensionalize(values[p_muEff], pressureScale);
+ values[p_lambdaEff] =
+ _normalizer->dimensionalize(values[p_lambdaEff], pressureScale);
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
- _normalizer->dimensionalize(&values[_GenMaxwellIsotropic3D::pidMaxwellTime],
+ _normalizer->dimensionalize(&values[p_maxwellTime],
numMaxwellModels, timeScale);
PetscLogFlops(3+1*numMaxwellModels);
} // _dimProperties
// ----------------------------------------------------------------------
-// Nondimensionalize initial state.
+// Compute initial state variables from values in spatial database.
void
-pylith::materials::GenMaxwellIsotropic3D::_nondimInitState(double* const values,
- const int nvalues) const
-{ // _nondimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _GenMaxwellIsotropic3D::numInitialStateDBValues);
+pylith::materials::GenMaxwellIsotropic3D::_dbToStateVars(
+ double* const stateValues,
+ const double_array& dbValues) const
+{ // _dbToStateVars
+ assert(0 != stateValues);
+ const int numDBValues = dbValues.size();
+ assert(_GenMaxwellIsotropic3D::numDBStateVars == numDBValues);
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->nondimensionalize(values, nvalues, pressureScale);
+ const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+ const int totalSize = (1 + numMaxwellModels) * _tensorSize;
+ assert(totalSize == _numVarsQuadPt);
+ assert(totalSize == numDBValues);
+ for (int i=0; i < totalSize; ++i)
+ stateValues[i] = dbValues[i];
- PetscLogFlops(nvalues);
-} // _nondimInitState
+ PetscLogFlops(0);
+} // _dbToStateVars
// ----------------------------------------------------------------------
-// Dimensionalize initial state.
-void
-pylith::materials::GenMaxwellIsotropic3D::_dimInitState(double* const values,
- const int nvalues) const
-{ // _dimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _GenMaxwellIsotropic3D::numInitialStateDBValues);
-
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->dimensionalize(values, nvalues, pressureScale);
-
- PetscLogFlops(nvalues);
-} // _dimInitState
-
-// ----------------------------------------------------------------------
// Compute density at location from properties.
void
-pylith::materials::GenMaxwellIsotropic3D::_calcDensity(
- double* const density,
- const double* properties,
- const int numProperties)
+pylith::materials::GenMaxwellIsotropic3D::_calcDensity(double* const density,
+ const double* properties,
+ 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[_GenMaxwellIsotropic3D::pidDensity];
+ density[0] = properties[p_density];
} // _calcDensity
// ----------------------------------------------------------------------
-// Compute viscous strain for current time step.
-void
-pylith::materials::GenMaxwellIsotropic3D::_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(_GenMaxwellIsotropic3D::tensorSize == strainSize);
- assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
-
- const int tensorSize = _GenMaxwellIsotropic3D::tensorSize;
- const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
-
- const double muRatio[] = {
- properties[_GenMaxwellIsotropic3D::pidShearRatio],
- properties[_GenMaxwellIsotropic3D::pidShearRatio+1],
- properties[_GenMaxwellIsotropic3D::pidShearRatio+2]};
- const double maxwellTime[] = {
- properties[_GenMaxwellIsotropic3D::pidMaxwellTime],
- properties[_GenMaxwellIsotropic3D::pidMaxwellTime+1],
- properties[_GenMaxwellIsotropic3D::pidMaxwellTime+2]};
-
- 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 diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
-
- const double meanStrainT =
- (properties[_GenMaxwellIsotropic3D::pidStrainT+0] +
- properties[_GenMaxwellIsotropic3D::pidStrainT+1] +
- properties[_GenMaxwellIsotropic3D::pidStrainT+2])/3.0;
-
- PetscLogFlops(6);
-
- // Compute Prony series terms
- double dq[] = { 0.0, 0.0, 0.0};
- for (int model=0; model < numMaxwellModels; ++model) {
- if (muRatio[model] != 0.0)
- dq[model] =
- ViscoelasticMaxwell::computeVisStrain(_dt, maxwellTime[model]);
- } // for
-
- // Compute new viscous strains
- double devStrainTpdt = 0.0;
- double devStrainT = 0.0;
- double deltaStrain = 0.0;
-
- for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
- devStrainT = properties[_GenMaxwellIsotropic3D::pidStrainT+iComp] -
- diag[iComp]*meanStrainT;
- deltaStrain = devStrainTpdt - devStrainT;
- PetscLogFlops(5);
- for (int model=0; model < numMaxwellModels; ++model) {
- if (muRatio[model] != 0.0) {
- int pindex = iComp + model * tensorSize;
- _visStrain[pindex] = exp(-_dt/maxwellTime[model])*
- properties[_GenMaxwellIsotropic3D::pidVisStrain + pindex] +
- dq[model] * deltaStrain;
- PetscLogFlops(6);
- } // if
- } // for
- } // for
-} // _computeStateVars
-
-// ----------------------------------------------------------------------
// Compute stress tensor at location from properties as an elastic
// material.
void
pylith::materials::GenMaxwellIsotropic3D::_calcStressElastic(
- double* const stress,
- const int stressSize,
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize,
- const bool computeStateVars)
+ double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars)
{ // _calcStressElastic
assert(0 != stress);
assert(_GenMaxwellIsotropic3D::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
- assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStrainSize);
- const double mu = properties[_GenMaxwellIsotropic3D::pidMuTot];
- const double lambda = properties[_GenMaxwellIsotropic3D::pidLambdaTot];
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
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 e23 = totalStrain[4];
- const double e13 = totalStrain[5];
+ // :TODO: Need to consider initial 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 e12 = totalStrain[3] - initialStrain[3];
+ const double e23 = totalStrain[4] - initialStrain[4];
+ const double e13 = totalStrain[5] - initialStrain[5];
- const double traceStrainTpdt = e11 + e22 + e33;
- const double s123 = lambda * traceStrainTpdt;
+ const double s123 = lambda * (e11 + e22 + e33);
- stress[0] = s123 + mu2*e11 + initialState[0];
- stress[1] = s123 + mu2*e22 + initialState[1];
- stress[2] = s123 + mu2*e33 + initialState[2];
- stress[3] = mu2 * e12 + initialState[3];
- stress[4] = mu2 * e23 + initialState[4];
- stress[5] = mu2 * e13 + initialState[5];
- // std::cout << " _calcStressElastic: " << std::endl;
- // std::cout << " totalStrain: " << std::endl;
- // for (int iComp=0; iComp < _GenMaxwellIsotropic3D::tensorSize; ++iComp)
- // std::cout << " " << totalStrain[iComp];
- // std::cout << std::endl << " stress: " << std::endl;
- // for (int iComp=0; iComp < _GenMaxwellIsotropic3D::tensorSize; ++iComp)
- // std::cout << " " << (*stress)[iComp];
- // std::cout << std::endl;
+ 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];
- PetscLogFlops(19);
+ PetscLogFlops(25);
} // _calcStressElastic
@@ -436,104 +459,102 @@
// material.
void
pylith::materials::GenMaxwellIsotropic3D::_calcStressViscoelastic(
- double* const stress,
- const int stressSize,
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize,
- const bool computeStateVars)
+ double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars)
{ // _calcStressViscoelastic
assert(0 != stress);
assert(_GenMaxwellIsotropic3D::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
- assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStrainSize);
- const int tensorSize = _GenMaxwellIsotropic3D::tensorSize;
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+ const int tensorSize = _tensorSize;
- const double mu = properties[_GenMaxwellIsotropic3D::pidMuTot];
- const double lambda = properties[_GenMaxwellIsotropic3D::pidLambdaTot];
- const double muRatio[] = {
- properties[_GenMaxwellIsotropic3D::pidShearRatio],
- properties[_GenMaxwellIsotropic3D::pidShearRatio+1],
- properties[_GenMaxwellIsotropic3D::pidShearRatio+2]};
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
+ const double muRatio[numMaxwellModels] = {
+ properties[p_shearRatio ],
+ properties[p_shearRatio+1],
+ properties[p_shearRatio+2]
+ };
const double mu2 = 2.0 * mu;
const double bulkModulus = lambda + mu2/3.0;
- 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];
+ // :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 traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
- const double meanStressTpdt = bulkModulus * traceStrainTpdt;
+ const double e123 = e11 + e22 + e33;
+ const double meanStrainTpdt = e123 / 3.0;
+ const double meanStressTpdt = bulkModulus * e123;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
double visFrac = 0.0;
- for (int model = 0; model < numMaxwellModels; ++model)
- visFrac += muRatio[model];
+ for (int imodel = 0; imodel < numMaxwellModels; ++imodel)
+ visFrac += muRatio[imodel];
+ assert(visFrac <= 1.0);
+ const double elasFrac = 1.0 - visFrac;
- if (visFrac > 1.0) {
- std::ostringstream msg;
- msg << "Shear modulus ratios sum to a value greater than 1.0 for\n"
- << "Generalized Maxwell model.\n"
- << "Ratio 1: " << muRatio[0] << "\n"
- << "Ratio 2: " << muRatio[1] << "\n"
- << "Ratio 3: " << muRatio[2] << "\n"
- << "Total: " << visFrac << "\n";
- throw std::runtime_error(msg.str());
- } // if
- double elasFrac = 1.0 - visFrac;
-
PetscLogFlops(8 + numMaxwellModels);
// Get viscous strains
if (computeStateVars) {
- pylith::materials::GenMaxwellIsotropic3D::_computeStateVars(properties,
- numProperties,
- totalStrain,
- strainSize,
- initialState,
- initialStateSize);
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
} else {
- memcpy(&_visStrain[0], &properties[_GenMaxwellIsotropic3D::pidVisStrain],
- numMaxwellModels * tensorSize * sizeof(double));
+ int index = 0;
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ _viscousStrain[index++] = stateVars[s_totalStrain+iComp];
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ _viscousStrain[index++] = stateVars[s_viscousStrain1+iComp];
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ _viscousStrain[index++] = stateVars[s_viscousStrain2+iComp];
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ _viscousStrain[index++] = stateVars[s_viscousStrain3+iComp];
} // else
-
// Compute new stresses
double devStrainTpdt = 0.0;
double devStressTpdt = 0.0;
- // std::cout << " _calcStressViscoelastic: " << std::endl;
- // std::cout << " stress totalStrain _visStrain: " << std::endl;
for (int iComp=0; iComp < tensorSize; ++iComp) {
devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
devStressTpdt = elasFrac*devStrainTpdt;
for (int model=0; model < numMaxwellModels; ++model) {
if (muRatio[model] != 0.0) {
int pindex = iComp + model * tensorSize;
- devStressTpdt += muRatio[model] * _visStrain[pindex];
+ devStressTpdt += muRatio[model] * _viscousStrain[pindex];
} // if
} // for
devStressTpdt = mu2*devStressTpdt;
stress[iComp] =diag[iComp] * meanStressTpdt + devStressTpdt +
- initialState[iComp];
- // std::cout << devStressTpdt << " " << stress[iComp] << std::endl;
-
- // Temporary to get stresses and strains.
- // std::cout << " " << stress[iComp] << " " << totalStrain[iComp] << " " << _visStrain << std:: endl;
+ initialStress[iComp];
} // for
PetscLogFlops((7 + 2*numMaxwellModels) * tensorSize);
@@ -543,25 +564,34 @@
// Compute derivative of elasticity matrix at location from properties.
void
pylith::materials::GenMaxwellIsotropic3D::_calcElasticConstsElastic(
- double* const elasticConsts,
- const int numElasticConsts,
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ 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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConstsElastic
assert(0 != elasticConsts);
assert(_GenMaxwellIsotropic3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
- assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStrainSize);
- const double mu = properties[_GenMaxwellIsotropic3D::pidMuTot];
- const double lambda = properties[_GenMaxwellIsotropic3D::pidLambdaTot];
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
@@ -596,42 +626,51 @@
// as a viscoelastic material.
void
pylith::materials::GenMaxwellIsotropic3D::_calcElasticConstsViscoelastic(
- double* const elasticConsts,
- const int numElasticConsts,
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ 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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConstsViscoelastic
assert(0 != elasticConsts);
assert(_GenMaxwellIsotropic3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
- assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
-
- const double mu = properties[_GenMaxwellIsotropic3D::pidMuTot];
- const double lambda = properties[_GenMaxwellIsotropic3D::pidLambdaTot];
+ assert(0 != initialStress);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStrainSize);
+
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
const double mu2 = 2.0 * mu;
- const double bulkModulus = lambda + mu2/3.0;
+ const double bulkModulus = lambda + mu2 / 3.0;
// Compute viscous contribution.
double visFac = 0.0;
double visFrac = 0.0;
double shearRatio = 0.0;
- for (int model = 0; model < numMaxwellModels; ++model) {
- shearRatio = properties[_GenMaxwellIsotropic3D::pidShearRatio + model];
+ for (int imodel = 0; imodel < numMaxwellModels; ++imodel) {
+ shearRatio = properties[p_shearRatio + imodel];
double maxwellTime = 1.0e30;
visFrac += shearRatio;
if (shearRatio != 0.0) {
- maxwellTime = properties[_GenMaxwellIsotropic3D::pidMaxwellTime + model];
+ maxwellTime = properties[p_maxwellTime + imodel];
visFac +=
- shearRatio*ViscoelasticMaxwell::computeVisStrain(_dt, maxwellTime);
+ shearRatio*ViscoelasticMaxwell::viscousStrainParam(_dt, maxwellTime);
} // if
} // for
double elasFrac = 1.0 - visFrac;
@@ -673,112 +712,247 @@
} // _calcElasticConstsViscoelastic
// ----------------------------------------------------------------------
-// Get stable time step for implicit time integration.
-double
-pylith::materials::GenMaxwellIsotropic3D::_stableTimeStepImplicit(const double* properties,
- const int numProperties) const
-{ // _stableTimeStepImplicit
- assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
-
- double dtStable = pylith::PYLITH_MAXDOUBLE;
-
- const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
- for (int i=0; i < numMaxwellModels; ++i) {
- const double maxwellTime =
- properties[_GenMaxwellIsotropic3D::pidMaxwellTime+i];
- const double dt = 0.1*maxwellTime;
- if (dt < dtStable)
- dtStable = dt;
- } // for
-
- return dtStable;
-} // _stableTimeStepImplicit
-
-// ----------------------------------------------------------------------
// Update state variables.
void
-pylith::materials::GenMaxwellIsotropic3D::_updatePropertiesElastic(
- double* const properties,
+pylith::materials::GenMaxwellIsotropic3D::_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(_GenMaxwellIsotropic3D::tensorSize == strainSize);
- assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStrainSize);
- const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+ 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 traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
-
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+ // Update total strain
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_totalStrain+iComp] = totalStrain[iComp];
+
// Initialize all viscous strains to deviatoric elastic strains.
double devStrain = 0.0;
double shearRatio = 0.0;
-
- for (int iComp=0; iComp < _GenMaxwellIsotropic3D::tensorSize; ++iComp) {
- properties[_GenMaxwellIsotropic3D::pidStrainT+iComp] = totalStrain[iComp];
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
devStrain = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
- for (int model = 0; model < numMaxwellModels; ++model) {
- shearRatio = properties[_GenMaxwellIsotropic3D::pidShearRatio + model];
- properties[_GenMaxwellIsotropic3D::pidVisStrain+
- iComp+model*_GenMaxwellIsotropic3D::tensorSize] =
- devStrain;
- } // for
+ // Maxwell model 1
+ shearRatio = properties[p_shearRatio + 0];
+ stateVars[s_viscousStrain1+iComp] = devStrain;
+ // Maxwell model 2
+ shearRatio = properties[p_shearRatio + 1];
+ stateVars[s_viscousStrain2+iComp] = devStrain;
+ // Maxwell model 3
+ shearRatio = properties[p_shearRatio + 2];
+ stateVars[s_viscousStrain3+iComp] = devStrain;
} // for
- PetscLogFlops(3 + 2 * _GenMaxwellIsotropic3D::tensorSize);
+ PetscLogFlops(3 + 2 * tensorSize);
_needNewJacobian = true;
-} // _updatePropertiesElastic
+} // _updateStateVarsElastic
// ----------------------------------------------------------------------
// Update state variables.
void
-pylith::materials::GenMaxwellIsotropic3D::_updatePropertiesViscoelastic(
- double* const properties,
+pylith::materials::GenMaxwellIsotropic3D::_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 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(_GenMaxwellIsotropic3D::tensorSize == strainSize);
- assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStrainSize);
- const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
- const int tensorSize = _GenMaxwellIsotropic3D::tensorSize;
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
- pylith::materials::GenMaxwellIsotropic3D::_computeStateVars(properties,
- numProperties,
- totalStrain,
- strainSize,
- initialState,
- initialStateSize);
+ const int tensorSize = _tensorSize;
- memcpy(&properties[_GenMaxwellIsotropic3D::pidVisStrain],
- &_visStrain[0],
- numMaxwellModels * tensorSize * sizeof(double));
- memcpy(&properties[_GenMaxwellIsotropic3D::pidStrainT],
- &totalStrain[0],
- tensorSize * sizeof(double));
+ // Total strain
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_totalStrain+iComp] = totalStrain[iComp];
+ // Viscous strain
+ int index = 0;
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_viscousStrain1+iComp] = _viscousStrain[index++];
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_viscousStrain2+iComp] = _viscousStrain[index++];
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_viscousStrain3+iComp] = _viscousStrain[index++];
+
_needNewJacobian = false;
+} // _updateStateVarsViscoelastic
-} // _updatePropertiesViscoelastic
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+pylith::materials::GenMaxwellIsotropic3D::_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);
+ double dtStable = pylith::PYLITH_MAXDOUBLE;
+ const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+ for (int i=0; i < numMaxwellModels; ++i) {
+ const double maxwellTime = properties[p_maxwellTime+i];
+ const double dt = 0.1*maxwellTime;
+ if (dt < dtStable)
+ dtStable = dt;
+ } // for
+
+ return dtStable;
+} // _stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Compute viscous strain for current time step.
+void
+pylith::materials::GenMaxwellIsotropic3D::_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(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStrainSize);
+
+ const int tensorSize = _tensorSize;
+ const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+
+ const double muRatio[numMaxwellModels] = {
+ properties[p_shearRatio ],
+ properties[p_shearRatio+1],
+ properties[p_shearRatio+2]
+ };
+ const double maxwellTime[numMaxwellModels] = {
+ properties[p_maxwellTime ],
+ properties[p_maxwellTime+1],
+ properties[p_maxwellTime+2]
+ };
+
+ // :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 diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
+ const double meanStrainT =
+ ( stateVars[s_totalStrain+0] +
+ stateVars[s_totalStrain+1] +
+ stateVars[s_totalStrain+2] ) / 3.0;
+
+ PetscLogFlops(6);
+
+ // Compute Prony series terms
+ double_array dq(numMaxwellModels);
+ dq = 0.0;
+ for (int i=0; i < numMaxwellModels; ++i)
+ if (muRatio[i] != 0.0)
+ dq[i] = ViscoelasticMaxwell::viscousStrainParam(_dt, maxwellTime[i]);
+
+ // Compute new viscous strains
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+ 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;
+ deltaStrain = devStrainTpdt - devStrainT;
+
+ // Maxwell model 1
+ int imodel = 0;
+ int vindex = imodel*tensorSize + iComp;
+ if (0.0 != muRatio[imodel]) {
+ _viscousStrain[vindex] = exp(-_dt/maxwellTime[imodel])*
+ stateVars[s_viscousStrain1 + iComp] + dq[imodel] * deltaStrain;
+ PetscLogFlops(6);
+ } // if
+
+ // Maxwell model 2
+ imodel = 1;
+ if (0.0 != muRatio[imodel]) {
+ _viscousStrain[vindex] = exp(-_dt/maxwellTime[imodel])*
+ stateVars[s_viscousStrain2 + iComp] + dq[imodel] * deltaStrain;
+ PetscLogFlops(6);
+ } // if
+
+ // Maxwell model 3
+ imodel = 2;
+ if (0.0 != muRatio[imodel]) {
+ _viscousStrain[vindex] = exp(-_dt/maxwellTime[imodel])*
+ stateVars[s_viscousStrain3 + iComp] + dq[imodel] * deltaStrain;
+ PetscLogFlops(6);
+ } // if
+
+ } // for
+
+ PetscLogFlops(5*tensorSize);
+} // _computeStateVars
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh 2009-06-12 23:35:46 UTC (rev 15221)
@@ -39,17 +39,10 @@
#if !defined(pylith_materials_genmaxwellisotropic3d_hh)
#define pylith_materials_genmaxwellisotropic3d_hh
-#include "ElasticMaterial.hh"
+// Include directives ---------------------------------------------------
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
-/// Namespace for pylith package
-namespace pylith {
- namespace materials {
- class GenMaxwellIsotropic3D;
- class TestGenMaxwellIsotropic3D; // unit testing
- } // materials
-} // pylith
-
-/// 3-D, isotropic, generalized linear Maxwell viscoelastic material.
+// GenMaxwellIsotropic3D ------------------------------------------------
class pylith::materials::GenMaxwellIsotropic3D : public ElasticMaterial
{ // class GenMaxwellIsotropic3D
friend class TestGenMaxwellIsotropic3D; // unit testing
@@ -75,13 +68,6 @@
*/
void useElasticBehavior(const bool flag);
- /** Get flag indicating whether material implements an empty
- * _updateProperties() method.
- *
- * @returns False if _updateProperties() is empty, true otherwise.
- */
- bool usesUpdateProperties(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -112,55 +98,63 @@
void _dimProperties(double* const values,
const int nvalues) const;
- /** Nondimensionalize initial state.
+ /** Compute initial state variables from values in spatial database.
*
- * @param values Array of initial state values.
- * @param nvalues Number of values.
+ * @param stateValues Array of state variable values.
+ * @param dbValues Array of database values.
*/
- void _nondimInitState(double* const values,
- const int nvalues) const;
+ void _dbToStateVars(double* const stateValues,
+ const double_array& dbValues) const;
- /** Dimensionalize initial state.
- *
- * @param values Array of initial state values.
- * @param nvalues Number of values.
- */
- void _dimInitState(double* const values,
- const int nvalues) const;
+ // Note: We do not need to dimensionalize or nondimensionalize state
+ // variables because there are strains, which are dimensionless.
+
/** Compute density from properties.
*
* @param density Array for density.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
*/
void _calcDensity(double* const density,
const double* properties,
- const int numProperties);
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars);
- /** Compute stress tensor from properties. If the state variables
- * are from the previous time step, then the computeStateVars flag
- * should be set to true so that the state variables are updated
- * (but not stored) when computing the stresses.
+ /** Compute stress tensor from properties and state variables. If
+ * the state variables are from the previous time step, then the
+ * computeStateVars flag should be set to true so that the state
+ * variables are updated (but not stored) when computing the stresses.
*
* @param stress Array for stress tensor.
* @param stressSize Size of stress tensor.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
- * @param computeStateVars Flag indicating to compute updated state vars.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
+ * @param computeStateVars Flag indicating to compute updated state variables.
*/
void _calcStress(double* const stress,
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);
/** Compute derivatives of elasticity matrix from properties.
@@ -169,42 +163,65 @@
* @param numElasticConsts Number of elastic constants.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
*/
void _calcElasticConsts(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);
- /** Get stable time step for implicit time integration.
+ /** Update state variables (for next time step).
*
- * @returns Time step
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
*/
- double _stableTimeStepImplicit(const double* properties,
- const int numProperties) const;
+ void _updateStateVars(double* const 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);
- /** Update properties (for next time step).
+ /** Get stable time step for implicit time integration.
*
* @param properties Properties at location.
* @param numProperties Number of properties.
- * @param totalStrain Total strain at location.
- * @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Time step
*/
- void _updateProperties(double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
+ double _stableTimeStepImplicit(const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const;
// PRIVATE TYPEDEFS ///////////////////////////////////////////////////
private :
@@ -219,6 +236,10 @@
const int,
const double*,
const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
const bool);
/// Member prototype for _calcElasticConsts()
@@ -230,56 +251,56 @@
const double*,
const int,
const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
const int);
- /// Member prototype for _updateProperties()
- typedef void (pylith::materials::GenMaxwellIsotropic3D::*updateProperties_fn_type)
+ /// Member prototype for _updateStateVars()
+ typedef void (pylith::materials::GenMaxwellIsotropic3D::*updateStateVars_fn_type)
(double* const,
const int,
const double*,
const int,
const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
const int);
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
- /** Compute viscous strains (state variables) for the current time step.
- *
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param totalStrain Total strain at location.
- * @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
- */
- void _computeStateVars(const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
-
/** Compute stress tensor from properties as an elastic material.
*
* @param stress Array for stress tensor.
* @param stressSize Size of stress tensor.
* @param properties Properties at locations.
* @param numProperties Number of properties.
+ * @param stateVars State variables at locations.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at locations.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
* @param computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStressElastic(double* const stress,
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);
/** Compute stress tensor from properties as an viscoelastic material.
@@ -288,20 +309,28 @@
* @param stressSize Size of stress tensor.
* @param properties Properties at locations.
* @param numProperties Number of properties.
+ * @param stateVars State variables at locations.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at locations.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
* @param computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStressViscoelastic(double* const stress,
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);
/** Compute derivatives of elasticity matrix from properties as an
@@ -311,19 +340,27 @@
* @param numElasticConsts Number of elastic constants.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars State variables at locations.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
*/
void _calcElasticConstsElastic(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);
/** Compute derivatives of elasticity matrix from properties as a
* viscoelastic material.
@@ -332,35 +369,51 @@
* @param numElasticConsts Number of elastic constants.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars State variables at locations.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
*/
void _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);
/** Update state variables after solve as an elastic material.
*
+ * @param stateVars State variables at locations.
+ * @param numStateVars Number of state variables.
* @param properties Properties at location.
* @param numProperties Number of properties.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
*/
- void _updatePropertiesElastic(double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
+ void _updateStateVarsElastic(double* const 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);
/** Update state variables after solve as a viscoelastic material.
*
@@ -371,27 +424,47 @@
* @param initialState Initial state values.
* @param initialStateSize Size of initial state array.
*/
- void _updatePropertiesViscoelastic(double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
+ void _updateStateVarsViscoelastic(double* const 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);
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
-private :
+ /** Compute viscous strains (state variables) for the current time
+ * step.
+ *
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _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);
- /// Not implemented
- GenMaxwellIsotropic3D(const GenMaxwellIsotropic3D& m);
-
- /// Not implemented
- const GenMaxwellIsotropic3D& operator=(const GenMaxwellIsotropic3D& m);
-
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
/// Viscous strain array.
- double_array _visStrain;
+ double_array _viscousStrain;
/// Method to use for _calcElasticConsts().
calcElasticConsts_fn_type _calcElasticConstsFn;
@@ -399,9 +472,41 @@
/// Method to use for _calcStress().
calcStress_fn_type _calcStressFn;
- /// Method to use for _updateProperties().
- updateProperties_fn_type _updatePropertiesFn;
+ /// Method to use for _updateStateVars().
+ updateStateVars_fn_type _updateStateVarsFn;
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ static const int p_density;
+ static const int p_muEff;
+ static const int p_lambdaEff;
+ static const int p_shearRatio;
+ static const int p_maxwellTime;
+ static const int db_density;
+ static const int db_vs;
+ static const int db_vp;
+ static const int db_shearRatio;
+ static const int db_viscosity;
+
+ static const int s_totalStrain;
+ static const int s_viscousStrain1;
+ static const int s_viscousStrain2;
+ static const int s_viscousStrain3;
+ static const int db_totalStrain;
+ static const int db_viscousStrain1;
+ static const int db_viscousStrain2;
+ static const int db_viscousStrain3;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ GenMaxwellIsotropic3D(const GenMaxwellIsotropic3D&);
+
+ /// Not implemented
+ const GenMaxwellIsotropic3D& operator=(const GenMaxwellIsotropic3D&);
+
}; // class GenMaxwellIsotropic3D
#include "GenMaxwellIsotropic3D.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc 2009-06-12 23:35:46 UTC (rev 15221)
@@ -27,52 +27,30 @@
_dt = dt;
} // timeStep
-// Set whether elastic or inelastic constitutive relations are used.
-inline
-void
-pylith::materials::GenMaxwellIsotropic3D::useElasticBehavior(const bool flag) {
- if (flag) {
- _calcStressFn =
- &pylith::materials::GenMaxwellIsotropic3D::_calcStressElastic;
- _calcElasticConstsFn =
- &pylith::materials::GenMaxwellIsotropic3D::_calcElasticConstsElastic;
- _updatePropertiesFn =
- &pylith::materials::GenMaxwellIsotropic3D::_updatePropertiesElastic;
- } else {
- _calcStressFn =
- &pylith::materials::GenMaxwellIsotropic3D::_calcStressViscoelastic;
- _calcElasticConstsFn =
- &pylith::materials::GenMaxwellIsotropic3D::_calcElasticConstsViscoelastic;
- _updatePropertiesFn =
- &pylith::materials::GenMaxwellIsotropic3D::_updatePropertiesViscoelastic;
- } // if/else
-} // useElasticBehavior
-
-// Get flag indicating whether material implements an empty
-inline
-bool
-pylith::materials::GenMaxwellIsotropic3D::usesUpdateProperties(void) const {
- return true;
-} // usesUpdateProperties
-
// Compute stress tensor from parameters.
inline
void
pylith::materials::GenMaxwellIsotropic3D::_calcStress(
- double* const stress,
- const int stressSize,
- const double* parameters,
- const int numParams,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize,
- const bool computeStateVars) {
+ double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars) {
assert(0 != _calcStressFn);
CALL_MEMBER_FN(*this, _calcStressFn)(stress, stressSize,
- parameters, numParams,
+ properties, numProperties,
+ stateVars, numStateVars,
totalStrain, strainSize,
- initialState, initialStateSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize,
computeStateVars);
} // _calcStress
@@ -80,35 +58,48 @@
inline
void
pylith::materials::GenMaxwellIsotropic3D::_calcElasticConsts(
- double* const elasticConsts,
- const int numElasticConsts,
- const double* parameters,
- const int numParams,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize) {
+ 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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize) {
assert(0 != _calcElasticConstsFn);
CALL_MEMBER_FN(*this, _calcElasticConstsFn)(elasticConsts, numElasticConsts,
- parameters, numParams,
+ properties, numProperties,
+ stateVars, numStateVars,
totalStrain, strainSize,
- initialState, initialStateSize);
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
} // _calcElasticConsts
// Update state variables after solve.
inline
void
-pylith::materials::GenMaxwellIsotropic3D::_updateProperties(double* const parameters,
- const int numParams,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize) {
- assert(0 != _updatePropertiesFn);
- CALL_MEMBER_FN(*this, _updatePropertiesFn)(parameters, numParams,
- totalStrain, strainSize,
- initialState, initialStateSize);
-} // _updateProperties
+pylith::materials::GenMaxwellIsotropic3D::_updateStateVars(
+ double* const 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) {
+ assert(0 != _updateStateVarsFn);
+ CALL_MEMBER_FN(*this, _updateStateVarsFn)(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+} // _updateStateVars
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2009-06-12 23:35:46 UTC (rev 15221)
@@ -310,93 +310,6 @@
} // _calcDensity
// ----------------------------------------------------------------------
-// Get stable time step for implicit time integration.
-double
-pylith::materials::MaxwellIsotropic3D::_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
-
-#include <iostream>
-// ----------------------------------------------------------------------
-// Compute viscous strain for current time step.
-// material.
-void
-pylith::materials::MaxwellIsotropic3D::_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(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(0 != initialStress);
- assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
- assert(0 != initialStrain);
- assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
-
- const int tensorSize = _tensorSize;
- const double maxwellTime = properties[p_maxwellTime];
-
- // :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 diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
-
- const double meanStrainT =
- ( stateVars[s_totalStrain+0] +
- stateVars[s_totalStrain+1] +
- stateVars[s_totalStrain+2] ) / 3.0;
-
- // Time integration.
- double dq = ViscoelasticMaxwell::viscousStrainParam(_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 = stateVars[s_totalStrain+iComp] - diag[iComp] * meanStrainT;
- _viscousStrain[iComp] = expFac * stateVars[s_viscousStrain+iComp] +
- dq * (devStrainTpdt - devStrainT);
- } // for
-
- PetscLogFlops(8 + 7 * tensorSize);
-} // _computeStateVars
-
-// ----------------------------------------------------------------------
// Compute stress tensor at location from properties as an elastic
// material.
void
@@ -748,5 +661,91 @@
_needNewJacobian = false;
} // _updateStateVarsViscoelastic
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+pylith::materials::MaxwellIsotropic3D::_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.
+// material.
+void
+pylith::materials::MaxwellIsotropic3D::_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(_MaxwellIsotropic3D::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
+
+ const int tensorSize = _tensorSize;
+ const double maxwellTime = properties[p_maxwellTime];
+
+ // :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 diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
+ const double meanStrainT =
+ ( stateVars[s_totalStrain+0] +
+ stateVars[s_totalStrain+1] +
+ stateVars[s_totalStrain+2] ) / 3.0;
+
+ // Time integration.
+ double dq = ViscoelasticMaxwell::viscousStrainParam(_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 = stateVars[s_totalStrain+iComp] - diag[iComp] * meanStrainT;
+ _viscousStrain[iComp] = expFac * stateVars[s_viscousStrain+iComp] +
+ dq * (devStrainTpdt - devStrainT);
+ } // for
+
+ PetscLogFlops(8 + 7 * tensorSize);
+} // _computeStateVars
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2009-06-12 23:35:46 UTC (rev 15221)
@@ -258,31 +258,6 @@
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
- /** Compute viscous strains (state variables) for the current time
- * step.
- *
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param totalStrain Total strain at location.
- * @param strainSize Size of strain tensor.
- * @param initialStress Initial stress tensor at location.
- * @param initialStressSize Size of initial stress array.
- * @param initialStrain Initial strain tensor at location.
- * @param initialStrainSize Size of initial strain array.
- */
- void _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);
-
/** Compute stress tensor from properties as an elastic material.
*
* @param stress Array for stress tensor.
@@ -449,6 +424,31 @@
const double* initialStrain,
const int initialStrainSize);
+ /** Compute viscous strains (state variables) for the current time
+ * step.
+ *
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _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);
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Added: short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellIsotropic3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellIsotropic3D.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellIsotropic3D.i 2009-06-12 23:35:46 UTC (rev 15221)
@@ -0,0 +1,203 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file modulesrc/materials/GenMaxwellIsotropic3D.ii
+ *
+ * @brief Python interface to C++ GenMaxwellIsotropic3D object.
+ */
+
+namespace pylith {
+ namespace materials {
+
+ class pylith::materials::GenMaxwellIsotropic3D : public ElasticMaterial
+ { // class GenMaxwellIsotropic3D
+
+ // PUBLIC METHODS /////////////////////////////////////////////////
+ public :
+
+ /// Default constructor
+ GenMaxwellIsotropic3D(void);
+
+ /// Destructor
+ ~GenMaxwellIsotropic3D(void);
+
+ /** Set current time step.
+ *
+ * @param dt Current time step.
+ */
+ void timeStep(const double dt);
+
+ /** Set whether elastic or inelastic constitutive relations are used.
+ *
+ * @param flag True to use elastic, false to use inelastic.
+ */
+ void useElasticBehavior(const bool flag);
+
+ // PROTECTED METHODS //////////////////////////////////////////////
+ protected :
+
+ /** Compute properties from values in spatial database.
+ *
+ * Order of values in arrays matches order used in dbValues() and
+ * parameterNames().
+ *
+ * @param propValues Array of property values.
+ * @param dbValues Array of database values.
+ */
+ void _dbToProperties(double* const propValues,
+ const double_array& dbValues) const;
+
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _nondimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _dimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Compute initial state variables from values in spatial database.
+ *
+ * @param stateValues Array of state variable values.
+ * @param dbValues Array of database values.
+ */
+ void _dbToStateVars(double* const stateValues,
+ const double_array& dbValues) const;
+
+ /** Compute density from properties.
+ *
+ * @param density Array for density.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ */
+ void _calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars);
+
+ /** Compute stress tensor from properties and state variables. If
+ * the state variables are from the previous time step, then the
+ * computeStateVars flag should be set to true so that the state
+ * variables are updated (but not stored) when computing the stresses.
+ *
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
+ * @param computeStateVars Flag indicating to compute updated
+ * state variables.
+ */
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars);
+
+ /** Compute derivatives of elasticity matrix from properties.
+ *
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _calcElasticConsts(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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** Update state variables (for next time step).
+ *
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _updateStateVars(double* const 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);
+
+ /** Get stable time step for implicit time integration.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Time step
+ */
+ double _stableTimeStepImplicit(const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const;
+
+ }; // class GenMaxwellIsotropic3D
+
+ } // materials
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2009-06-12 23:35:46 UTC (rev 15221)
@@ -27,6 +27,7 @@
ElasticPlaneStress.i \
ElasticIsotropic3D.i \
MaxwellIsotropic3D.i \
+ GenMaxwellIsotropic3D.i \
PowerLaw3D.i
swig_generated = \
Modified: short/3D/PyLith/trunk/modulesrc/materials/materials.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/materials.i 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/modulesrc/materials/materials.i 2009-06-12 23:35:46 UTC (rev 15221)
@@ -27,7 +27,7 @@
#include "pylith/materials/ElasticPlaneStress.hh"
#include "pylith/materials/ElasticIsotropic3D.hh"
#include "pylith/materials/MaxwellIsotropic3D.hh"
-//#include "pylith/materials/GenMaxwellIsotropic3D.hh"
+#include "pylith/materials/GenMaxwellIsotropic3D.hh"
#include "pylith/materials/PowerLaw3D.hh"
#include "pylith/utils/arrayfwd.hh"
@@ -63,7 +63,7 @@
%include "ElasticPlaneStress.i"
%include "ElasticIsotropic3D.i"
%include "MaxwellIsotropic3D.i"
-//%include "GenMaxwellIsotropic3D.i"
+%include "GenMaxwellIsotropic3D.i"
%include "PowerLaw3D.i"
Modified: short/3D/PyLith/trunk/pylith/materials/GenMaxwellIsotropic3D.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/GenMaxwellIsotropic3D.py 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/pylith/materials/GenMaxwellIsotropic3D.py 2009-06-12 23:35:46 UTC (rev 15221)
@@ -12,16 +12,19 @@
## @file pylith/materials/GenMaxwellIsotropic3D.py
##
-## @brief Python object implementing 3-D generalized isotropic linear Maxwell viscoelastic material.
+## @brief Python object implementing 3-D isotropic linear GenMaxwell
+## viscoelastic material.
##
## Factory: material.
from ElasticMaterial import ElasticMaterial
+from materials import GenMaxwellIsotropic3D as ModuleGenMaxwellIsotropic3D
# GenMaxwellIsotropic3D class
-class GenMaxwellIsotropic3D(ElasticMaterial):
+class GenMaxwellIsotropic3D(ElasticMaterial, ModuleGenMaxwellIsotropic3D):
"""
- Python object implementing 3-D generalized isotropic linear Maxwell viscoelastic material.
+ Python object implementing 3-D isotropic linear GenMaxwell viscoelastic
+ material.
Factory: material.
"""
@@ -38,20 +41,22 @@
{'info': [],
'data': []},
'cell': \
- {'info': ["mu", "lambda", "density", "shear_ratio", "maxwell_time"],
- 'data': ["total_strain", "viscous_strain", "stress"]}}
- self._loggingPrefix = "MaGM3D "
+ {'info': ["mu", "lambda", "density",
+ "shear_ratio", "maxwell_time"],
+ 'data': ["total_strain", "stress",
+ "viscous_strain_1",
+ "viscous_strain_2",
+ "viscous_strain_3",
+ ]}}
+ self._loggingPrefix = "MaMx3D "
return
- def _createCppHandle(self):
+ def _createModuleObj(self):
"""
- Create handle to corresponding C++ object.
+ Call constructor for module object for access to C++ object.
"""
- if None == self.cppHandle:
- import pylith.materials.materials as bindings
- self.cppHandle = bindings.GenMaxwellIsotropic3D()
- self.dimension = self.cppHandle.dimension
+ ModuleGenMaxwellIsotropic3D.__init__(self)
return
Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestGenMaxwellIsotropic3D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestGenMaxwellIsotropic3D.py 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestGenMaxwellIsotropic3D.py 2009-06-12 23:35:46 UTC (rev 15221)
@@ -24,37 +24,54 @@
Unit testing of GenMaxwellIsotropic3D object.
"""
+ def setUp(self):
+ """
+ Setup test subject.
+ """
+ self.material = GenMaxwellIsotropic3D()
+ return
+
+
def test_constructor(self):
"""
Test constructor.
"""
- from pylith.materials.GenMaxwellIsotropic3D import GenMaxwellIsotropic3D
- material = GenMaxwellIsotropic3D()
- material._createCppHandle()
- self.assertNotEqual(None, material.cppHandle)
+ self.assertEqual(3, self.material.dimension())
return
- def test_dimension(self):
+ def test_useElasticBehavior(self):
"""
- Test dimension().
+ Test useElasticBehavior().
"""
- material = GenMaxwellIsotropic3D()
- material._createCppHandle()
- self.assertEqual(3, material.dimension)
+ self.material.useElasticBehavior(False)
return
- def test_useElasticBehavior(self):
+ def testHasStateVars(self):
+ self.failUnless(self.material.hasStateVars())
+ return
+
+
+ def testTensorSize(self):
+ self.assertEqual(6, self.material.tensorSize())
+ return
+
+
+ def testNeedNewJacobian(self):
"""
- Test useElasticBehavior().
+ Test needNewJacobian().
"""
- material = GenMaxwellIsotropic3D()
- material._createCppHandle()
- material.useElasticBehavior(False)
+ # Default should be False.
+ self.failIf(self.material.needNewJacobian())
+
+ # Changing time step should require new Jacobian.
+ self.material.timeStep(1.0)
+ self.material.timeStep(2.0)
+ self.failUnless(self.material.needNewJacobian())
return
+
-
def test_factory(self):
"""
Test factory method.
Modified: short/3D/PyLith/trunk/unittests/pytests/materials/testmaterials.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/testmaterials.py 2009-06-12 21:04:23 UTC (rev 15220)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/testmaterials.py 2009-06-12 23:35:46 UTC (rev 15221)
@@ -74,8 +74,8 @@
from TestMaxwellIsotropic3D import TestMaxwellIsotropic3D
suite.addTest(unittest.makeSuite(TestMaxwellIsotropic3D))
- #from TestGenMaxwellIsotropic3D import TestGenMaxwellIsotropic3D
- #suite.addTest(unittest.makeSuite(TestGenMaxwellIsotropic3D))
+ from TestGenMaxwellIsotropic3D import TestGenMaxwellIsotropic3D
+ suite.addTest(unittest.makeSuite(TestGenMaxwellIsotropic3D))
from TestPowerLaw3D import TestPowerLaw3D
suite.addTest(unittest.makeSuite(TestPowerLaw3D))
More information about the CIG-COMMITS
mailing list