[cig-commits] r15689 - in short/3D/PyLith/trunk/libsrc: . materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Sep 23 19:35:40 PDT 2009
Author: willic3
Date: 2009-09-23 19:35:39 -0700 (Wed, 23 Sep 2009)
New Revision: 15689
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc
short/3D/PyLith/trunk/libsrc/materials/materialsfwd.hh
Log:
Finished fixing and updating plane strain Maxwell.
Need to create unit tests, python, and bindings.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2009-09-24 02:35:39 UTC (rev 15689)
@@ -82,6 +82,7 @@
materials/ViscoelasticMaxwell.cc \
materials/GenMaxwellIsotropic3D.cc \
materials/MaxwellIsotropic3D.cc \
+ materials/MaxwellPlaneStrain.cc \
materials/PowerLaw3D.cc \
meshio/BinaryIO.cc \
meshio/GMVFile.cc \
Modified: short/3D/PyLith/trunk/libsrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2009-09-24 02:35:39 UTC (rev 15689)
@@ -25,6 +25,8 @@
GenMaxwellIsotropic3D.icc \
MaxwellIsotropic3D.hh \
MaxwellIsotropic3D.icc \
+ MaxwellPlaneStrain.hh \
+ MaxwellPlaneStrain.icc \
PowerLaw3D.hh \
PowerLaw3D.icc \
Metadata.hh \
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2009-09-24 02:35:39 UTC (rev 15689)
@@ -762,7 +762,7 @@
dq * (devStrainTpdt - devStrainT);
} // for
- PetscLogFlops(8 + 7 * tensorSize);
+ PetscLogFlops(9 + 7 * tensorSize);
} // _computeStateVars
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2009-09-24 02:35:39 UTC (rev 15689)
@@ -101,6 +101,8 @@
* @param density Array for density.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars Number of state variables.
+ * @param numStateVars Number of state variables.
*/
void _calcDensity(double* const density,
const double* properties,
@@ -288,7 +290,7 @@
const int initialStrainSize,
const bool computeStateVars);
- /** Compute stress tensor from properties as an viscoelastic material.
+ /** Compute stress tensor from properties as a viscoelastic material.
*
* @param stress Array for stress tensor.
* @param stressSize Size of stress tensor.
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2009-09-24 02:35:39 UTC (rev 15689)
@@ -34,7 +34,7 @@
namespace _MaxwellPlaneStrain{
// Dimension of material
- const int dimension = 3;
+ const int dimension = 2;
/// Number of entries in stress/strain tensors.
const int tensorSize = 3;
@@ -46,7 +46,7 @@
const int numProperties = 4;
/// Physical properties.
- const Material::PropMetaData properties[] = {
+ const Metadata::ParamDescription properties[] = {
{ "density", 1, pylith::topology::FieldBase::SCALAR },
{ "mu", 1, pylith::topology::FieldBase::SCALAR },
{ "lambda", 1, pylith::topology::FieldBase::SCALAR },
@@ -60,18 +60,14 @@
/// Number of state variables.
const int numStateVars = 2;
- /// Size of viscous_strain state variable.
- const int viscousStrainSize = 4;
-
/// State variables.
const Metadata::ParamDescription stateVars[] = {
{ "total_strain", tensorSize, pylith::topology::FieldBase::TENSOR },
- { "viscous_strain", viscousStrainSize,
- pylith::topology::FieldBase::TENSOR },
+ { "viscous_strain", 4, pylith::topology::FieldBase::TENSOR },
};
/// Values expected in state variables spatial database
- const int numDBStateVars = tensorSize + viscousStrainSize;
+ const int numDBStateVars = tensorSize + 4;
const char* dbStateVars[] = { "total-strain-xx",
"total-strain-yy",
"total-strain-xy",
@@ -139,7 +135,7 @@
_updateStateVarsFn(0)
{ // constructor
useElasticBehavior(true);
- _viscousStrain.resize(tensorSize);
+ _viscousStrain.resize(4);
} // constructor
// ----------------------------------------------------------------------
@@ -375,12 +371,19 @@
const int initialStrainSize,
const bool computeStateVars)
{ // _calcStressViscoelastic
+ // Note that there is a problem with the way things are presently set up.
+ // Since the stress tensor is required to have a dimension of 3 for a 2D
+ // problem, the ZZ component of stress is never computed, even though it is
+ // part of the solution. The only way around this at present is to compute
+ // this component as part of postprocessing, which will require knowledge of
+ // the current stress, strain, and viscous strain values as well as the
+ // initial values.
assert(0 != stress);
assert(_MaxwellPlaneStrain::tensorSize == stressSize);
assert(0 != properties);
assert(_numPropsQuadPt == numProperties);
assert(0 != stateVars);
- assert(_numPropsQuadPt == numStateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
assert(0 != initialStress);
@@ -397,16 +400,20 @@
const double mu2 = 2.0 * mu;
const double bulkModulus = lambda + mu2 / 3.0;
- const double e11 = totalStrain[0] - initialStrain[0];
- const double e22 = totalStrain[1] - initialStrain[1];
- const double e12 = totalStrain[2] - initialStrain[2];
- const double e33 = 0.0;
-
- const double traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
- const double meanStressTpdt = bulkModulus * traceStrainTpdt;
+ // Initial stress and strain values
+ const double meanStrainInitial = (initialStrain[0] + initialStrain[1]) / 3.0;
+ const double meanStressInitial = (initialStress[0] + initialStress[1]) / 3.0;
+ const double devStrainInitial[] = {initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ initialStrain[2]};
+ const double devStressInitial[] = {initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ initialStress[2]};
- const double diag[] = { 1.0, 1.0, 1.0, 0.0 };
+ // Mean stress and strain for t + dt
+ const double meanStrainTpdt = (totalStrain[0] + totalStrain[1]) / 3.0;
+ const double meanStressTpdt = 3.0 * bulkModulus *
+ (meanStrainTpdt - meanStrainInitial) + meanStressInitial;
// Get viscous strains
if (computeStateVars) {
@@ -421,16 +428,14 @@
} // else
// Compute new stresses
- double devStressTpdt = 0.0;
+ stress[0] = meanStressTpdt + mu2 * (_viscousStrain[0] - devStrainInitial[0]) +
+ devStressInitial[0];
+ stress[1] = meanStressTpdt + mu2 * (_viscousStrain[1] - devStrainInitial[1]) +
+ devStressInitial[1];
+ stress[2] = mu2 * (_viscousStrain[2] - devStrainInitial[2]) +
+ devStressInitial[2];
- for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStressTpdt = mu2 * _viscousStrain[iComp];
-
- stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
- initialState[iComp];
- } // for
-
- PetscLogFlops(10 + 4 * tensorSize);
+ PetscLogFlops(28);
} // _calcStressViscoelastic
// ----------------------------------------------------------------------
@@ -455,7 +460,7 @@
assert(0 != properties);
assert(_numPropsQuadPt == numProperties);
assert(0 != stateVars);
- assert(_numPropsQuadPt == numStateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
assert(0 != initialStress);
@@ -502,7 +507,7 @@
assert(0 != properties);
assert(_numPropsQuadPt == numProperties);
assert(0 != stateVars);
- assert(_numPropsQuadPt == numStateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_MaxwellPlaneStrain::tensorSize == strainSize);
assert(0 != initialStress);
@@ -593,8 +598,8 @@
const double* totalStrain,
const int strainSize,
const double* initialStress,
- const double* initialStress,
- const int initialStrainSize,
+ const int initialStressSize,
+ const double* initialStrain,
const int initialStrainSize)
{ // _updateStateVarsViscoelastic
assert(0 != stateVars);
@@ -639,7 +644,7 @@
assert(_numVarsQuadPt == numStateVars);
const double maxwellTime = properties[p_maxwellTime];
- const double dtStable = 0.1*maxwellTime;
+ const double dtStable = 0.2*maxwellTime;
return dtStable;
} // _stableTimeStepImplicit
@@ -673,37 +678,35 @@
const int tensorSize = _tensorSize;
const double maxwellTime = properties[p_maxwellTime];
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
- const double e12 = totalStrain[3];
+ const double strainTpdt[] = {totalStrain[0],
+ totalStrain[1],
+ 0.0,
+ totalStrain[2]};
+ const double strainT[] = {stateVars[s_totalStrain+0],
+ stateVars[s_totalStrain+1],
+ 0.0,
+ stateVars[s_totalStrain+2]};
- const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
+ const double meanStrainTpdt = (strainTpdt[0] + strainTpdt[1])/3.0;
+ const double meanStrainT = (strainT[0] + strainT[1])/3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0 };
- const double meanStrainT =
- (properties[_MaxwellPlaneStrain::pidStrainT+0] +
- properties[_MaxwellPlaneStrain::pidStrainT+1] +
- properties[_MaxwellPlaneStrain::pidStrainT+2])/3.0;
-
// Time integration.
- double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
- const double expFac = exp(-_dt/maxwelltime);
+ double 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[_MaxwellPlaneStrain::pidStrainT+iComp] -
- diag[iComp] * meanStrainT;
- _visStrain[iComp] = expFac *
- properties[_MaxwellPlaneStrain::pidVisStrain + iComp] +
+ for (int iComp=0; iComp < 4; ++iComp) {
+ devStrainTpdt = strainTpdt[iComp] - diag[iComp] * meanStrainTpdt;
+ devStrainT = strainT[iComp] - diag[iComp] * meanStrainT;
+ _viscousStrain[iComp] = expFac * stateVars[s_viscousStrain+iComp] +
dq * (devStrainTpdt - devStrainT);
} // for
- PetscLogFlops(8 + 7 * tensorSize);
+ PetscLogFlops(39);
} // _computeStateVars
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh 2009-09-24 02:35:39 UTC (rev 15689)
@@ -25,17 +25,10 @@
#if !defined(pylith_materials_maxwellplanestrain_hh)
#define pylith_materials_maxwellplanestrain_hh
-#include "ElasticMaterial.hh"
+// Include directives ---------------------------------------------------
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
-/// Namespace for pylith package
-namespace pylith {
- namespace materials {
- class MaxwellPlaneStrain;
- class TestMaxwellPlaneStrain; // unit testing
- } // materials
-} // pylith
-
-/// 2-D, isotropic, linear Maxwell viscoelastic material.
+// MaxwellPlaneStrain ---------------------------------------------------
class pylith::materials::MaxwellPlaneStrain : public ElasticMaterial
{ // class MaxwellPlaneStrain
friend class TestMaxwellPlaneStrain; // 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,54 +84,64 @@
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 Number of state variables.
+ * @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 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 +150,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 +223,10 @@
const int,
const double*,
const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
const bool);
/// Member prototype for _calcElasticConsts()
@@ -215,78 +238,86 @@
const double*,
const int,
const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
const int);
- /// Member prototype for _updateProperties()
- typedef void (pylith::materials::MaxwellPlaneStrain::*updateProperties_fn_type)
+ /// Member prototype for _updateStateVars()
+ typedef void (pylith::materials::MaxwellPlaneStrain::*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 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 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.
+ /** Compute stress tensor from properties as a viscoelastic 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 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 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
@@ -296,19 +327,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 +356,105 @@
* @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 :
+ /** 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
- MaxwellPlaneStrain(const MaxwellPlaneStrain& m);
-
- /// Not implemented
- const MaxwellPlaneStrain& operator=(const MaxwellPlaneStrain& 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 +462,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
+ MaxwellPlaneStrain(const MaxwellPlaneStrain&);
+
+ /// Not implemented
+ const MaxwellPlaneStrain& operator=(const MaxwellPlaneStrain&);
+
}; // class MaxwellPlaneStrain
#include "MaxwellPlaneStrain.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc 2009-09-24 02:35:39 UTC (rev 15689)
@@ -27,51 +27,29 @@
_dt = dt;
} // timeStep
-// Set whether elastic or inelastic constitutive relations are used.
-inline
-void
-pylith::materials::MaxwellPlaneStrain::useElasticBehavior(const bool flag) {
- if (flag) {
- _calcStressFn =
- &pylith::materials::MaxwellPlaneStrain::_calcStressElastic;
- _calcElasticConstsFn =
- &pylith::materials::MaxwellPlaneStrain::_calcElasticConstsElastic;
- _updatePropertiesFn =
- &pylith::materials::MaxwellPlaneStrain::_updatePropertiesElastic;
- } else {
- _calcStressFn =
- &pylith::materials::MaxwellPlaneStrain::_calcStressViscoelastic;
- _calcElasticConstsFn =
- &pylith::materials::MaxwellPlaneStrain::_calcElasticConstsViscoelastic;
- _updatePropertiesFn =
- &pylith::materials::MaxwellPlaneStrain::_updatePropertiesViscoelastic;
- } // if/else
-} // useElasticBehavior
-
-// Get flag indicating whether material implements an empty
-inline
-bool
-pylith::materials::MaxwellPlaneStrain::usesUpdateProperties(void) const {
- return true;
-} // usesUpdateProperties
-
// Compute stress tensor from parameters.
inline
void
pylith::materials::MaxwellPlaneStrain::_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 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 +59,44 @@
pylith::materials::MaxwellPlaneStrain::_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::MaxwellPlaneStrain::_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,
+pylith::materials::MaxwellPlaneStrain::_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,
- initialState, initialStateSize);
-} // _updateProperties
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+} // _updateStateVars
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/materialsfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/materialsfwd.hh 2009-09-24 00:30:14 UTC (rev 15688)
+++ short/3D/PyLith/trunk/libsrc/materials/materialsfwd.hh 2009-09-24 02:35:39 UTC (rev 15689)
@@ -36,6 +36,7 @@
class ElasticPlaneStress;
class ElasticIsotropic3D;
class MaxwellIsotropic3D;
+ class MaxwellPlaneStrain;
class GenMaxwellIsotropic3D;
class PowerLaw3D;
More information about the CIG-COMMITS
mailing list