[cig-commits] r12465 - in short/3D/PyLith/trunk: libsrc/materials modulesrc/materials pylith/materials unittests/libtests/materials unittests/libtests/materials/data
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Jul 22 09:56:12 PDT 2008
Author: willic3
Date: 2008-07-22 09:56:12 -0700 (Tue, 22 Jul 2008)
New Revision: 12465
Modified:
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/materials/Material.hh
short/3D/PyLith/trunk/libsrc/materials/Material.icc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src
short/3D/PyLith/trunk/pylith/materials/Material.py
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.hh
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3D.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrain.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStress.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1D.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1D.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
Log:
Updated materials and unit tests to use initial stresses.
Things have been set up to eventually have a more general initial state,
primarily to allow specification of an initial strain.
Code compiles and unit tests pass, but initial stresses have not yet been
tested, and they are not yet tested in the unit tests.
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -57,6 +57,12 @@
const int didVs = 1;
const int didVp = 2;
+ /// 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" };
+
} // _ElasticIsotropic3D
} // materials
} // pylith
@@ -68,6 +74,7 @@
ElasticMaterial(_ElasticIsotropic3D::tensorSize,
_ElasticIsotropic3D::numElasticConsts,
_ElasticIsotropic3D::namesDBValues,
+ _ElasticIsotropic3D::namesInitialStateDBValues,
_ElasticIsotropic3D::numDBValues,
_ElasticIsotropic3D::properties,
_ElasticIsotropic3D::numProperties)
@@ -140,6 +147,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStress
assert(0 != stress);
@@ -148,6 +157,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticIsotropic3D::tensorSize == strainSize);
+ assert(_ElasticIsotropic3D::tensorSize == initialStateSize);
const double density = properties[_ElasticIsotropic3D::pidDensity];
const double mu = properties[_ElasticIsotropic3D::pidMu];
@@ -164,14 +174,14 @@
const double s123 = lambda * (e11 + e22 + e33);
- stress[0] = s123 + mu2*e11;
- stress[1] = s123 + mu2*e22;
- stress[2] = s123 + mu2*e33;
- stress[3] = mu2 * e12;
- stress[4] = mu2 * e23;
- stress[5] = mu2 * e13;
+ 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];
- PetscLogFlops(13);
+ PetscLogFlops(19);
} // _calcStress
// ----------------------------------------------------------------------
@@ -183,7 +193,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _calcElasticConsts
assert(0 != elasticConsts);
assert(_ElasticIsotropic3D::numElasticConsts == numElasticConsts);
@@ -191,6 +203,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticIsotropic3D::tensorSize == strainSize);
+ assert(_ElasticIsotropic3D::tensorSize == initialStateSize);
const double density = properties[_ElasticIsotropic3D::pidDensity];
const double mu = properties[_ElasticIsotropic3D::pidMu];
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -84,6 +84,8 @@
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStress(double* const stress,
@@ -92,6 +94,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -102,13 +106,17 @@
* @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 _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Get stable time step for implicit time integration.
*
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -26,12 +26,12 @@
pylith::materials::ElasticMaterial::ElasticMaterial(const int tensorSize,
const int numElasticConsts,
const char** dbValues,
+ const char** initialStateDBValues,
const int numDBValues,
const PropMetaData* properties,
const int numProperties) :
- Material(dbValues, numDBValues, properties, numProperties),
+ Material(tensorSize, dbValues, initialStateDBValues, numDBValues, properties, numProperties),
_numQuadPts(0),
- _tensorSize(tensorSize),
_numElasticConsts(numElasticConsts)
{ // constructor
} // constructor
@@ -70,12 +70,14 @@
assert(_propertiesCell.size() == numQuadPts*totalPropsQuadPt);
assert(totalStrain.size() == numQuadPts*_tensorSize);
assert(_stress.size() == numQuadPts*_tensorSize);
+ assert(_initialStateCell.size() == numQuadPts*_initialStateSize);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
_calcStress(&_stress[iQuad*_tensorSize], _tensorSize,
&_propertiesCell[iQuad*totalPropsQuadPt], totalPropsQuadPt,
&totalStrain[iQuad*_tensorSize], _tensorSize,
- computeStateVars);
+ &_initialStateCell[iQuad*_initialStateSize],
+ _initialStateSize, computeStateVars);
return _stress;
} // calcStress
@@ -91,13 +93,16 @@
assert(_propertiesCell.size() == numQuadPts*totalPropsQuadPt);
assert(totalStrain.size() == numQuadPts*_tensorSize);
assert(_elasticConsts.size() == numQuadPts*_numElasticConsts);
+ assert(_initialStateCell.size() == numQuadPts*_initialStateSize);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
_calcElasticConsts(&_elasticConsts[iQuad*_numElasticConsts],
_numElasticConsts,
&_propertiesCell[iQuad*totalPropsQuadPt],
totalPropsQuadPt,
- &totalStrain[iQuad*_tensorSize], _tensorSize);
+ &totalStrain[iQuad*_tensorSize], _tensorSize,
+ &_initialStateCell[iQuad*_initialStateSize],
+ _initialStateSize);
return _elasticConsts;
} // calcDerivElastic
@@ -138,6 +143,7 @@
_density.resize(numQuadPts*1);
_stress.resize(numQuadPts*_tensorSize);
_elasticConsts.resize(numQuadPts*_numElasticConsts);
+ _initialStateCell.resize(numQuadPts*_tensorSize);
} // if
_getProperties(cell);
@@ -157,7 +163,9 @@
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
_updateProperties(&_propertiesCell[iQuad*totalPropsQuadPt],
totalPropsQuadPt,
- &totalStrain[iQuad*_tensorSize], _tensorSize);
+ &totalStrain[iQuad*_tensorSize], _tensorSize,
+ &_initialStateCell[iQuad*_initialStateSize],
+ _initialStateSize);
_properties->updatePoint(cell, &_propertiesCell[0]);
} // updateProperties
@@ -173,10 +181,21 @@
const int totalPropsQuadPt = _totalPropsQuadPt;
assert(_propertiesCell.size() == numQuadPts*totalPropsQuadPt);
assert(_properties->getFiberDimension(cell) == numQuadPts*totalPropsQuadPt);
+ assert(_initialStateCell.size() == numQuadPts*_initialStateSize);
const real_section_type::value_type* parameterCell =
_properties->restrictPoint(cell);
memcpy(&_propertiesCell[0], parameterCell,
numQuadPts*totalPropsQuadPt*sizeof(double));
+ if (0 != _initialState) {
+ assert(_initialState->getFiberDimension(cell) == numQuadPts*_initialStateSize);
+ const real_section_type::value_type* initialStateValuesCell =
+ _initialState->restrictPoint(cell);
+ memcpy(&_initialStateCell[0], initialStateValuesCell,
+ numQuadPts*_initialStateSize*sizeof(double));
+ } else {
+ for (int iVal=0; iVal < _initialStateSize*numQuadPts; ++iVal)
+ _initialStateCell[iVal] = 0.0;
+ } // if
} // _getProperties
// ----------------------------------------------------------------------
@@ -185,7 +204,9 @@
pylith::materials::ElasticMaterial::_updateProperties(double* const properties,
const int totalPropsQuadPt,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _updateProperties
} // _updateProperties
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -44,6 +44,7 @@
* @param tensorSize Number of entries in stress tensor.
* @param numElasticConsts Number of elastic constants.
* @param dbValues Array of names of database values for material.
+ * @param InitialStateDBValues Names of initial state database values for material.
* @param numDBValues Number of database values.
* @param properties Array of physical property meta data.
* @param numProperties Number of physical properties for material.
@@ -51,6 +52,7 @@
ElasticMaterial(const int tensorSize,
const int numElasticConsts,
const char** dbValues,
+ const char** initialStateDBValues,
const int numDBValues,
const PropMetaData* properties,
const int numProperties);
@@ -183,6 +185,8 @@
* @param numProperties Number of properties.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
+ * @param initialState Initial state variable values.
+ * @param initialStateSize Size of initial state array.
* @param computeStateVars Flag indicating to compute updated state vars.
*/
virtual
@@ -192,6 +196,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars) = 0;
/** Compute derivatives of elasticity matrix from properties.
@@ -202,6 +208,8 @@
* @param numProperties Number of properties.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
+ * @param initialState Initial state variable values.
+ * @param initialStateSize Size of initial state array.
*/
virtual
void _calcElasticConsts(double* const elasticConsts,
@@ -209,7 +217,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize) = 0;
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) = 0;
/** Get stable time step for implicit time integration.
*
@@ -225,12 +235,16 @@
* @param numProperties Number of properties.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
+ * @param initialState Initial state variable values.
+ * @param initialStateSize Size of initial state array.
*/
virtual
void _updateProperties(double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -254,7 +268,6 @@
private :
int _numQuadPts; ///< Number of quadrature points
- int _tensorSize; ///< Number of entries in stress tensor.
int _numElasticConsts; ///< Number of elastic constants.
/** Properties at quadrature points for current cell.
@@ -264,6 +277,13 @@
*/
double_array _propertiesCell;
+ /** Initial state values at quadrature points for current cell.
+ *
+ * size = numQuadPts*initialStateSize
+ * index = iQuadPt*iValue
+ */
+ double_array _initialStateCell;
+
/** Density value at quadrature points for current cell.
*
* size = numQuadPts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -57,6 +57,11 @@
const int didVs = 1;
const int didVp = 2;
+ /// Initial state values expected in spatial database
+ const int numInitialStateDBValues = tensorSize;
+ const char* namesInitialStateDBValues[] = { "stress_xx", "stress_yy",
+ "stress_xy" };
+
} // _ElasticPlaneStrain
} // materials
} // pylith
@@ -68,6 +73,7 @@
ElasticMaterial(_ElasticPlaneStrain::tensorSize,
_ElasticPlaneStrain::numElasticConsts,
_ElasticPlaneStrain::namesDBValues,
+ _ElasticPlaneStrain::namesInitialStateDBValues,
_ElasticPlaneStrain::numDBValues,
_ElasticPlaneStrain::properties,
_ElasticPlaneStrain::numProperties)
@@ -140,6 +146,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStress
assert(0 != stress);
@@ -148,6 +156,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticPlaneStrain::tensorSize == strainSize);
+ assert(_ElasticPlaneStrain::tensorSize == initialStateSize);
const double density = properties[_ElasticPlaneStrain::pidDensity];
const double mu = properties[_ElasticPlaneStrain::pidMu];
@@ -161,15 +170,15 @@
const double s12 = lambda * (e11 + e22);
- stress[0] = s12 + mu2*e11;
- stress[1] = s12 + mu2*e22;
- stress[2] = mu2 * e12;
+ stress[0] = s12 + mu2*e11 + initialState[0];
+ stress[1] = s12 + mu2*e22 + initialState[1];
+ stress[2] = mu2 * e12 + initialState[2];
- PetscLogFlops(8);
+ PetscLogFlops(11);
} // _calcStress
// ----------------------------------------------------------------------
-// Compute density at location from properties.
+// Compute elastic constants at location from properties.
void
pylith::materials::ElasticPlaneStrain::_calcElasticConsts(
double* const elasticConsts,
@@ -177,7 +186,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // calcElasticConsts
assert(0 != elasticConsts);
assert(_ElasticPlaneStrain::numElasticConsts == numElasticConsts);
@@ -185,6 +196,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticPlaneStrain::tensorSize == strainSize);
+ assert(_ElasticPlaneStrain::tensorSize == initialStateSize);
const double density = properties[_ElasticPlaneStrain::pidDensity];
const double mu = properties[_ElasticPlaneStrain::pidMu];
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -84,6 +84,8 @@
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStress(double* const stress,
@@ -92,6 +94,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -102,13 +106,17 @@
* @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 _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Get stable time step for implicit time integration.
*
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -56,6 +56,11 @@
const int didVs = 1;
const int didVp = 2;
+ /// Initial state values expected in spatial database
+ const int numInitialStateDBValues = tensorSize;
+ const char* namesInitialStateDBValues[] = { "stress_xx", "stress_yy",
+ "stress_xy" };
+
} // _ElasticPlaneStress
} // materials
} // pylith
@@ -66,6 +71,7 @@
ElasticMaterial(_ElasticPlaneStress::tensorSize,
_ElasticPlaneStress::numElasticConsts,
_ElasticPlaneStress::namesDBValues,
+ _ElasticPlaneStress::namesInitialStateDBValues,
_ElasticPlaneStress::numDBValues,
_ElasticPlaneStress::properties,
_ElasticPlaneStress::numProperties)
@@ -136,6 +142,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStress
assert(0 != stress);
@@ -144,6 +152,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticPlaneStress::tensorSize == strainSize);
+ assert(_ElasticPlaneStress::tensorSize == initialStateSize);
const double density = properties[_ElasticPlaneStress::pidDensity];
const double mu = properties[_ElasticPlaneStress::pidMu];
@@ -157,9 +166,9 @@
const double e22 = totalStrain[1];
const double e12 = totalStrain[2];
- stress[0] = (2.0*mu2*lambdamu * e11 + mu2*lambda * e22) / lambda2mu;
- stress[1] = (mu2*lambda * e11 + 2.0*mu2*lambdamu * e22) / lambda2mu;
- stress[2] = mu2 * e12;
+ stress[0] = (2.0*mu2*lambdamu * e11 + mu2*lambda * e22) / lambda2mu + initialState[0];
+ stress[1] = (mu2*lambda * e11 + 2.0*mu2*lambdamu * e22) / lambda2mu + initialState[1];
+ stress[2] = mu2 * e12 + initialState[2];
PetscLogFlops(18);
} // _calcStress
@@ -173,7 +182,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // calcElasticConsts
assert(0 != elasticConsts);
assert(_ElasticPlaneStress::numElasticConsts == numElasticConsts);
@@ -181,6 +192,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticPlaneStress::tensorSize == strainSize);
+ assert(_ElasticPlaneStress::tensorSize == initialStateSize);
const double density = properties[_ElasticPlaneStress::pidDensity];
const double mu = properties[_ElasticPlaneStress::pidMu];
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -84,6 +84,8 @@
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStress(double* const stress,
@@ -92,6 +94,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -102,13 +106,17 @@
* @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 _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Get stable time step for implicit time integration.
*
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -53,6 +53,10 @@
const int didDensity = 0;
const int didVs = 1;
const int didVp = 2;
+
+ /// Initial state values expected in spatial database.
+ const int numInitialStateDBValues = tensorSize;
+ const char* namesInitialStateDBValues[] = { "stress_xx" };
} // _ElasticStrain1D
} // materials
@@ -65,6 +69,7 @@
ElasticMaterial(_ElasticStrain1D::tensorSize,
_ElasticStrain1D::numElasticConsts,
_ElasticStrain1D::namesDBValues,
+ _ElasticStrain1D::namesInitialStateDBValues,
_ElasticStrain1D::numDBValues,
_ElasticStrain1D::properties,
_ElasticStrain1D::numProperties)
@@ -136,6 +141,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStress
assert(0 != stress);
@@ -144,15 +151,16 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticStrain1D::tensorSize == strainSize);
+ assert(_ElasticStrain1D::tensorSize == initialStateSize);
const double density = properties[_ElasticStrain1D::pidDensity];
const double lambda = properties[_ElasticStrain1D::pidLambda];
const double mu = properties[_ElasticStrain1D::pidMu];
const double e11 = totalStrain[0];
- stress[0] = (lambda + 2.0*mu) * e11;
+ stress[0] = (lambda + 2.0*mu) * e11 + initialState[0];
- PetscLogFlops(3);
+ PetscLogFlops(4);
} // _calcStress
// ----------------------------------------------------------------------
@@ -164,7 +172,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _calcElasticConsts
assert(0 != elasticConsts);
assert(_ElasticStrain1D::numElasticConsts == numElasticConsts);
@@ -172,6 +182,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticStrain1D::tensorSize == strainSize);
+ assert(_ElasticStrain1D::tensorSize == initialStateSize);
const double density = properties[_ElasticStrain1D::pidDensity];
const double lambda = properties[_ElasticStrain1D::pidLambda];
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -83,6 +83,8 @@
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStress(double* const stress,
@@ -91,6 +93,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -101,13 +105,17 @@
* @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 _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Get stable time step for implicit time integration.
*
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -53,6 +53,11 @@
const int didDensity = 0;
const int didVs = 1;
const int didVp = 2;
+
+ /// Initial state values expected in spatial database
+ const int numInitialStateDBValues = tensorSize;
+ const char* namesInitialStateDBValues[] = { "stress_xx" };
+
} // _ElasticStress1D
} // materials
} // pylith
@@ -63,6 +68,7 @@
ElasticMaterial(_ElasticStress1D::tensorSize,
_ElasticStress1D::numElasticConsts,
_ElasticStress1D::namesDBValues,
+ _ElasticStress1D::namesInitialStateDBValues,
_ElasticStress1D::numDBValues,
_ElasticStress1D::properties,
_ElasticStress1D::numProperties)
@@ -133,6 +139,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStress
assert(0 != stress);
@@ -141,13 +149,14 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticStress1D::tensorSize == strainSize);
+ assert(_ElasticStress1D::tensorSize == initialStateSize);
const double density = properties[_ElasticStress1D::pidDensity];
const double mu = properties[_ElasticStress1D::pidMu];
const double lambda = properties[_ElasticStress1D::pidLambda];
const double e11 = totalStrain[0];
- stress[0] = mu * (3.0*lambda+2.0*mu) / (lambda + mu) * e11;
+ stress[0] = mu * (3.0*lambda+2.0*mu) / (lambda + mu) * e11 + initialState[0];
PetscLogFlops(7);
} // _calcStress
@@ -161,7 +170,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _calcElasticConsts
assert(0 != elasticConsts);
assert(_ElasticStress1D::numElasticConsts == numElasticConsts);
@@ -169,6 +180,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_ElasticStress1D::tensorSize == strainSize);
+ assert(_ElasticStress1D::tensorSize == initialStateSize);
const double density = properties[_ElasticStress1D::pidDensity];
const double mu = properties[_ElasticStress1D::pidMu];
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -84,6 +84,8 @@
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStress(double* const stress,
@@ -92,6 +94,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -102,13 +106,17 @@
* @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 _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Get stable time step for implicit time integration.
*
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -86,6 +86,12 @@
const int didViscosity2 = 7;
const int didViscosity3 = 8;
+ /// Initial state values expected in spatial database
+ const int numInitialStateDBValues = tensorSize;
+ const char* namesInitialStateDBValues[] = { "stress_xx", "stress_yy",
+ "stress_zz", "stress_xy",
+ "stress_yz", "stress_xz" };
+
} // _GenMaxwellIsotropic3D
} // materials
} // pylith
@@ -96,6 +102,7 @@
ElasticMaterial(_GenMaxwellIsotropic3D::tensorSize,
_GenMaxwellIsotropic3D::numElasticConsts,
_GenMaxwellIsotropic3D::namesDBValues,
+ _GenMaxwellIsotropic3D::namesInitialStateDBValues,
_GenMaxwellIsotropic3D::numDBValues,
_GenMaxwellIsotropic3D::properties,
_GenMaxwellIsotropic3D::numProperties),
@@ -185,12 +192,15 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _computeStateVars
assert(0 != properties);
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
const int tensorSize = _GenMaxwellIsotropic3D::tensorSize;
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
@@ -264,6 +274,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStressElastic
assert(0 != stress);
@@ -272,6 +284,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
const double mu = properties[_GenMaxwellIsotropic3D::pidMuTot];
const double lambda = properties[_GenMaxwellIsotropic3D::pidLambdaTot];
@@ -287,12 +300,12 @@
const double traceStrainTpdt = e11 + e22 + e33;
const double s123 = lambda * traceStrainTpdt;
- stress[0] = s123 + mu2*e11;
- stress[1] = s123 + mu2*e22;
- stress[2] = s123 + mu2*e33;
- stress[3] = mu2 * e12;
- stress[4] = mu2 * e23;
- stress[5] = mu2 * e13;
+ stress[0] = s123 + mu2*e11 + initialState[0];
+ stress[1] = s123 + mu2*e22 + initialState[1];
+ stress[2] = s123 + mu2*e33 + initialState[2];
+ stress[3] = mu2 * e12 + initialState[3];
+ stress[4] = mu2 * e23 + initialState[4];
+ stress[5] = mu2 * e13 + initialState[5];
// std::cout << " _calcStressElastic: " << std::endl;
// std::cout << " totalStrain: " << std::endl;
// for (int iComp=0; iComp < _GenMaxwellIsotropic3D::tensorSize; ++iComp)
@@ -302,7 +315,7 @@
// std::cout << " " << (*stress)[iComp];
// std::cout << std::endl;
- PetscLogFlops(13);
+ PetscLogFlops(19);
} // _calcStressElastic
@@ -317,6 +330,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStressViscoelastic
assert(0 != stress);
@@ -325,6 +340,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
const int tensorSize = _GenMaxwellIsotropic3D::tensorSize;
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
@@ -375,7 +391,9 @@
pylith::materials::GenMaxwellIsotropic3D::_computeStateVars(properties,
numProperties,
totalStrain,
- strainSize);
+ strainSize,
+ initialState,
+ initialStateSize);
} else {
memcpy(&_visStrain[0], &properties[_GenMaxwellIsotropic3D::pidVisStrain],
numMaxwellModels * tensorSize * sizeof(double));
@@ -398,15 +416,15 @@
} // for
devStressTpdt = mu2*devStressTpdt;
- // Later I will want to put in initial stresses.
- stress[iComp] =diag[iComp]*meanStressTpdt+devStressTpdt;
+ stress[iComp] =diag[iComp] * meanStressTpdt + devStressTpdt +
+ initialState[iComp];
// std::cout << devStressTpdt << " " << stress[iComp] << std::endl;
// Temporary to get stresses and strains.
// std::cout << " " << stress[iComp] << " " << totalStrain[iComp] << " " << _visStrain << std:: endl;
} // for
- PetscLogFlops((3 + 2*numMaxwellModels) * tensorSize);
+ PetscLogFlops((7 + 2*numMaxwellModels) * tensorSize);
} // _calcStressViscoelastic
// ----------------------------------------------------------------------
@@ -418,7 +436,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _calcElasticConstsElastic
assert(0 != elasticConsts);
assert(_GenMaxwellIsotropic3D::numElasticConsts == numElasticConsts);
@@ -426,6 +446,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
const double mu = properties[_GenMaxwellIsotropic3D::pidMuTot];
const double lambda = properties[_GenMaxwellIsotropic3D::pidLambdaTot];
@@ -468,7 +489,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _calcElasticConstsViscoelastic
assert(0 != elasticConsts);
assert(_GenMaxwellIsotropic3D::numElasticConsts == numElasticConsts);
@@ -476,6 +499,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
const double mu = properties[_GenMaxwellIsotropic3D::pidMuTot];
const double lambda = properties[_GenMaxwellIsotropic3D::pidLambdaTot];
@@ -566,12 +590,15 @@
double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _updatePropertiesElastic
assert(0 != properties);
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
@@ -610,12 +637,15 @@
double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _updatePropertiesViscoelastic
assert(0 != properties);
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_GenMaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_GenMaxwellIsotropic3D::tensorSize == initialStateSize);
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
const int tensorSize = _GenMaxwellIsotropic3D::tensorSize;
@@ -623,7 +653,9 @@
pylith::materials::GenMaxwellIsotropic3D::_computeStateVars(properties,
numProperties,
totalStrain,
- strainSize);
+ strainSize,
+ initialState,
+ initialStateSize);
memcpy(&properties[_GenMaxwellIsotropic3D::pidVisStrain],
&_visStrain[0],
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -113,6 +113,8 @@
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStress(double* const stress,
@@ -121,6 +123,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -131,13 +135,17 @@
* @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 _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Get stable time step for implicit time integration.
*
@@ -152,11 +160,15 @@
* @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 _updateProperties(double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
// PRIVATE TYPEDEFS ///////////////////////////////////////////////////
private :
@@ -169,6 +181,8 @@
const int,
const double*,
const int,
+ const double*,
+ const int,
const bool);
/// Member prototype for _calcElasticConsts()
@@ -178,6 +192,8 @@
const double*,
const int,
const double*,
+ const int,
+ const double*,
const int);
/// Member prototype for _updateProperties()
@@ -185,6 +201,8 @@
(double* const,
const int,
const double*,
+ const int,
+ const double*,
const int);
// PRIVATE METHODS ////////////////////////////////////////////////////
@@ -196,11 +214,15 @@
* @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 int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Compute stress tensor from properties as an elastic material.
*
@@ -210,6 +232,8 @@
* @param numProperties Number of properties.
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStressElastic(double* const stress,
@@ -218,6 +242,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute stress tensor from properties as an viscoelastic material.
@@ -228,6 +254,8 @@
* @param numProperties Number of properties.
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStressViscoelastic(double* const stress,
@@ -236,6 +264,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties as an
@@ -247,13 +277,17 @@
* @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 _calcElasticConstsElastic(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Compute derivatives of elasticity matrix from properties as a
* viscoelastic material.
@@ -264,13 +298,17 @@
* @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 _calcElasticConstsViscoelastic(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Update state variables after solve as an elastic material.
*
@@ -278,11 +316,15 @@
* @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 _updatePropertiesElastic(double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Update state variables after solve as a viscoelastic material.
*
@@ -290,11 +332,15 @@
* @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 _updatePropertiesViscoelastic(double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.icc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -65,11 +65,14 @@
const int numParams,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars) {
assert(0 != _calcStressFn);
CALL_MEMBER_FN(*this, _calcStressFn)(stress, stressSize,
parameters, numParams,
totalStrain, strainSize,
+ initialState, initialStateSize,
computeStateVars);
} // _calcStress
@@ -82,11 +85,14 @@
const double* parameters,
const int numParams,
const double* totalStrain,
- const int strainSize) {
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) {
assert(0 != _calcElasticConstsFn);
CALL_MEMBER_FN(*this, _calcElasticConstsFn)(elasticConsts, numElasticConsts,
parameters, numParams,
- totalStrain, strainSize);
+ totalStrain, strainSize,
+ initialState, initialStateSize);
} // _calcElasticConsts
// Update state variables after solve.
@@ -95,10 +101,13 @@
pylith::materials::GenMaxwellIsotropic3D::_updateProperties(double* const parameters,
const int numParams,
const double* totalStrain,
- const int strainSize) {
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) {
assert(0 != _updatePropertiesFn);
CALL_MEMBER_FN(*this, _updatePropertiesFn)(parameters, numParams,
- totalStrain, strainSize);
+ totalStrain, strainSize,
+ initialState, initialStateSize);
} // _updateProperties
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -29,15 +29,21 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::materials::Material::Material(const char** dbValues,
+pylith::materials::Material::Material(const int tensorSize,
+ const char** dbValues,
+ const char** initialStateDBValues,
const int numDBValues,
const PropMetaData* properties,
const int numProperties) :
_dt(0.0),
_totalPropsQuadPt(0),
_dimension(0),
+ _tensorSize(tensorSize),
+ _initialStateSize(tensorSize),
+ _initialStateDBValues(initialStateDBValues),
_needNewJacobian(false),
_db(0),
+ _initialStateDB(0),
_id(0),
_label(""),
_propMetaData(properties),
@@ -56,10 +62,11 @@
{ // destructor
// Python db object owns database, so just set pointer to null
_db = 0;
+ _initialStateDB = 0;
} // destructor
// ----------------------------------------------------------------------
-// Get physical property parameters from database.
+// Get physical property parameters and initial state (if used) from database.
void
pylith::materials::Material::initialize(const ALE::Obj<Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs,
@@ -94,6 +101,34 @@
_properties->setFiberDimension(cells, fiberDim);
mesh->allocate(_properties);
+ // Fiber dimension for initial stresses is number of quadrature points times
+ // initial state size.
+ const int initialStateFiberDim = _initialStateSize * numQuadPts;
+
+ // Container for data returned in query of initial state database
+ const int initialStateSize = _initialStateSize;
+ double_array initialStateQueryData(initialStateSize);
+
+ // Container of initial state values at cell's quadrature points
+ double_array initialStateCellData(initialStateFiberDim);
+
+ // If initial state is being used, create a section to hold it.
+ if (0 == _initialStateDB)
+ _initialState = NULL;
+ else {
+ _initialState = new real_section_type(mesh->comm(), mesh->debug());
+ assert(!_initialState.isNull());
+ _initialState->setChart(real_section_type::chart_type(*std::min_element(cells->begin(), cells->end()),
+ *std::max_element(cells->begin(), cells->end())+1));
+
+ _initialState->setFiberDimension(cells, initialStateFiberDim);
+ mesh->allocate(_initialState);
+
+ // Setup database for querying
+ _initialStateDB->open();
+ _initialStateDB->queryVals(_initialStateDBValues, initialStateSize);
+ } // if
+
// Setup database for querying
const int numValues = _numDBValues;
_db->open();
@@ -132,13 +167,38 @@
} // if
_dbToProperties(&cellData[totalPropsQuadPt*iQuadPt], queryData);
+ if (0 != _initialStateDB) {
+ const int err2 = _initialStateDB->query(&initialStateQueryData[0],
+ initialStateSize,
+ &quadPts[index],
+ spaceDim, cs);
+
+ if (err2) {
+ std::ostringstream msg;
+ msg << "Could not find initial state values at \n"
+ << "(";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPts[index+i];
+ msg << ") in material " << _label << "\n"
+ << "using spatial database '" << _initialStateDB->label() << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+ memcpy(&initialStateCellData[initialStateSize * iQuadPt],
+ &initialStateQueryData[0],
+ initialStateSize * sizeof(double));
+ } // if
+
} // for
// Insert cell contribution into fields
_properties->updatePoint(*c_iter, &cellData[0]);
+ if (0 != _initialStateDB)
+ _initialState->updatePoint(*c_iter, &initialStateCellData[0]);
} // for
- // Close database
+ // Close databases
_db->close();
+ if (0 != _initialStateDB)
+ _initialStateDB->close();
} // initialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -69,12 +69,16 @@
/** Default constructor.
*
+ * @param tensorSize Array of names of database values for material.
+ * @param initialStateDBValues Names of initial state database values for material.
* @param dbValues Array of names of database values for material.
* @param numDBValues Number of database values.
* @param properties Array of physical property meta data.
* @param numProperties Number of physical properties for material.
*/
- Material(const char** dbValues,
+ Material(const int tensorSize,
+ const char** dbValues,
+ const char** initialStateDBValues,
const int numDBValues,
const PropMetaData* properties,
const int numProperties);
@@ -95,6 +99,12 @@
*/
void db(spatialdata::spatialdb::SpatialDB* value);
+ /** Set database for initial state variables.
+ *
+ * @param value Pointer to database.
+ */
+ void initialStateDB(spatialdata::spatialdb::SpatialDB* value);
+
/** Set identifier of material.
*
* @param value Material identifier
@@ -203,9 +213,14 @@
/// Section containing physical properties of material.
ALE::Obj<real_section_type> _properties;
+
+ /// Section containing the initial state variables for the material.
+ ALE::Obj<real_section_type> _initialState;
int _totalPropsQuadPt; ///< Total # of property values per quad point.
int _dimension; ///< Spatial dimension associated with material.
+ int _tensorSize; ///< Tensor size for material.
+ int _initialStateSize; ///< Initial state size for material.
bool _needNewJacobian; ///< True if need to reform Jacobian, false otherwise.
// PRIVATE MEMBERS ////////////////////////////////////////////////////
@@ -214,6 +229,9 @@
/// Database of parameters for physical properties of material
spatialdata::spatialdb::SpatialDB* _db;
+ /// Database of initial state values for the material
+ spatialdata::spatialdb::SpatialDB* _initialStateDB;
+
int _id; ///< Material identifier
std::string _label; ///< Label of material
@@ -223,6 +241,8 @@
const char** _dbValues; ///< Names of database values
const int _numDBValues; ///< Number of database values
+ const char** _initialStateDBValues; ///< Names of initial state database values
+
}; // class Material
#include "Material.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.icc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.icc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -28,6 +28,13 @@
_db = value;
}
+// Set database for initial state variables.
+inline
+void
+pylith::materials::Material::initialStateDB(spatialdata::spatialdb::SpatialDB* value) {
+ _initialStateDB = value;
+}
+
// Set identifier of material.
inline
void
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -66,6 +66,12 @@
const int didVp = 2;
const int didViscosity = 3;
+ /// 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
@@ -76,6 +82,7 @@
ElasticMaterial(_MaxwellIsotropic3D::tensorSize,
_MaxwellIsotropic3D::numElasticConsts,
_MaxwellIsotropic3D::namesDBValues,
+ _MaxwellIsotropic3D::namesInitialStateDBValues,
_MaxwellIsotropic3D::numDBValues,
_MaxwellIsotropic3D::properties,
_MaxwellIsotropic3D::numProperties),
@@ -155,12 +162,15 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _computeStateVars
assert(0 != properties);
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
const int tensorSize = _MaxwellIsotropic3D::tensorSize;
const double maxwelltime = properties[_MaxwellIsotropic3D::pidMaxwellTime];
@@ -211,6 +221,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStressElastic
assert(0 != stress);
@@ -219,6 +231,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
const double mu = properties[_MaxwellIsotropic3D::pidMu];
const double lambda = properties[_MaxwellIsotropic3D::pidLambda];
@@ -234,14 +247,14 @@
const double traceStrainTpdt = e11 + e22 + e33;
const double s123 = lambda * traceStrainTpdt;
- stress[0] = s123 + mu2*e11;
- stress[1] = s123 + mu2*e22;
- stress[2] = s123 + mu2*e33;
- stress[3] = mu2 * e12;
- stress[4] = mu2 * e23;
- stress[5] = mu2 * e13;
+ 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];
- PetscLogFlops(13);
+ PetscLogFlops(19);
} // _calcStressElastic
// ----------------------------------------------------------------------
@@ -255,6 +268,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars)
{ // _calcStressViscoelastic
assert(0 != stress);
@@ -263,6 +278,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
const int tensorSize = _MaxwellIsotropic3D::tensorSize;
@@ -288,7 +304,9 @@
pylith::materials::MaxwellIsotropic3D::_computeStateVars(properties,
numProperties,
totalStrain,
- strainSize);
+ strainSize,
+ initialState,
+ initialStateSize);
} else {
memcpy(&_visStrain[0], &properties[_MaxwellIsotropic3D::pidVisStrain],
tensorSize * sizeof(double));
@@ -301,10 +319,11 @@
devStressTpdt = mu2 * _visStrain[iComp];
// Later I will want to put in initial stresses.
- stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt;
+ stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
+ initialState[iComp];
} // for
- PetscLogFlops(7 + 3 * tensorSize);
+ PetscLogFlops(7 + 4 * tensorSize);
} // _calcStressViscoelastic
// ----------------------------------------------------------------------
@@ -316,7 +335,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _calcElasticConstsElastic
assert(0 != elasticConsts);
assert(_MaxwellIsotropic3D::numElasticConsts == numElasticConsts);
@@ -324,6 +345,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
const double mu = properties[_MaxwellIsotropic3D::pidMu];
const double lambda = properties[_MaxwellIsotropic3D::pidLambda];
@@ -366,7 +388,9 @@
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _calcElasticConstsViscoelastic
assert(0 != elasticConsts);
assert(_MaxwellIsotropic3D::numElasticConsts == numElasticConsts);
@@ -374,6 +398,7 @@
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
const double mu = properties[_MaxwellIsotropic3D::pidMu];
const double lambda = properties[_MaxwellIsotropic3D::pidLambda];
@@ -433,7 +458,9 @@
double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _updatePropertiesElastic
assert(0 != properties);
assert(_totalPropsQuadPt == numProperties);
@@ -468,19 +495,24 @@
double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize)
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
{ // _updatePropertiesViscoelastic
assert(0 != properties);
assert(_totalPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_MaxwellIsotropic3D::tensorSize == strainSize);
+ assert(_MaxwellIsotropic3D::tensorSize == initialStateSize);
const int tensorSize = _MaxwellIsotropic3D::tensorSize;
pylith::materials::MaxwellIsotropic3D::_computeStateVars(properties,
numProperties,
totalStrain,
- strainSize);
+ strainSize,
+ initialState,
+ initialStateSize);
memcpy(&properties[_MaxwellIsotropic3D::pidVisStrain],
&_visStrain[0],
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -103,6 +103,8 @@
* @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 _calcStress(double* const stress,
const int stressSize,
@@ -110,6 +112,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -120,13 +124,17 @@
* @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 _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Get stable time step for implicit time integration.
*
@@ -141,11 +149,15 @@
* @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 _updateProperties(double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
// PRIVATE TYPEDEFS ///////////////////////////////////////////////////
private :
@@ -158,6 +170,8 @@
const int,
const double*,
const int,
+ const double*,
+ const int,
const bool);
/// Member prototype for _calcElasticConsts()
@@ -167,6 +181,8 @@
const double*,
const int,
const double*,
+ const int,
+ const double*,
const int);
/// Member prototype for _updateProperties()
@@ -174,6 +190,8 @@
(double* const,
const int,
const double*,
+ const int,
+ const double*,
const int);
// PRIVATE METHODS ////////////////////////////////////////////////////
@@ -185,11 +203,15 @@
* @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 int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Compute stress tensor from properties as an elastic material.
*
@@ -199,6 +221,8 @@
* @param numProperties Number of properties.
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStressElastic(double* const stress,
@@ -207,6 +231,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute stress tensor from properties as an viscoelastic material.
@@ -217,6 +243,8 @@
* @param numProperties Number of properties.
* @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 computeStateVars Flag indicating to compute updated state vars.
*/
void _calcStressViscoelastic(double* const stress,
@@ -225,6 +253,8 @@
const int numProperties,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties as an
@@ -236,13 +266,17 @@
* @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 _calcElasticConstsElastic(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Compute derivatives of elasticity matrix from properties as a
* viscoelastic material.
@@ -253,13 +287,17 @@
* @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 _calcElasticConstsViscoelastic(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Update state variables after solve as an elastic material.
*
@@ -267,11 +305,15 @@
* @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 _updatePropertiesElastic(double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
/** Update state variables after solve as a viscoelastic material.
*
@@ -279,11 +321,15 @@
* @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 _updatePropertiesViscoelastic(double* const properties,
const int numProperties,
const double* totalStrain,
- const int strainSize);
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -64,11 +64,14 @@
const int numParams,
const double* totalStrain,
const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
const bool computeStateVars) {
assert(0 != _calcStressFn);
CALL_MEMBER_FN(*this, _calcStressFn)(stress, stressSize,
parameters, numParams,
totalStrain, strainSize,
+ initialState, initialStateSize,
computeStateVars);
} // _calcStress
@@ -81,11 +84,14 @@
const double* parameters,
const int numParams,
const double* totalStrain,
- const int strainSize) {
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) {
assert(0 != _calcElasticConstsFn);
CALL_MEMBER_FN(*this, _calcElasticConstsFn)(elasticConsts, numElasticConsts,
parameters, numParams,
- totalStrain, strainSize);
+ totalStrain, strainSize,
+ initialState, initialStateSize);
} // _calcElasticConsts
// Update state variables after solve.
@@ -94,10 +100,13 @@
pylith::materials::MaxwellIsotropic3D::_updateProperties(double* const parameters,
const int numParams,
const double* totalStrain,
- const int strainSize) {
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) {
assert(0 != _updatePropertiesFn);
CALL_MEMBER_FN(*this, _updatePropertiesFn)(parameters, numParams,
- totalStrain, strainSize);
+ totalStrain, strainSize,
+ initialState, initialStateSize);
} // _updateProperties
// End of file
Modified: short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src 2008-07-22 16:56:12 UTC (rev 12465)
@@ -156,6 +156,40 @@
Material_db_set(self.thisptr, ptrFromHandle(value))
+ property initialStateDB:
+ def __set__(self, value):
+ """
+ Set database for initial state variables.
+ """
+ # create shim for method 'initialStateDB'
+ #embed{ void Material_initialStateDB_set(void* objVptr, void* isVptr)
+ try {
+ assert(0 != objVptr);
+ spatialdata::spatialdb::SpatialDB* initialStateDB =
+ (spatialdata::spatialdb::SpatialDB*) isVptr;
+ ((pylith::materials::Material*) objVptr)->initialStateDB(initialStateDB);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if not value.name == "spatialdata_spatialdb_SpatialDB":
+ raise TypeError, \
+ "Argument must be extension module type " \
+ "'spatialdata::spatialdb::SpatialDB'."
+ if None == value:
+ Material_initialStateDB_set(self.thisptr, NULL)
+ else:
+ Material_initialStateDB_set(self.thisptr, ptrFromHandle(value))
+
+
+
property id:
def __set__(self, value):
"""
Modified: short/3D/PyLith/trunk/pylith/materials/Material.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/Material.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/pylith/materials/Material.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -51,13 +51,18 @@
## \b Properties
## @li \b id Material identifier (from mesh generator)
## @li \b name Name of material
+ ## @li \b useInitialStress Use initial stresses (true) or not (false).
##
## \b Facilities
## @li \b db Database of material property parameters
## @li \b quadrature Quadrature object for numerical integration
+ ## @li \b initialStressDB Database for initial stresses.
import pyre.inventory
+ useInitialStress = pyre.inventory.bool("use_initial_stress", default=False)
+ useInitialStress.meta['tip'] = "Use initial stresses for material."
+
id = pyre.inventory.int("id", default=0)
id.meta['tip'] = "Material identifier (from mesh generator)."
@@ -68,6 +73,11 @@
db = pyre.inventory.facility("db", family="spatial_database",
factory=SimpleDB)
db.meta['tip'] = "Database of material property parameters."
+
+ initialStressDB = pyre.inventory.facility("initial_stress_db",
+ family="spatial_database",
+ factory=SimpleDB)
+ initialStressDB.meta['tip'] = "Database used for initial stresses."
from pylith.feassemble.quadrature.Quadrature import Quadrature
quadrature = pyre.inventory.facility("quadrature", factory=Quadrature)
@@ -133,6 +143,10 @@
assert(None != self.cppHandle)
self.db.initialize()
self.cppHandle.db = self.db.cppHandle
+ if self.initialStressDB != None:
+ self._info.log("Initializing initial stress database.")
+ self.initialStressDB.initialize()
+ self.cppHandle.initialStressDB = self.initialStressDB.cppHandle
self.cppHandle.initialize(mesh.cppHandle, mesh.coordsys.cppHandle,
self.quadrature.cppHandle)
@@ -158,6 +172,10 @@
self.label = self.inventory.label
self.db = self.inventory.db
self.quadrature = self.inventory.quadrature
+ if self.inventory.useInitialStress:
+ self.initialStressDB = self.inventory.initialStressDB
+ else:
+ self.initialStressDB = None
return
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -59,11 +59,16 @@
} // for
const int tensorSize = 9;
+ const int initialStateSize = 9;
double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
+ double_array initialState(initialStateSize);
+ for (int i=0; i < tensorSize; ++i) {
totalStrain[i] = i;
+ initialState[i] = 0;
+ } // for
- material._updateProperties(¶meters[0], numParams, &totalStrain[0], tensorSize);
+ material._updateProperties(¶meters[0], numParams, &totalStrain[0],
+ tensorSize, &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -143,6 +143,9 @@
const int strainSize = material._tensorSize;
double_array strain(data.strain, numQuadPts*strainSize);
+ material._initialState = NULL;
+ material._initialStateSize = material._tensorSize;
+
material.getPropertiesCell(*c_iter, numQuadPts);
const double_array& stress = material.calcStress(strain);
@@ -214,6 +217,9 @@
const int strainSize = material._tensorSize;
double_array strain(data.strain, numQuadPts*strainSize);
+ material._initialState = NULL;
+ material._initialStateSize = material._tensorSize;
+
material.getPropertiesCell(*c_iter, numQuadPts);
const double_array& elasticConsts = material.calcDerivElastic(strain);
@@ -276,6 +282,9 @@
const int numQuadPts = 2;
const int numParams = data.numParameters;
const int numParamsQuadPt = data.numParamsQuadPt;
+
+ material._initialState = NULL;
+ material._initialStateSize = material._tensorSize;
Mesh::label_sequence::iterator c_iter = cells->begin();
@@ -349,7 +358,9 @@
double_array stress(stressSize);
_matElastic->_calcStress(&stress[0], stress.size(),
data->parameterData, data->numParamsQuadPt,
- data->strain, stressSize, computeStateVars);
+ data->strain, stressSize,
+ data->initialState, data->numInitialStateValues,
+ computeStateVars);
const double* stressE = data->stress;
CPPUNIT_ASSERT(0 != stressE);
@@ -393,7 +404,9 @@
double_array elasticConsts(numConsts);
_matElastic->_calcElasticConsts(&elasticConsts[0], numConsts,
data->parameterData, data->numParamsQuadPt,
- data->strain, strainSize);
+ data->strain, strainSize,
+ data->initialState,
+ data->numInitialStateValues);
const double* elasticConstsE = data->elasticConsts;
CPPUNIT_ASSERT(0 != elasticConstsE);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -59,11 +59,16 @@
} // for
const int tensorSize = 9;
+ const int initialStateSize = 9;
double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
+ double_array initialState(initialStateSize);
+ for (int i=0; i < tensorSize; ++i) {
totalStrain[i] = i;
+ initialState[i] = 0;
+ } // for
- material._updateProperties(¶meters[0], numParams, &totalStrain[0], tensorSize);
+ material._updateProperties(¶meters[0], numParams, &totalStrain[0],
+ tensorSize, &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -59,11 +59,16 @@
} // for
const int tensorSize = 9;
+ const int initialStateSize = 9;
double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
+ double_array initialState(initialStateSize);
+ for (int i=0; i < tensorSize; ++i) {
totalStrain[i] = i;
+ initialState[i] = 0;
+ } // for
- material._updateProperties(¶meters[0], numParams, &totalStrain[0], tensorSize);
+ material._updateProperties(¶meters[0], numParams, &totalStrain[0],
+ tensorSize, &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -59,11 +59,16 @@
} // for
const int tensorSize = 9;
+ const int initialStateSize = 9;
+ double_array initialState(initialStateSize);
double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
+ for (int i=0; i < tensorSize; ++i) {
totalStrain[i] = i;
+ initialState[i] = 0;
+ } // for
- material._updateProperties(¶meters[0], numParams, &totalStrain[0], tensorSize);
+ material._updateProperties(¶meters[0], numParams, &totalStrain[0],
+ tensorSize, &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -59,11 +59,16 @@
} // for
const int tensorSize = 9;
+ const int initialStateSize = 9;
double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
+ double_array initialState(initialStateSize);
+ for (int i=0; i < tensorSize; ++i) {
totalStrain[i] = i;
+ initialState[i] = i;
+ } // for
- material._updateProperties(¶meters[0], numParams, &totalStrain[0], tensorSize);
+ material._updateProperties(¶meters[0], numParams, &totalStrain[0],
+ tensorSize, &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -120,9 +120,13 @@
const int numParamsQuadPt = data.numParamsQuadPt;
const int tensorSize = 6;
+ const int initialStateSize = 6;
double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
+ double_array initialState(initialStateSize);
+ for (int i=0; i < tensorSize; ++i) {
totalStrain[i] = i;
+ initialState[i] = 0;
+ } // for
const double meanStrain =
(totalStrain[0] + totalStrain[1] + totalStrain[2]) / 3.0;
@@ -154,7 +158,8 @@
} // for
material._updateProperties(¶meters[0], numParamsQuadPt,
- &totalStrain[0], tensorSize);
+ &totalStrain[0], tensorSize,
+ &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParamsQuadPt; ++i)
@@ -204,6 +209,7 @@
const int numMaxwellModels = 3;
const int tensorSize = 6;
+ const int initialStateSize = 6;
const double mu = 3.0e10;
material.useElasticBehavior(false);
@@ -218,12 +224,14 @@
double_array totalStrainTpdt(tensorSize);
double_array totalStrainT(tensorSize);
double_array visStrainT(numMaxwellModels * tensorSize);
+ double_array initialState(initialStateSize);
for (int i=0; i < tensorSize; ++i) {
totalStrainTpdt[i] = i;
totalStrainT[i] = totalStrainTpdt[i] / 2.0;
visStrainT[i] = totalStrainTpdt[i] / 4.0;
visStrainT[i + tensorSize] = totalStrainTpdt[i] / 4.0;
visStrainT[i + 2 * tensorSize] = totalStrainTpdt[i] / 4.0;
+ initialState[i] = 0;
} // for
const double meanStrainTpdt =
@@ -281,7 +289,8 @@
} // for
material._updateProperties(¶meters[0], numParamsQuadPt,
- &totalStrainTpdt[0], tensorSize);
+ &totalStrainTpdt[0], tensorSize,
+ &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParamsQuadPt; ++i)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -51,6 +51,22 @@
} // testDB
// ----------------------------------------------------------------------
+// Test initialStateDB()
+void
+pylith::materials::TestMaterial::testInitialStateDB(void)
+{ // testInitialStateDB
+ const char* label = "my_database";
+ spatialdata::spatialdb::SimpleDB initialStateDB;
+ initialStateDB.label(label);
+
+ ElasticIsotropic3D material;
+ material.initialStateDB(&initialStateDB);
+
+ CPPUNIT_ASSERT(0 != material._initialStateDB);
+ CPPUNIT_ASSERT(0 == strcmp(label, material._initialStateDB->label()));
+} // testInitialStateDB
+// ----------------------------------------------------------------------
+
// Test id()
void
pylith::materials::TestMaterial::testID(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -40,6 +40,7 @@
CPPUNIT_TEST_SUITE( TestMaterial );
CPPUNIT_TEST( testDB );
+ CPPUNIT_TEST( testInitialStateDB );
CPPUNIT_TEST( testID );
CPPUNIT_TEST( testLabel );
CPPUNIT_TEST( testTimeStep );
@@ -54,6 +55,9 @@
/// Test db()
void testDB(void);
+ /// Test initialStateDB()
+ void testInitialStateDB(void);
+
/// Test id()
void testID(void);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -120,9 +120,13 @@
const int numParamsQuadPt = data.numParamsQuadPt;
const int tensorSize = 6;
+ const int initialStateSize = 6;
double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
+ double_array initialState(initialStateSize);
+ for (int i=0; i < tensorSize; ++i) {
totalStrain[i] = i;
+ initialState[i] = 0;
+ } // for
const double meanStrain =
(totalStrain[0] + totalStrain[1] + totalStrain[2]) / 3.0;
@@ -149,7 +153,8 @@
} // for
material._updateProperties(¶meters[0], numParamsQuadPt,
- &totalStrain[0], tensorSize);
+ &totalStrain[0], tensorSize,
+ &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParamsQuadPt; ++i)
@@ -205,13 +210,16 @@
const double maxwelltime = viscosity / mu;
const int tensorSize = 6;
+ const int initialStateSize = 6;
double_array totalStrainTpdt(tensorSize);
double_array totalStrainT(tensorSize);
double_array visStrainT(tensorSize);
+ double_array initialState(initialStateSize);
for (int i=0; i < tensorSize; ++i) {
totalStrainTpdt[i] = i;
totalStrainT[i] = totalStrainTpdt[i] / 2.0;
visStrainT[i] = totalStrainTpdt[i] / 4.0;
+ initialState[i] = 0;
} // for
const double meanStrainTpdt =
@@ -252,7 +260,8 @@
} //for
material._updateProperties(¶meters[0], numParamsQuadPt,
- &totalStrainTpdt[0], tensorSize);
+ &totalStrainTpdt[0], tensorSize,
+ &initialState[0], initialStateSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParamsQuadPt; ++i)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3D.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3D.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -38,7 +38,10 @@
self.dimension = 3
self.numDBValues = 3
+ self.numInitialStateValues = 6
self.dbValues = ["density", "vs", "vp"]
+ self.initialStateDBValues = ["stress_xx", "stress_yy", "stress_zz",
+ "stress_xy", "stress_yz", "stress_xz" ]
self.numParameters = 3
self.numParamValues = [1, 1, 1]
self.parameterNames = ["density", "mu", "lambda"]
@@ -47,11 +50,13 @@
vsA = 3000.0
vpA = vsA*3**0.5
strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.6e-4]
+ initialStateA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
+ initialStateB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
self.dbData = numpy.array([ [densityA, vsA, vpA],
[densityB, vsB, vpB] ],
@@ -63,6 +68,11 @@
self.parameterData = numpy.array([ [densityA, muA, lambdaA],
[densityB, muB, lambdaB] ],
dtype=numpy.float64)
+
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.numLocs = 2
numElasticConsts = 21
@@ -71,18 +81,21 @@
self.strain = numpy.array([strainA, strainB],
dtype=numpy.float64)
+
self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, densityA, muA, lambdaA)
+ self._calcStress(strainA, densityA, muA, lambdaA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, densityB, muB, lambdaB)
+ self._calcStress(strainB, densityB, muB, lambdaB,
+ initialStateB)
return
- def _calcStress(self, strainV, densityV, muV, lambdaV):
+ def _calcStress(self, strainV, densityV, muV, lambdaV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -115,6 +128,7 @@
C1313], dtype=numpy.float64)
strain = numpy.reshape(strainV, (6,1))
+ initialState = numpy.reshape(initialStateV, (6,1))
elastic = numpy.array([ [C1111, C1122, C1133, C1112, C1123, C1113],
[C1122, C2222, C2233, C2212, C2223, C2213],
[C1133, C2233, C3333, C3312, C3323, C3313],
@@ -122,7 +136,7 @@
[C1123, C2223, C3323, C1223, C2323, C2313],
[C1113, C2213, C3313, C1213, C2313, C1313] ],
dtype=numpy.float64)
- stress = numpy.dot(elastic, strain)
+ stress = numpy.dot(elastic, strain) + initialState
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::ElasticIsotropic3DData::_numDBValues = 3;
+const int pylith::materials::ElasticIsotropic3DData::_numInitialStateValues = 6;
+
const int pylith::materials::ElasticIsotropic3DData::_numParameters = 3;
const int pylith::materials::ElasticIsotropic3DData::_numParamsQuadPt = 3;
@@ -39,6 +41,15 @@
"vp",
};
+const char* pylith::materials::ElasticIsotropic3DData::_initialStateDBValues[] = {
+"stress_xx",
+"stress_yy",
+"stress_zz",
+"stress_xy",
+"stress_yz",
+"stress_xz",
+};
+
const double pylith::materials::ElasticIsotropic3DData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -48,6 +59,21 @@
2.07846097e+03,
};
+const double pylith::materials::ElasticIsotropic3DData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticIsotropic3DData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -57,6 +83,21 @@
2.88000000e+09,
};
+const double pylith::materials::ElasticIsotropic3DData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticIsotropic3DData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -141,14 +182,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -60,12 +60,16 @@
# Material information
self.dimension = None
self.numDBValues = None
+ self.numInitialStateValues = None
self.numParameters = None
self.numParamsQuadPt = None
self.numParamValues = None
self.dbValues = None
+ self.initialStateDBValues = None
self.dbData = None
+ self.initialStateDBData = None
self.parameterData = None
+ self.initialState = None
# Elastic material information
self.numLocs = None
@@ -106,6 +110,9 @@
self.data.addScalar(vtype="int", name="_numDBValues",
value=self.numDBValues,
format="%d")
+ self.data.addScalar(vtype="int", name="_numInitialStateValues",
+ value=self.numInitialStateValues,
+ format="%d")
self.data.addScalar(vtype="int", name="_numParameters",
value=self.numParameters,
format="%d")
@@ -117,11 +124,20 @@
format="%d", ncols=1)
self.data.addArray(vtype="char*", name="_dbValues", values=self.dbValues,
format="\"%s\"", ncols=1)
+ self.data.addArray(vtype="char*", name="_initialStateDBValues",
+ values=self.initialStateDBValues,
+ format="\"%s\"", ncols=1)
self.data.addArray(vtype="double", name="_dbData", values=self.dbData,
format="%16.8e", ncols=1)
+ self.data.addArray(vtype="double", name="_initialStateDBData",
+ values=self.initialStateDBData,
+ format="%16.8e", ncols=1)
self.data.addArray(vtype="double", name="_parameterData",
values=self.parameterData,
format="%16.8e", ncols=1)
+ self.data.addArray(vtype="double", name="_initialState",
+ values=self.initialState,
+ format="%16.8e", ncols=1)
self.data.addScalar(vtype="int", name="_numLocs", value=self.numLocs,
format="%d")
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrain.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrain.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -38,7 +38,9 @@
self.dimension = 2
self.numDBValues = 3
+ self.numInitialStateValues = 3
self.dbValues = ["density", "vs", "vp"]
+ self.initialStateDBValues = ["stress_xx", "stress_yy", "stress_xy"]
self.numParameters = 3
self.numParamValues = [1, 1, 1]
self.parameterNames = ["density", "mu", "lambda"]
@@ -47,11 +49,13 @@
vsA = 3000.0
vpA = vsA*3**0.5
strainA = [1.1e-4, 2.2e-4, 3.3e-4]
+ initialStateA = [0.0, 0.0, 0.0]
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
strainB = [1.2e-4, 2.3e-4, 3.4e-4]
+ initialStateB = [0.0, 0.0, 0.0]
self.dbData = numpy.array([ [densityA, vsA, vpA],
[densityB, vsB, vpB] ],
@@ -63,6 +67,12 @@
self.parameterData = numpy.array([ [densityA, muA, lambdaA],
[densityB, muB, lambdaB] ],
dtype=numpy.float64)
+
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.numLocs = 2
numElasticConsts = 6
@@ -76,13 +86,15 @@
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, densityA, muA, lambdaA)
+ self._calcStress(strainA, densityA, muA, lambdaA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, densityB, muB, lambdaB)
+ self._calcStress(strainB, densityB, muB, lambdaB,
+ initialStateB)
return
- def _calcStress(self, strainV, densityV, muV, lambdaV):
+ def _calcStress(self, strainV, densityV, muV, lambdaV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -101,7 +113,8 @@
[C1122, C2222, C2212],
[C1112, C2212, C1212] ],
dtype=numpy.float64)
- stress = numpy.dot(elastic, strain)
+ initialState = numpy.reshape(initialStateV, (3,1))
+ stress = numpy.dot(elastic, strain) + initialState
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::ElasticPlaneStrainData::_numDBValues = 3;
+const int pylith::materials::ElasticPlaneStrainData::_numInitialStateValues = 3;
+
const int pylith::materials::ElasticPlaneStrainData::_numParameters = 3;
const int pylith::materials::ElasticPlaneStrainData::_numParamsQuadPt = 3;
@@ -39,6 +41,12 @@
"vp",
};
+const char* pylith::materials::ElasticPlaneStrainData::_initialStateDBValues[] = {
+"stress_xx",
+"stress_yy",
+"stress_xy",
+};
+
const double pylith::materials::ElasticPlaneStrainData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -48,6 +56,15 @@
2.07846097e+03,
};
+const double pylith::materials::ElasticPlaneStrainData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticPlaneStrainData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -57,6 +74,15 @@
2.88000000e+09,
};
+const double pylith::materials::ElasticPlaneStrainData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticPlaneStrainData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -99,14 +125,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStress.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStress.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStress.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -38,7 +38,9 @@
self.dimension = 2
self.numDBValues = 3
+ self.numInitialStateValues = 3
self.dbValues = ["density", "vs", "vp"]
+ self.initialStateDBValues = ["stress_xx", "stress_yy", "stress_xy"]
self.numParameters = 3
self.numParamValues = [1, 1, 1]
self.parameterNames = ["density", "mu", "lambda"]
@@ -47,11 +49,13 @@
vsA = 3000.0
vpA = vsA*3**0.5
strainA = [1.1e-4, 2.2e-4, 3.3e-4]
+ initialStateA = [0.0, 0.0, 0.0]
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
strainB = [1.2e-4, 2.3e-4, 3.4e-4]
+ initialStateB = [0.0, 0.0, 0.0]
self.dbData = numpy.array([ [densityA, vsA, vpA],
[densityB, vsB, vpB] ],
@@ -63,6 +67,11 @@
self.parameterData = numpy.array([ [densityA, muA, lambdaA],
[densityB, muB, lambdaB] ],
dtype=numpy.float64)
+
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.numLocs = 2
numElasticConsts = 6
@@ -76,13 +85,15 @@
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, densityA, muA, lambdaA)
+ self._calcStress(strainA, densityA, muA, lambdaA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, densityB, muB, lambdaB)
+ self._calcStress(strainB, densityB, muB, lambdaB,
+ initialStateB)
return
- def _calcStress(self, strainV, densityV, muV, lambdaV):
+ def _calcStress(self, strainV, densityV, muV, lambdaV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -97,11 +108,12 @@
C1212], dtype=numpy.float64)
strain = numpy.reshape(strainV, (3,1))
+ initialState = numpy.reshape(initialStateV, (3,1))
elastic = numpy.array([ [C1111, C1122, C1112],
[C1122, C2222, C2212],
[C1112, C2212, C1212] ],
dtype=numpy.float64)
- stress = numpy.dot(elastic, strain)
+ stress = numpy.dot(elastic, strain) + initialState
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::ElasticPlaneStressData::_numDBValues = 3;
+const int pylith::materials::ElasticPlaneStressData::_numInitialStateValues = 3;
+
const int pylith::materials::ElasticPlaneStressData::_numParameters = 3;
const int pylith::materials::ElasticPlaneStressData::_numParamsQuadPt = 3;
@@ -39,6 +41,12 @@
"vp",
};
+const char* pylith::materials::ElasticPlaneStressData::_initialStateDBValues[] = {
+"stress_xx",
+"stress_yy",
+"stress_xy",
+};
+
const double pylith::materials::ElasticPlaneStressData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -48,6 +56,15 @@
2.07846097e+03,
};
+const double pylith::materials::ElasticPlaneStressData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticPlaneStressData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -57,6 +74,15 @@
2.88000000e+09,
};
+const double pylith::materials::ElasticPlaneStressData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticPlaneStressData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -99,14 +125,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1D.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1D.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -38,7 +38,9 @@
self.dimension = 1
self.numDBValues = 3
+ self.numInitialStateValues = 1
self.dbValues = ["density", "vs", "vp"]
+ self.initialStateDBValues = ["stress_xx"]
self.numParameters = 3
self.numParamValues = [1, 1, 1]
self.parameterNames = ["density", "mu", "lambda"]
@@ -47,11 +49,13 @@
vsA = 3000.0
vpA = vsA*3**0.5
strainA = [1.1e-4]
+ initialStateA = [0.0]
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
strainB = [1.2e-4]
+ initialStateB = [0.0]
self.dbData = numpy.array([ [densityA, vsA, vpA],
[densityB, vsB, vpB] ],
@@ -65,6 +69,11 @@
self.parameterData = numpy.array([ [densityA, muA, lambdaA],
[densityB, muB, lambdaB] ],
dtype=numpy.float64)
+
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.numLocs = 2
numElasticConsts = 1
@@ -78,13 +87,15 @@
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, densityA, lambda2muA)
+ self._calcStress(strainA, densityA, lambda2muA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, densityB, lambda2muB)
+ self._calcStress(strainB, densityB, lambda2muB,
+ initialStateB)
return
- def _calcStress(self, strainV, densityV, lambda2muV):
+ def _calcStress(self, strainV, densityV, lambda2muV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -92,9 +103,10 @@
elasticConsts = numpy.array([C1111], dtype=numpy.float64)
strain = numpy.reshape(strainV, (1,1))
+ initialState = numpy.reshape(initialStateV, (1,1))
elastic = numpy.array([ [C1111] ],
dtype=numpy.float64)
- stress = numpy.dot(elastic, strain)
+ stress = numpy.dot(elastic, strain) + initialState
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::ElasticStrain1DData::_numDBValues = 3;
+const int pylith::materials::ElasticStrain1DData::_numInitialStateValues = 1;
+
const int pylith::materials::ElasticStrain1DData::_numParameters = 3;
const int pylith::materials::ElasticStrain1DData::_numParamsQuadPt = 3;
@@ -39,6 +41,10 @@
"vp",
};
+const char* pylith::materials::ElasticStrain1DData::_initialStateDBValues[] = {
+"stress_xx",
+};
+
const double pylith::materials::ElasticStrain1DData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -48,6 +54,11 @@
2.07846097e+03,
};
+const double pylith::materials::ElasticStrain1DData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticStrain1DData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -57,6 +68,11 @@
2.88000000e+09,
};
+const double pylith::materials::ElasticStrain1DData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticStrain1DData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -81,14 +97,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1D.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1D.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -38,7 +38,9 @@
self.dimension = 1
self.numDBValues = 3
+ self.numInitialStateValues = 1
self.dbValues = ["density", "vs", "vp"]
+ self.initialStateDBValues = ["stress_xx"]
self.numParameters = 3
self.numParamValues = [1, 1, 1]
self.parameterNames = ["density", "mu", "lambda"]
@@ -47,11 +49,13 @@
vsA = 3000.0
vpA = vsA*3**0.5
strainA = [1.1e-4]
+ initialStateA = [0.0]
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
strainB = [1.2e-4]
+ initialStateB = [0.0]
self.dbData = numpy.array([ [densityA, vsA, vpA],
[densityB, vsB, vpB] ],
@@ -63,6 +67,11 @@
self.parameterData = numpy.array([ [densityA, muA, lambdaA],
[densityB, muB, lambdaB] ],
dtype=numpy.float64)
+
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.numLocs = 2
numElasticConsts = 1
@@ -76,13 +85,15 @@
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, densityA, muA, lambdaA)
+ self._calcStress(strainA, densityA, muA, lambdaA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, densityB, muB, lambdaB)
+ self._calcStress(strainB, densityB, muB, lambdaB,
+ initialStateB)
return
- def _calcStress(self, strainV, densityV, muV, lambdaV):
+ def _calcStress(self, strainV, densityV, muV, lambdaV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -90,9 +101,10 @@
elasticConsts = numpy.array([C1111], dtype=numpy.float64)
strain = numpy.reshape(strainV, (1,1))
+ initialState = numpy.reshape(initialStateV, (1,1))
elastic = numpy.array([ [C1111] ],
dtype=numpy.float64)
- stress = numpy.dot(elastic, strain)
+ stress = numpy.dot(elastic, strain) + initialState
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::ElasticStress1DData::_numDBValues = 3;
+const int pylith::materials::ElasticStress1DData::_numInitialStateValues = 1;
+
const int pylith::materials::ElasticStress1DData::_numParameters = 3;
const int pylith::materials::ElasticStress1DData::_numParamsQuadPt = 3;
@@ -39,6 +41,10 @@
"vp",
};
+const char* pylith::materials::ElasticStress1DData::_initialStateDBValues[] = {
+"stress_xx",
+};
+
const double pylith::materials::ElasticStress1DData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -48,6 +54,11 @@
2.07846097e+03,
};
+const double pylith::materials::ElasticStress1DData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticStress1DData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -57,6 +68,11 @@
2.88000000e+09,
};
+const double pylith::materials::ElasticStress1DData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::ElasticStress1DData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -81,14 +97,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElastic.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElastic.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -40,9 +40,12 @@
self.numMaxwellModels = 3
self.numDBValues = 3 + 2*self.numMaxwellModels
+ self.numInitialStateValues = self.tensorSize
self.dbValues = ["density", "vs", "vp",
"shear_ratio_1", "shear_ratio_2", "shear_ratio_3",
"viscosity_1", "viscosity_2", "viscosity_3"]
+ self.initialStateDBValues = ["stress_xx", "stress_yy", "stress_zz",
+ "stress_xy", "stress_yz", "stress_xy"]
self.numParameters = 7
self.numParamValues = [1, 1, 1,
self.numMaxwellModels, self.numMaxwellModels,
@@ -68,6 +71,7 @@
lambdaA = vpA*vpA*densityA - 2.0*muA
elasDataA = [densityA, vsA, vpA]
matDbA = elasDataA + shearRatioA + viscosityA
+ initialStateA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
densityB = 2000.0
vsB = 1200.0
@@ -83,6 +87,7 @@
lambdaB = vpB*vpB*densityB - 2.0*muB
elasDataB = [densityB, vsB, vpB]
matDbB = elasDataB + shearRatioB + viscosityB
+ initialStateB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
matDb = matDbA + matDbB
@@ -109,18 +114,24 @@
self.density = numpy.array([densityA, densityB], dtype=numpy.float64)
self.strain = numpy.array([strainA, strainB], dtype=numpy.float64)
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, muA, lambdaA)
+ self._calcStress(strainA, muA, lambdaA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, muB, lambdaB)
+ self._calcStress(strainB, muB, lambdaB,
+ initialStateB)
return
- def _calcStress(self, strainV, muV, lambdaV):
+ def _calcStress(self, strainV, muV, lambdaV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -153,6 +164,7 @@
C1313], dtype=numpy.float64)
strain = numpy.reshape(strainV, (6,1))
+ initialState = numpy.reshape(initialStateV, (6,1))
elastic = numpy.array([ [C1111, C1122, C1133, C1112, C1123, C1113],
[C1122, C2222, C2233, C2212, C2223, C2213],
[C1133, C2233, C3333, C3312, C3323, C3313],
@@ -160,7 +172,7 @@
[C1123, C2223, C3323, C1223, C2323, C2313],
[C1113, C2213, C3313, C1213, C2313, C1313] ],
dtype=numpy.float64)
- stress = numpy.dot(elastic, strain)
+ stress = numpy.dot(elastic, strain) + initialState
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::GenMaxwellIsotropic3DElasticData::_numDBValues = 9;
+const int pylith::materials::GenMaxwellIsotropic3DElasticData::_numInitialStateValues = 6;
+
const int pylith::materials::GenMaxwellIsotropic3DElasticData::_numParameters = 7;
const int pylith::materials::GenMaxwellIsotropic3DElasticData::_numParamsQuadPt = 33;
@@ -49,6 +51,15 @@
"viscosity_3",
};
+const char* pylith::materials::GenMaxwellIsotropic3DElasticData::_initialStateDBValues[] = {
+"stress_xx",
+"stress_yy",
+"stress_zz",
+"stress_xy",
+"stress_yz",
+"stress_xy",
+};
+
const double pylith::materials::GenMaxwellIsotropic3DElasticData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -70,6 +81,21 @@
1.00000000e+20,
};
+const double pylith::materials::GenMaxwellIsotropic3DElasticData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::GenMaxwellIsotropic3DElasticData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -139,6 +165,21 @@
0.00000000e+00,
};
+const double pylith::materials::GenMaxwellIsotropic3DElasticData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::GenMaxwellIsotropic3DElasticData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -223,14 +264,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DElasticData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -44,9 +44,12 @@
self.numMaxwellModels = 3
self.numDBValues = 3+2*self.numMaxwellModels
+ self.numInitialStateValues = self.tensorSize
self.dbValues = ["density", "vs", "vp" ,
"shear_ratio_1", "shear_ratio_2", "shear_ratio_3",
"viscosity_1", "viscosity_2", "viscosity_3"]
+ self.initialStateDBValues = ["stress_xx", "stress_yy", "stress_zz",
+ "stress_xy", "stress_yz", "stress_xy"]
self.numParameters = 7
self.numParamValues = [1, 1, 1,
self.numMaxwellModels, self.numMaxwellModels,
@@ -72,6 +75,7 @@
matDbA = elasDataA + shearRatioA + viscosityA
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
+ initialStateA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
densityB = 2000.0
vsB = 1200.0
@@ -87,6 +91,7 @@
matDbB = elasDataB + shearRatioB + viscosityB
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
+ initialStateB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
matDb = matDbA + matDbB
self.dbData = numpy.array(matDb, dtype=numpy.float64)
@@ -125,6 +130,10 @@
self.density = numpy.array([densityA, densityB], dtype=numpy.float64)
self.strain = numpy.array([strainA, strainB], dtype=numpy.float64)
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
dtype=numpy.float64)
@@ -132,18 +141,21 @@
(self.elasticConsts[0,:], self.stress[0,:]) = \
self._calcStress(strainA, muA, lambdaA,
shearRatioA, maxwellTimeA,
- strainTA, visStrainA)
+ strainTA, visStrainA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
self._calcStress(strainB, muB, lambdaB,
shearRatioB, maxwellTimeB,
- strainTB, visStrainB)
+ strainTB, visStrainB,
+ initialStateB)
self.dtStableImplicit = 0.1*min(min(maxwellTimeA),
min(maxwellTimeB))
return
def _calcStress(self, strainV, muV, lambdaV,
- shearRatioV, maxwellTimeV, strainTV, visStrainV):
+ shearRatioV, maxwellTimeV, strainTV, visStrainV,
+ initialStateV):
"""
Compute stress and derivative of elasticity matrix.
This assumes behavior is always viscoelastic.
@@ -236,7 +248,8 @@
dq[model] * deltaStrain
devStressTpdt += shearRatioV[model] * visStrain
devStressTpdt = elasFac * devStressTpdt
- stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt
+ stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt + \
+ initialStateV[iComp]
stress = numpy.reshape(stressV, (6,1))
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numDBValues = 9;
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numInitialStateValues = 6;
+
const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numParameters = 7;
const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numParamsQuadPt = 33;
@@ -49,6 +51,15 @@
"viscosity_3",
};
+const char* pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStateDBValues[] = {
+"stress_xx",
+"stress_yy",
+"stress_zz",
+"stress_xy",
+"stress_yz",
+"stress_xy",
+};
+
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -70,6 +81,21 @@
1.00000000e+20,
};
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -139,6 +165,21 @@
6.70000000e-04,
};
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -223,14 +264,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -39,13 +39,17 @@
int numDBValues; ///< Number of database values
int numParameters; ///< Number of parameters
+ int numInitialStateValues; ///< Number of initial state values.
int numParamsQuadPt; ///< Number of parameters per quadrature point.
int* numParamValues; ///< Number of values for each parameter
char** dbValues; ///< Aray of names of database values;
+ char** initialStateDBValues; ///< Names of initial state database values;
double* dbData; ///< Array of database values at locations
+ double* initialStateDBData; ///< Initial state database values at locations
double* parameterData; ///< Array of parameter values at locations
+ double* initialState; ///< Initial state values at locations
};
#endif // pylith_materials_materialdata_hh
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -38,7 +38,10 @@
self.dimension = 3
self.numDBValues = 4
+ self.numInitialStateValues = 6
self.dbValues = ["density", "vs", "vp", "viscosity"]
+ self.initialStateDBValues = ["stress_xx", "stress_yy", "stress_zz",
+ "stress_xy", "stress_yz", "stress_xy"]
self.numParameters = 6
self.numParamValues = [1, 1, 1, 1, 6, 6]
self.parameterNames = ["density", "mu", "lambda", "maxwellTime", "strainT", "visStrain"]
@@ -50,12 +53,14 @@
vpA = vsA*3**0.5
viscosityA = 1.0e18
strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.6e-4]
+ initialStateA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
viscosityB = 1.0e18
strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
+ initialStateB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
self.dbData = numpy.array([ [densityA, vsA, vpA, viscosityA],
[densityB, vsB, vpB, viscosityB] ],
@@ -87,18 +92,24 @@
self.strain = numpy.array([strainA, strainB],
dtype=numpy.float64)
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, muA, lambdaA)
+ self._calcStress(strainA, muA, lambdaA,
+ initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, muB, lambdaB)
+ self._calcStress(strainB, muB, lambdaB,
+ initialStateB)
return
- def _calcStress(self, strainV, muV, lambdaV):
+ def _calcStress(self, strainV, muV, lambdaV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
"""
@@ -131,6 +142,7 @@
C1313], dtype=numpy.float64)
strain = numpy.reshape(strainV, (6,1))
+ initialState = numpy.reshape(initialStateV, (6,1))
elastic = numpy.array([ [C1111, C1122, C1133, C1112, C1123, C1113],
[C1122, C2222, C2233, C2212, C2223, C2213],
[C1133, C2233, C3333, C3312, C3323, C3313],
@@ -138,7 +150,7 @@
[C1123, C2223, C3323, C1223, C2323, C2313],
[C1113, C2213, C3313, C1213, C2313, C1313] ],
dtype=numpy.float64)
- stress = numpy.dot(elastic, strain)
+ stress = numpy.dot(elastic, strain) + initialState
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::MaxwellIsotropic3DElasticData::_numDBValues = 4;
+const int pylith::materials::MaxwellIsotropic3DElasticData::_numInitialStateValues = 6;
+
const int pylith::materials::MaxwellIsotropic3DElasticData::_numParameters = 6;
const int pylith::materials::MaxwellIsotropic3DElasticData::_numParamsQuadPt = 16;
@@ -43,6 +45,15 @@
"viscosity",
};
+const char* pylith::materials::MaxwellIsotropic3DElasticData::_initialStateDBValues[] = {
+"stress_xx",
+"stress_yy",
+"stress_zz",
+"stress_xy",
+"stress_yz",
+"stress_xy",
+};
+
const double pylith::materials::MaxwellIsotropic3DElasticData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -54,6 +65,21 @@
1.00000000e+18,
};
+const double pylith::materials::MaxwellIsotropic3DElasticData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::MaxwellIsotropic3DElasticData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -89,6 +115,21 @@
0.00000000e+00,
};
+const double pylith::materials::MaxwellIsotropic3DElasticData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::MaxwellIsotropic3DElasticData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -173,14 +214,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py 2008-07-22 16:56:12 UTC (rev 12465)
@@ -38,7 +38,10 @@
self.dimension = 3
self.numDBValues = 4
+ self.numInitialStateValues = 6
self.dbValues = ["density", "vs", "vp" , "viscosity"]
+ self.initialStateDBValues = ["stress_xx", "stress_yy", "stress_zz",
+ "stress_xy", "stress_yz", "stress_xy"]
self.numParameters = 6
self.numParamValues = [1, 1, 1, 1, 6, 6]
self.parameterNames = ["density", "mu", "lambda", "maxwellTime", "strainT", "visStrain"]
@@ -51,6 +54,7 @@
viscosityA = 1.0e18
strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.6e-4]
meanStrainA = (strainA[1] + strainA[2] + strainA[3])/3.0
+ initialStateA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
densityB = 2000.0
vsB = 1200.0
@@ -58,6 +62,7 @@
viscosityB = 1.0e19
strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
meanStrainB = (strainB[1] + strainB[2] + strainB[3])/3.0
+ initialStateB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
@@ -98,22 +103,29 @@
self.strain = numpy.array([strainA, strainB],
dtype=numpy.float64)
+ self.initialStateDBData = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
+ self.initialState = numpy.array([initialStateA, initialStateB],
+ dtype=numpy.float64)
self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
self._calcStress(strainA, muA, lambdaA,
- maxwellTimeA, strainTA, visStrainA)
+ maxwellTimeA, strainTA,
+ visStrainA, initialStateA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
self._calcStress(strainB, muB, lambdaB,
- maxwellTimeB, strainTB, visStrainB)
+ maxwellTimeB, strainTB,
+ visStrainB, initialStateB)
self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
return
- def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, strainTV, visStrainV):
+ def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, strainTV,
+ visStrainV, initialStateV):
"""
Compute stress and derivative of elasticity matrix.
This assumes behavior is always viscoelastic.
@@ -187,7 +199,8 @@
devStrainT = strainTV[iComp] - diag[iComp]*meanStrainT
visStrain = expFac*visStrainV[iComp] + dq*(devStrainTpdt - devStrainT)
devStressTpdt = elasFac*visStrain
- stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt
+ stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt + \
+ initialStateV[iComp]
stress = numpy.reshape(stressV, (6,1))
return (elasticConsts, numpy.ravel(stress))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2008-07-22 16:56:12 UTC (rev 12465)
@@ -19,6 +19,8 @@
const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numDBValues = 4;
+const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numInitialStateValues = 6;
+
const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numParameters = 6;
const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numParamsQuadPt = 16;
@@ -43,6 +45,15 @@
"viscosity",
};
+const char* pylith::materials::MaxwellIsotropic3DTimeDepData::_initialStateDBValues[] = {
+"stress_xx",
+"stress_yy",
+"stress_zz",
+"stress_xy",
+"stress_yz",
+"stress_xy",
+};
+
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -54,6 +65,21 @@
1.00000000e+19,
};
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_initialStateDBData[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_parameterData[] = {
2.50000000e+03,
2.25000000e+10,
@@ -89,6 +115,21 @@
6.70000000e-04,
};
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_initialState[] = {
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_density[] = {
2.50000000e+03,
2.00000000e+03,
@@ -173,14 +214,18 @@
{ // constructor
dimension = _dimension;
numDBValues = _numDBValues;
+ numInitialStateValues = _numInitialStateValues;
numParameters = _numParameters;
numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
dtStableImplicit = _dtStableImplicit;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
+ initialStateDBValues = const_cast<char**>(_initialStateDBValues);
dbData = const_cast<double*>(_dbData);
+ initialStateDBData = const_cast<double*>(_initialStateDBData);
parameterData = const_cast<double*>(_parameterData);
+ initialState = const_cast<double*>(_initialState);
density = const_cast<double*>(_density);
strain = const_cast<double*>(_strain);
stress = const_cast<double*>(_stress);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2008-07-22 13:49:06 UTC (rev 12464)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2008-07-22 16:56:12 UTC (rev 12465)
@@ -41,6 +41,8 @@
static const int _numDBValues;
+ static const int _numInitialStateValues;
+
static const int _numParameters;
static const int _numParamsQuadPt;
@@ -53,10 +55,16 @@
static const char* _dbValues[];
+ static const char* _initialStateDBValues[];
+
static const double _dbData[];
+ static const double _initialStateDBData[];
+
static const double _parameterData[];
+ static const double _initialState[];
+
static const double _density[];
static const double _strain[];
More information about the cig-commits
mailing list