[cig-commits] r14138 - in short/3D/PyLith/branches/pylith-swig/libsrc: . materials
brad at geodynamics.org
brad at geodynamics.org
Tue Feb 24 15:45:34 PST 2009
Author: brad
Date: 2009-02-24 15:45:33 -0800 (Tue, 24 Feb 2009)
New Revision: 14138
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh
short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.hh
short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.icc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.hh
Log:
Updated Maxwell viscoelastic model.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-02-24 23:45:33 UTC (rev 14138)
@@ -58,6 +58,8 @@
materials/ElasticPlaneStrain.cc \
materials/ElasticPlaneStress.cc \
materials/ElasticIsotropic3D.cc \
+ materials/ViscoelasticMaxwell.cc \
+ materials/MaxwellIsotropic3D.cc \
meshio/BinaryIO.cc \
meshio/GMVFile.cc \
meshio/GMVFileAscii.cc \
@@ -90,9 +92,7 @@
# feassemble/ElasticityExplicit.cc \
# feassemble/ElasticityImplicit.cc \
# feassemble/IntegratorElasticity.cc \
-# materials/MaxwellIsotropic3D.cc \
# materials/GenMaxwellIsotropic3D.cc \
-# materials/ViscoelasticMaxwell.cc \
# meshio/CellFilter.cc \
# meshio/CellFilterAvg.cc \
# meshio/DataWriter.cc \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc 2009-02-24 23:45:33 UTC (rev 14138)
@@ -515,7 +515,7 @@
pylith::materials::ElasticMaterial::_updateStateVars(
double* const stateVars,
const int numStateVars,
- double* const properties,
+ const double* properties,
const int numProperties,
const double* totalStrain,
const int strainSize,
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh 2009-02-24 23:45:33 UTC (rev 14138)
@@ -265,7 +265,7 @@
virtual
void _updateStateVars(double* const stateVars,
const int numStateVars,
- double* const properties,
+ const double* properties,
const int numProperties,
const double* totalStrain,
const int strainSize,
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.cc 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.cc 2009-02-24 23:45:33 UTC (rev 14138)
@@ -15,6 +15,7 @@
#include "MaxwellIsotropic3D.hh" // implementation of object methods
#include "ViscoelasticMaxwell.hh" // USES computeVisStrain
+#include "Metadata.hh" // USES Metadata
#include "pylith/utils/array.hh" // USES double_array
@@ -27,11 +28,17 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
+// :QUESTION: Do we ignore initialStrain and only use initialStress
+// and state variables totalStrain and viscousStrain?
+
// ----------------------------------------------------------------------
namespace pylith {
namespace materials {
namespace _MaxwellIsotropic3D{
+ // Dimension of material.
+ const int dimension = 3;
+
/// Number of entries in stress/strain tensors.
const int tensorSize = 6;
@@ -39,61 +46,105 @@
const int numElasticConsts = 21;
/// 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", 6, OTHER_FIELD },
- { "viscous_strain", 6, OTHER_FIELD },
+ const Metadata::ParamDescription properties[] = {
+ { "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 properties spatial database
+ const int numDBProperties = 4;
+ const char* dbProperties[] = {"density", "vs", "vp" , "viscosity"};
- /// Values expected in spatial database
- const int numDBValues = 4;
- const char* namesDBValues[] = {"density", "vs", "vp" , "viscosity"};
+ /// Number of state variables.
+ const int numStateVars = 6;
+
+ /// State variables.
+ const Metadata::ParamDescription stateVars[] = {
+ { "total_strain", 6, pylith::topology::FieldBase::TENSOR },
+ { "viscous_strain", 6, pylith::topology::FieldBase::TENSOR },
+ };
- /// Indices (order) of database values
- const int didDensity = 0;
- const int didVs = 1;
- const int didVp = 2;
- const int didViscosity = 3;
+ // Values expected in state variables spatial database
+ const int numDBStateVars = 12;
+ const char* dbStateVars[] = {"total-strain-xx",
+ "total-strain-yy",
+ "total-strain-zz",
+ "total-strain-xy",
+ "total-strain-yz",
+ "total-strain-xz",
+ "viscous-strain-xx",
+ "viscous-strain-yy",
+ "viscous-strain-zz",
+ "viscous-strain-xy",
+ "viscous-strain-yz",
+ "viscous-strain-xz",
+ };
- /// 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" };
-
} // _MaxwellIsotropic3D
} // materials
} // pylith
+// Indices of physical properties
+const int pylith::materials::MaxwellIsotropic3D::p_density = 0;
+
+const int pylith::materials::MaxwellIsotropic3D::p_mu =
+ pylith::materials::MaxwellIsotropic3D::p_density + 1;
+
+const int pylith::materials::MaxwellIsotropic3D::p_lambda =
+ pylith::materials::MaxwellIsotropic3D::p_mu + 1;
+
+const int pylith::materials::MaxwellIsotropic3D::p_maxwellTime =
+ pylith::materials::MaxwellIsotropic3D::p_maxwellTime + 1;
+
+// Indices of database values (order must match dbProperties)
+const int pylith::materials::MaxwellIsotropic3D::db_density = 0;
+
+const int pylith::materials::MaxwellIsotropic3D::db_vs =
+ pylith::materials::MaxwellIsotropic3D::db_density + 1;
+
+const int pylith::materials::MaxwellIsotropic3D::db_vp =
+ pylith::materials::MaxwellIsotropic3D::db_vs + 1;
+
+const int pylith::materials::MaxwellIsotropic3D::db_viscosity =
+ pylith::materials::MaxwellIsotropic3D::db_vp + 1;
+
+// Indices of state variables
+const int pylith::materials::MaxwellIsotropic3D::s_totalStrain = 0;
+
+const int pylith::materials::MaxwellIsotropic3D::s_viscousStrain =
+ pylith::materials::MaxwellIsotropic3D::s_totalStrain + 6;
+
+// Indices of database values (order must match dbStateVars)
+const int pylith::materials::MaxwellIsotropic3D::db_totalStrain = 0;
+
+const int pylith::materials::MaxwellIsotropic3D::db_viscousStrain =
+ pylith::materials::MaxwellIsotropic3D::db_totalStrain + 6;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::MaxwellIsotropic3D::MaxwellIsotropic3D(void) :
- ElasticMaterial(_MaxwellIsotropic3D::tensorSize,
+ ElasticMaterial(_MaxwellIsotropic3D::dimension,
+ _MaxwellIsotropic3D::tensorSize,
_MaxwellIsotropic3D::numElasticConsts,
- _MaxwellIsotropic3D::namesDBValues,
- _MaxwellIsotropic3D::namesInitialStateDBValues,
- _MaxwellIsotropic3D::numDBValues,
- _MaxwellIsotropic3D::properties,
- _MaxwellIsotropic3D::numProperties),
- _calcElasticConstsFn(&pylith::materials::MaxwellIsotropic3D::_calcElasticConstsElastic),
- _calcStressFn(&pylith::materials::MaxwellIsotropic3D::_calcStressElastic),
- _updatePropertiesFn(&pylith::materials::MaxwellIsotropic3D::_updatePropertiesElastic)
+ Metadata(_MaxwellIsotropic3D::properties,
+ _MaxwellIsotropic3D::numProperties,
+ _MaxwellIsotropic3D::dbProperties,
+ _MaxwellIsotropic3D::numDBProperties,
+ _MaxwellIsotropic3D::stateVars,
+ _MaxwellIsotropic3D::numStateVars,
+ _MaxwellIsotropic3D::dbStateVars,
+ _MaxwellIsotropic3D::numDBStateVars)),
+ _calcElasticConstsFn(0),
+ _calcStressFn(0),
+ _updateStateVarsFn(0)
{ // constructor
- _dimension = 3;
- _visStrain.resize(_MaxwellIsotropic3D::tensorSize);
+ useElasticBehavior(true);
+ _viscousStrain.resize(_tensorSize);
} // constructor
// ----------------------------------------------------------------------
@@ -103,6 +154,29 @@
} // destructor
// ----------------------------------------------------------------------
+// Set whether elastic or inelastic constitutive relations are used.
+void
+pylith::materials::MaxwellIsotropic3D::useElasticBehavior(const bool flag)
+{ // useElasticBehavior
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsElastic;
+ _updateStateVarsFn =
+ &pylith::materials::MaxwellIsotropic3D::_updateStateVarsElastic;
+
+ } else {
+ _calcStressFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcStressViscoelastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsViscoelastic;
+ _updateStateVarsFn =
+ &pylith::materials::MaxwellIsotropic3D::_updateStateVarsViscoelastic;
+ } // if/else
+} // useElasticBehavior
+
+// ----------------------------------------------------------------------
// Compute properties from values in spatial database.
void
pylith::materials::MaxwellIsotropic3D::_dbToProperties(
@@ -111,12 +185,12 @@
{ // _dbToProperties
assert(0 != propValues);
const int numDBValues = dbValues.size();
- assert(_MaxwellIsotropic3D::numDBValues == numDBValues);
+ assert(_MaxwellIsotropic3D::numDBProperties == numDBValues);
- const double density = dbValues[_MaxwellIsotropic3D::didDensity];
- const double vs = dbValues[_MaxwellIsotropic3D::didVs];
- const double vp = dbValues[_MaxwellIsotropic3D::didVp];
- const double viscosity = dbValues[_MaxwellIsotropic3D::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;
@@ -142,13 +216,13 @@
} // if
assert(mu > 0);
- const double maxwelltime = viscosity / mu;
- assert(maxwelltime > 0.0);
+ const double maxwellTime = viscosity / mu;
+ assert(maxwellTime > 0.0);
- propValues[_MaxwellIsotropic3D::pidDensity] = density;
- propValues[_MaxwellIsotropic3D::pidMu] = mu;
- propValues[_MaxwellIsotropic3D::pidLambda] = lambda;
- propValues[_MaxwellIsotropic3D::pidMaxwellTime] = maxwelltime;
+ propValues[p_density] = density;
+ propValues[p_mu] = mu;
+ propValues[p_lambda] = lambda;
+ propValues[p_maxwellTime] = maxwellTime;
PetscLogFlops(7);
} // _dbToProperties
@@ -161,23 +235,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[_MaxwellIsotropic3D::pidDensity] =
- _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::pidDensity],
- densityScale);
- values[_MaxwellIsotropic3D::pidMu] =
- _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::pidMu],
- pressureScale);
- values[_MaxwellIsotropic3D::pidLambda] =
- _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::pidLambda],
- pressureScale);
- values[_MaxwellIsotropic3D::pidMaxwellTime] =
- _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::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
@@ -190,93 +260,107 @@
{ // _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[_MaxwellIsotropic3D::pidDensity] =
- _normalizer->dimensionalize(values[_MaxwellIsotropic3D::pidDensity],
- densityScale);
- values[_MaxwellIsotropic3D::pidMu] =
- _normalizer->dimensionalize(values[_MaxwellIsotropic3D::pidMu],
- pressureScale);
- values[_MaxwellIsotropic3D::pidLambda] =
- _normalizer->dimensionalize(values[_MaxwellIsotropic3D::pidLambda],
- pressureScale);
- values[_MaxwellIsotropic3D::pidMaxwellTime] =
- _normalizer->dimensionalize(values[_MaxwellIsotropic3D::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::MaxwellIsotropic3D::_nondimInitState(double* const values,
- const int nvalues) const
-{ // _nondimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _MaxwellIsotropic3D::numInitialStateDBValues);
+pylith::materials::MaxwellIsotropic3D::_dbToStateVars(
+ double* const stateValues,
+ const double_array& dbValues) const
+{ // _dbToStateVars
+ assert(0 != stateValues);
+ const int numDBValues = dbValues.size();
+ assert(_MaxwellIsotropic3D::numDBStateVars == numDBValues);
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->nondimensionalize(values, nvalues, pressureScale);
+ const int totalSize = 2 * _tensorSize;
+ 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::MaxwellIsotropic3D::_dimInitState(double* const values,
- const int nvalues) const
-{ // _dimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _MaxwellIsotropic3D::numInitialStateDBValues);
-
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->dimensionalize(values, nvalues, pressureScale);
-
- PetscLogFlops(nvalues);
-} // _dimInitState
-
-// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::MaxwellIsotropic3D::_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[_MaxwellIsotropic3D::pidDensity];
+ density[0] = properties[p_density];
} // _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
+
+// ----------------------------------------------------------------------
// Compute viscous strain for current time step.
// material.
void
pylith::materials::MaxwellIsotropic3D::_computeStateVars(
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ 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(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
- const int tensorSize = _MaxwellIsotropic3D::tensorSize;
- const double maxwelltime = properties[_MaxwellIsotropic3D::pidMaxwellTime];
+ const int tensorSize = _tensorSize;
+ const double maxwellTime = properties[p_maxwellTime];
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
@@ -290,23 +374,21 @@
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
const double meanStrainT =
- (properties[_MaxwellIsotropic3D::pidStrainT+0] +
- properties[_MaxwellIsotropic3D::pidStrainT+1] +
- properties[_MaxwellIsotropic3D::pidStrainT+2])/3.0;
+ ( stateVars[s_totalStrain+0] +
+ stateVars[s_totalStrain+1] +
+ stateVars[s_totalStrain+2] ) / 3.0;
// Time integration.
- double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
- const double expFac = exp(-_dt/maxwelltime);
+ 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 = properties[_MaxwellIsotropic3D::pidStrainT+iComp] -
- diag[iComp] * meanStrainT;
- _visStrain[iComp] = expFac *
- properties[_MaxwellIsotropic3D::pidVisStrain + iComp] +
+ devStrainT = stateVars[s_totalStrain+iComp] - diag[iComp] * meanStrainT;
+ _viscousStrain[iComp] = expFac * stateVars[s_viscousStrain+iComp] +
dq * (devStrainTpdt - devStrainT);
} // for
@@ -318,26 +400,35 @@
// material.
void
pylith::materials::MaxwellIsotropic3D::_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(_MaxwellIsotropic3D::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
- const double mu = properties[_MaxwellIsotropic3D::pidMu];
- const double lambda = properties[_MaxwellIsotropic3D::pidLambda];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0 * mu;
const double e11 = totalStrain[0];
@@ -350,12 +441,12 @@
const double traceStrainTpdt = e11 + e22 + e33;
const double s123 = lambda * traceStrainTpdt;
- 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];
+ 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);
} // _calcStressElastic
@@ -365,29 +456,38 @@
// material.
void
pylith::materials::MaxwellIsotropic3D::_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(_MaxwellIsotropic3D::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
const int tensorSize = _MaxwellIsotropic3D::tensorSize;
- const double mu = properties[_MaxwellIsotropic3D::pidMu];
- const double lambda = properties[_MaxwellIsotropic3D::pidLambda];
- const double maxwelltime = properties[_MaxwellIsotropic3D::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;
@@ -403,27 +503,26 @@
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
// Get viscous strains
- if (computeStateVars) {
- pylith::materials::MaxwellIsotropic3D::_computeStateVars(properties,
- numProperties,
- totalStrain,
- strainSize,
- initialState,
- initialStateSize);
- } else {
- memcpy(&_visStrain[0], &properties[_MaxwellIsotropic3D::pidVisStrain],
- tensorSize * sizeof(double));
- } // else
+ if (computeStateVars)
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+ else
+ memcpy(&_viscousStrain[0], &stateVars[s_viscousStrain],
+ tensorSize*sizeof(double));
// Compute new stresses
double devStressTpdt = 0.0;
for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStressTpdt = mu2 * _visStrain[iComp];
+ devStressTpdt = mu2 * _viscousStrain[iComp];
+ // :QUESTION: ASK CHARLES ABOUT COMMENT ON NEXT LINE
// Later I will want to put in initial stresses.
stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
- initialState[iComp];
+ initialStress[iComp];
} // for
PetscLogFlops(7 + 4 * tensorSize);
@@ -433,25 +532,34 @@
// Compute derivative of elasticity matrix at location from properties.
void
pylith::materials::MaxwellIsotropic3D::_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(_MaxwellIsotropic3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
- const double mu = properties[_MaxwellIsotropic3D::pidMu];
- const double lambda = properties[_MaxwellIsotropic3D::pidLambda];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
@@ -486,33 +594,42 @@
// as an elastic material.
void
pylith::materials::MaxwellIsotropic3D::_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(_MaxwellIsotropic3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
-
- const double mu = properties[_MaxwellIsotropic3D::pidMu];
- const double lambda = properties[_MaxwellIsotropic3D::pidLambda];
- const double maxwelltime = properties[_MaxwellIsotropic3D::pidMaxwellTime];
+ assert(0 != initialStress);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
+ 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;
+ 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
@@ -539,94 +656,94 @@
} // _calcElasticConstsViscoelastic
// ----------------------------------------------------------------------
-// Get stable time step for implicit time integration.
-double
-pylith::materials::MaxwellIsotropic3D::_stableTimeStepImplicit(const double* properties,
- const int numProperties) const
-{ // _stableTimeStepImplicit
- assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
-
- const double maxwellTime =
- properties[_MaxwellIsotropic3D::pidMaxwellTime];
- const double dtStable = 0.1*maxwellTime;
-
- return dtStable;
-} // _stableTimeStepImplicit
-
-// ----------------------------------------------------------------------
// Update state variables.
void
-pylith::materials::MaxwellIsotropic3D::_updatePropertiesElastic(
- double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _updatePropertiesElastic
+pylith::materials::MaxwellIsotropic3D::_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)
+{ // _updateStateVarsElastic
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ 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 double maxwelltime = properties[_MaxwellIsotropic3D::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 traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
+ const double meanStrainTpdt = traceStrainTpdt / 3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
- for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
- properties[_MaxwellIsotropic3D::pidStrainT+iComp] = totalStrain[iComp];
- properties[_MaxwellIsotropic3D::pidVisStrain+iComp] =
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ stateVars[s_totalStrain+iComp] = totalStrain[iComp];
+ stateVars[s_viscousStrain+iComp] =
totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
} // for
- PetscLogFlops(3 + 2 * _MaxwellIsotropic3D::tensorSize);
+ PetscLogFlops(3 + 2 * _tensorSize);
_needNewJacobian = true;
-} // _updatePropertiesElastic
+} // _updateStateVarsElastic
// ----------------------------------------------------------------------
// Update state variables.
void
-pylith::materials::MaxwellIsotropic3D::_updatePropertiesViscoelastic(
- double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _updatePropertiesViscoelastic
+pylith::materials::MaxwellIsotropic3D::_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)
+{ // _updateStateVarsViscoelastic
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStrainSize);
- const int tensorSize = _MaxwellIsotropic3D::tensorSize;
+ const int tensorSize = _tensorSize;
- pylith::materials::MaxwellIsotropic3D::_computeStateVars(properties,
- numProperties,
- totalStrain,
- strainSize,
- initialState,
- initialStateSize);
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
- memcpy(&properties[_MaxwellIsotropic3D::pidVisStrain],
- &_visStrain[0],
- tensorSize * sizeof(double));
- memcpy(&properties[_MaxwellIsotropic3D::pidStrainT],
- &totalStrain[0],
- tensorSize * sizeof(double));
+ memcpy(&stateVars[s_totalStrain], totalStrain, tensorSize*sizeof(double));
+ memcpy(&stateVars[s_viscousStrain], &_viscousStrain[0],
+ tensorSize*sizeof(double));
+
_needNewJacobian = false;
+} // _updateStateVarsViscoelastic
-} // _updatePropertiesViscoelastic
-
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.hh 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.hh 2009-02-24 23:45:33 UTC (rev 14138)
@@ -25,17 +25,10 @@
#if !defined(pylith_materials_maxwellisotropic3d_hh)
#define pylith_materials_maxwellisotropic3d_hh
-#include "ElasticMaterial.hh"
+// Include directives ---------------------------------------------------
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
-/// Namespace for pylith package
-namespace pylith {
- namespace materials {
- class MaxwellIsotropic3D;
- class TestMaxwellIsotropic3D; // unit testing
- } // materials
-} // pylith
-
-/// 3-D, isotropic, linear Maxwell viscoelastic material.
+// MaxwellIsotropic3D ---------------------------------------------------
class pylith::materials::MaxwellIsotropic3D : public ElasticMaterial
{ // class MaxwellIsotropic3D
friend class TestMaxwellIsotropic3D; // unit testing
@@ -61,13 +54,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 :
@@ -98,22 +84,18 @@
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.
@@ -122,30 +104,42 @@
*/
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 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.
+ * @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.
@@ -154,42 +148,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 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 _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 tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @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 :
@@ -204,6 +221,10 @@
const int,
const double*,
const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
const bool);
/// Member prototype for _calcElasticConsts()
@@ -215,78 +236,111 @@
const double*,
const int,
const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
const int);
- /// Member prototype for _updateProperties()
- typedef void (pylith::materials::MaxwellIsotropic3D::*updateProperties_fn_type)
+ /// Member prototype for _updateStateVars()
+ typedef void (pylith::materials::MaxwellIsotropic3D::*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.
+ /** 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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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* properties,
+ void _computeStateVars(const double* stateVars,
+ const int numStateVars,
+ const double* properties,
const int numProperties,
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 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 properties Properties at location.
* @param numProperties Number of properties.
- * @param totalStrain Total strain at locations.
+ * @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 tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ * @param computeStateVars Flag indicating to compute updated state variables.
*/
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.
*
* @param stress Array for stress tensor.
* @param stressSize Size of stress tensor.
- * @param properties Properties at locations.
+ * @param properties Properties at location.
* @param numProperties Number of properties.
- * @param totalStrain Total strain at locations.
+ * @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 tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ * @param computeStateVars Flag indicating to compute updated state variables.
*/
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
@@ -296,19 +350,27 @@
* @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 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 _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.
@@ -317,66 +379,80 @@
* @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 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 _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 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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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 _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.
*
+ * @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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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 _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 :
-
- /// Not implemented
- MaxwellIsotropic3D(const MaxwellIsotropic3D& m);
-
- /// Not implemented
- const MaxwellIsotropic3D& operator=(const MaxwellIsotropic3D& m);
-
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
- /// Viscous strain array.
- double_array _visStrain;
+ double_array _viscousStrain; ///< Array for viscous strain tensor
/// Method to use for _calcElasticConsts().
calcElasticConsts_fn_type _calcElasticConstsFn;
@@ -384,9 +460,35 @@
/// 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_mu;
+ static const int p_lambda;
+ static const int p_maxwellTime;
+ static const int db_density;
+ static const int db_vs;
+ static const int db_vp;
+ static const int db_viscosity;
+
+ static const int s_totalStrain;
+ static const int s_viscousStrain;
+ static const int db_totalStrain;
+ static const int db_viscousStrain;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ MaxwellIsotropic3D(const MaxwellIsotropic3D&);
+
+ /// Not implemented
+ const MaxwellIsotropic3D& operator=(const MaxwellIsotropic3D&);
+
}; // class MaxwellIsotropic3D
#include "MaxwellIsotropic3D.icc" // inline methods
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.icc 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/MaxwellIsotropic3D.icc 2009-02-24 23:45:33 UTC (rev 14138)
@@ -14,7 +14,7 @@
#error "MaxwellIsotropic3D.icc can only be included from MaxwellIsotropic3D.hh"
#endif
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
// Set current time step.
@@ -27,51 +27,30 @@
_dt = dt;
} // timeStep
-// Set whether elastic or inelastic constitutive relations are used.
-inline
-void
-pylith::materials::MaxwellIsotropic3D::useElasticBehavior(const bool flag) {
- if (flag) {
- _calcStressFn =
- &pylith::materials::MaxwellIsotropic3D::_calcStressElastic;
- _calcElasticConstsFn =
- &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsElastic;
- _updatePropertiesFn =
- &pylith::materials::MaxwellIsotropic3D::_updatePropertiesElastic;
- } else {
- _calcStressFn =
- &pylith::materials::MaxwellIsotropic3D::_calcStressViscoelastic;
- _calcElasticConstsFn =
- &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsViscoelastic;
- _updatePropertiesFn =
- &pylith::materials::MaxwellIsotropic3D::_updatePropertiesViscoelastic;
- } // if/else
-} // useElasticBehavior
-
-// Get flag indicating whether material implements an empty
-inline
-bool
-pylith::materials::MaxwellIsotropic3D::usesUpdateProperties(void) const {
- return true;
-} // usesUpdateProperties
-
// Compute stress tensor from parameters.
inline
void
pylith::materials::MaxwellIsotropic3D::_calcStress(double* const stress,
const int stressSize,
- const double* parameters,
- const int numParams,
+ 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 bool computeStateVars) {
+ 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
@@ -81,32 +60,45 @@
pylith::materials::MaxwellIsotropic3D::_calcElasticConsts(
double* const elasticConsts,
const int numElasticConsts,
- const double* parameters,
- const int numParams,
+ 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) {
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::MaxwellIsotropic3D::_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::MaxwellIsotropic3D::_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/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.cc 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.cc 2009-02-24 23:45:33 UTC (rev 14138)
@@ -21,45 +21,52 @@
// ----------------------------------------------------------------------
// Compute viscous strain parameter for a linear Maxwell model.
double
-pylith::materials::ViscoelasticMaxwell::computeVisStrain(const double dt,
- const double maxwelltime)
-{ // check parameters and define cutoff values
- if (maxwelltime <= 0.0)
- throw std::runtime_error("Maxwell time must be > 0.");
+pylith::materials::ViscoelasticMaxwell::viscousStrainParam(const double dt,
+ const double maxwellTime)
+{ // viscousStrainParam
+ // Check parameters
+ if (maxwellTime <= 0.0)
+ throw std::runtime_error("Maxwell time must be greater than 0.");
+
+ // Define cutoff values
const double timeFrac = 1.0e-10;
- const int numTerms = 5;
- // Compute viscous strain parameter.
- // The ratio of dt and maxwelltime should never approach timeFrac for any
- // reasonable computation, but I have put in alternative solutions just in
+ // Compute viscous strain parameter. The ratio of dt and
+ // maxwellTime should never approach timeFrac for any reasonable
+ // computation, but I have put in alternative solutions just in
// case.
+
double dq = 0.0;
- // Use series expansion if dt is very small, since default solution blows
- // up otherwise.
- if(dt < timeFrac*maxwelltime) {
+
+ // Use series expansion if dt is very small, since default solution
+ // blows up otherwise.
+
+ if (dt < timeFrac*maxwellTime) {
double fSign = 1.0;
double factorial = 1.0;
double fraction = 1.0;
dq = 1.0;
+
+ const int numTerms = 5;
for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
factorial *= iTerm;
fSign *= -1.0;
- fraction *= dt/maxwelltime;
- dq += fSign*fraction/factorial;
+ fraction *= dt / maxwellTime;
+ dq += fSign * fraction / factorial;
} // for
PetscLogFlops(8*(numTerms-1));
- // Throw away exponential term if maxwelltime is very small.
- } else if (maxwelltime < timeFrac*dt) {
- dq = maxwelltime/dt;
+ } else if (maxwellTime < timeFrac*dt) {
+ // Throw away exponential term if maxwellTime is very small.
+ dq = maxwellTime / dt;
PetscLogFlops(1);
- // Default solution.
} else{
- dq = maxwelltime*(1.0-exp(-dt/maxwelltime))/dt;
+ // Default solution.
+ dq = maxwellTime*(1.0-exp(-dt/maxwellTime))/dt;
PetscLogFlops(6);
} // else
return dq;
-} // computeVisStrain
+} // viscousStrainParam
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.hh 2009-02-24 22:59:50 UTC (rev 14137)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ViscoelasticMaxwell.hh 2009-02-24 23:45:33 UTC (rev 14138)
@@ -38,10 +38,13 @@
/** Compute viscous strain parameter.
*
+ * @param dt Time step.
+ * @param maxwellTime Maxwell time.
+ *
* @returns Viscous strain parameter.
*/
- static double computeVisStrain(const double dt,
- const double maxwelltime);
+ static double viscousStrainParam(const double dt,
+ const double maxwellTime);
}; // class ViscoelasticMaxwell
More information about the CIG-COMMITS
mailing list