[cig-commits] r7002 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed May 30 12:12:45 PDT 2007
Author: willic3
Date: 2007-05-30 12:12:45 -0700 (Wed, 30 May 2007)
New Revision: 7002
Modified:
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
Log:
More work on Maxwell isotropic material.
Put in calcElasticConsts and updateState.
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-05-30 18:56:06 UTC (rev 7001)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-05-30 19:12:45 UTC (rev 7002)
@@ -42,19 +42,18 @@
const int didViscosity = 3;
/// Parameters
- const int numParameters = 7;
- const int numParamValues[] = { 1, 1, 1, 1, 1, 6, 6};
+ const int numParameters = 6;
+ const int numParamValues[] = { 1, 1, 1, 1, 6, 6};
const char* namesParameters[] =
- {"density", "mu", "lambda" , "viscosity", "maxwellTime", "strainT", "visStrain"};
+ {"density", "mu", "lambda" , "maxwellTime", "strainT", "visStrain"};
/// Indices (order) of parameters
const int pidDensity = 0;
const int pidMu = 1;
const int pidLambda = 2;
- const int pidViscosity = 3;
- const int pidMaxwellTime = 4;
- const int pidStrainT = 5;
- const int pidVisStrain = 6;
+ const int pidMaxwellTime = 3;
+ const int pidStrainT = 4;
+ const int pidVisStrain = 5;
} // _MaxwellIsotropic3D
} // materials
} // pylith
@@ -145,7 +144,6 @@
(*paramVals)[_MaxwellIsotropic3D::pidDensity][0] = density;
(*paramVals)[_MaxwellIsotropic3D::pidMu][0] = mu;
(*paramVals)[_MaxwellIsotropic3D::pidLambda][0] = lambda;
- (*paramVals)[_MaxwellIsotropic3D::pidViscosity][0] = viscosity;
(*paramVals)[_MaxwellIsotropic3D::pidMaxwellTime][0] = maxwelltime;
} // _dbToParameters
@@ -199,7 +197,6 @@
const double density = parameters[_MaxwellIsotropic3D::pidDensity][0];
const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
- const double viscosity = parameters[_MaxwellIsotropic3D::pidViscosity][0];
const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
const double lambda2mu = lambda + 2.0 * mu;
@@ -272,7 +269,7 @@
void
pylith::materials::MaxwellIsotropic3D::_calcElasticConsts(
double_array* const elasticConsts,
- const double_array& parameters,
+ const std::vector<double_array>& parameters,
const double_array& totalStrain)
{ // _calcElasticConsts
assert(0 != elasticConsts);
@@ -283,31 +280,153 @@
const double density = parameters[_MaxwellIsotropic3D::pidDensity];
const double mu = parameters[_MaxwellIsotropic3D::pidMu];
const double lambda = parameters[_MaxwellIsotropic3D::pidLambda];
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
const double lambda2mu = lambda + 2.0 * mu;
+ const double bulkmodulus = lambda + 2.0 * mu/3.0;
+
+ if (useElasticBehavior()) {
+ (*elasticConsts)[ 0] = lambda2mu; // C1111
+ (*elasticConsts)[ 1] = lambda; // C1122
+ (*elasticConsts)[ 2] = lambda; // C1133
+ (*elasticConsts)[ 3] = 0; // C1112
+ (*elasticConsts)[ 4] = 0; // C1123
+ (*elasticConsts)[ 5] = 0; // C1113
+ (*elasticConsts)[ 6] = lambda2mu; // C2222
+ (*elasticConsts)[ 7] = lambda; // C2233
+ (*elasticConsts)[ 8] = 0; // C2212
+ (*elasticConsts)[ 9] = 0; // C2223
+ (*elasticConsts)[10] = 0; // C2213
+ (*elasticConsts)[11] = lambda2mu; // C3333
+ (*elasticConsts)[12] = 0; // C3312
+ (*elasticConsts)[13] = 0; // C3323
+ (*elasticConsts)[14] = 0; // C3313
+ (*elasticConsts)[15] = 2.0 * mu; // C1212
+ (*elasticConsts)[16] = 0; // C1223
+ (*elasticConsts)[17] = 0; // C1213
+ (*elasticConsts)[18] = 2.0 * mu; // C2323
+ (*elasticConsts)[19] = 0; // C2313
+ (*elasticConsts)[20] = 2.0 * mu; // C1313
+ } else {
+ const double timeFrac = 1.0e-5;
+ const int numTerms = 5;
+ double dq = 0.0;
+ if(maxwelltime < timeFrac*_dt) {
+ double fSign = 1.0;
+ double factorial = 1.0;
+ double fraction = 1.0;
+ dq = 1.0;
+ for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
+ factorial *= iTerm;
+ fSign *= -1.0;
+ fraction *= _dt/maxwelltime;
+ dq += fSign*fraction/factorial;
+ } // for
+ } else
+ dq = maxwelltime*(1.0-exp(-_dt/maxwelltime))/_dt;
+ const double visFac = mu*dq/3.0;
+ (*elasticConsts)[ 0] = bulkmodulus + 4.0*dq; // C1111
+ (*elasticConsts)[ 1] = bulkmodulus - 2.0*dq; // C1122
+ (*elasticConsts)[ 2] = (*elasticConsts)[1]; // C1133
+ (*elasticConsts)[ 3] = 0; // C1112
+ (*elasticConsts)[ 4] = 0; // C1123
+ (*elasticConsts)[ 5] = 0; // C1113
+ (*elasticConsts)[ 6] = (*elasticConsts)[0]; // C2222
+ (*elasticConsts)[ 7] = (*elasticConsts)[1]; // C2233
+ (*elasticConsts)[ 8] = 0; // C2212
+ (*elasticConsts)[ 9] = 0; // C2223
+ (*elasticConsts)[10] = 0; // C2213
+ (*elasticConsts)[11] = (*elasticConsts)[0]; // C3333
+ (*elasticConsts)[12] = 0; // C3312
+ (*elasticConsts)[13] = 0; // C3323
+ (*elasticConsts)[14] = 0; // C3313
+ (*elasticConsts)[15] = 3.0 * visFac; // C1212
+ (*elasticConsts)[16] = 0; // C1223
+ (*elasticConsts)[17] = 0; // C1213
+ (*elasticConsts)[18] = (*elasticConsts)[15]; // C2323
+ (*elasticConsts)[19] = 0; // C2313
+ (*elasticConsts)[20] = (*elasticConsts)[15]; // C1313
+ } // else
- (*elasticConsts)[ 0] = lambda2mu; // C1111
- (*elasticConsts)[ 1] = lambda; // C1122
- (*elasticConsts)[ 2] = lambda; // C1133
- (*elasticConsts)[ 3] = 0; // C1112
- (*elasticConsts)[ 4] = 0; // C1123
- (*elasticConsts)[ 5] = 0; // C1113
- (*elasticConsts)[ 6] = lambda2mu; // C2222
- (*elasticConsts)[ 7] = lambda; // C2233
- (*elasticConsts)[ 8] = 0; // C2212
- (*elasticConsts)[ 9] = 0; // C2223
- (*elasticConsts)[10] = 0; // C2213
- (*elasticConsts)[11] = lambda2mu; // C3333
- (*elasticConsts)[12] = 0; // C3312
- (*elasticConsts)[13] = 0; // C3323
- (*elasticConsts)[14] = 0; // C3313
- (*elasticConsts)[15] = 2.0 * mu; // C1212
- (*elasticConsts)[16] = 0; // C1223
- (*elasticConsts)[17] = 0; // C1213
- (*elasticConsts)[18] = 2.0 * mu; // C2323
- (*elasticConsts)[19] = 0; // C2313
- (*elasticConsts)[20] = 2.0 * mu; // C1313
} // _calcElasticConsts
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::MaxwellIsotropic3D::_updateState(
+ const std::vector<double_array>& parameters,
+ const double_array& totalStrain)
+{ // _updateState
+ assert(_MaxwellIsotropic3D::numParameters == parameters.size());
+ assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
+ assert(1 == parameters[_MaxwellIsotropic3D::pidDensity].size());
+ assert(1 == parameters[_MaxwellIsotropic3D::pidMu].size());
+ assert(1 == parameters[_MaxwellIsotropic3D::pidLambda].size());
+ assert(6 == parameters[_MaxwellIsotropic3D::pidStrainT].size());
+ assert(6 == parameters[_MaxwellIsotropic3D::pidVisStrain].size());
+ const double density = parameters[_MaxwellIsotropic3D::pidDensity][0];
+ const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
+ const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
+
+ const double lambda2mu = lambda + 2.0 * mu;
+ const double bulkmodulus = lambda + 2.0 * mu/3.0;
+
+ const double e11 = totalStrain[0];
+ const double e22 = totalStrain[1];
+ const double e33 = totalStrain[2];
+
+ const double traceStrainTpdt = e11 + e22 + e33;
+ const double meanStrainTpdt = traceStrainTpdt/3.0;
+ const double s123 = lambda * traceStrainTpdt;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
+ if (useElasticBehavior()) {
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ parameters[_MaxwellIsotropic3D::pidStrainT][iComp] = totalStrain[iComp];
+ parameters[_MaxwellIsotropic3D::pidVisStrain][iComp] =
+ totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
+ } // for
+ } else {
+ const double meanStrainT = parameters[_MaxwellIsotropic3D::pidStrainT][0] +
+ parameters[_MaxwellIsotropic3D::pidStrainT][1] +
+ parameters[_MaxwellIsotropic3D::pidStrainT][2];
+
+ // The code below should probably be in a separate function since it
+ // is used more than once. I should also probably cover the possibility
+ // that Maxwell time is zero (although this should never happen).
+ const double timeFrac = 1.0e-5;
+ const int numTerms = 5;
+ double dq = 0.0;
+ if(maxwelltime < timeFrac*_dt) {
+ double fSign = 1.0;
+ double factorial = 1.0;
+ double fraction = 1.0;
+ dq = 1.0;
+ for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
+ factorial *= iTerm;
+ fSign *= -1.0;
+ fraction *= _dt/maxwelltime;
+ dq += fSign*fraction/factorial;
+ } // for
+ } else
+ dq = maxwelltime*(1.0-exp(-_dt/maxwelltime))/_dt;
+ const double expFac = exp(-_dt/maxwelltime);
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+ for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
+ devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
+ devStrainT = parameters[_MaxwellIsotropic3D::pidStrainT][iComp] -
+ diag[iComp]*meanStrainT;
+ visStrain = expFac*parameters[_MaxwellIsotropic3D::pidVisStrain][iComp] +
+ dq*(devStrainTpdt - devStrainT);
+ parameters[_MaxwellIsotropic3D::pidVisStrain][iComp] = visStrain;
+ parameters[_MaxwellIsotropic3D::pidStrainT][iComp] = totalStrain[iComp];
+ } // for
+ } //else
+} // _calcStress
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2007-05-30 18:56:06 UTC (rev 7001)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2007-05-30 19:12:45 UTC (rev 7002)
@@ -10,44 +10,44 @@
// ----------------------------------------------------------------------
//
-/** @file libsrc/materials/ElasticIsotropic3D.h
+/** @file libsrc/materials/MaxwellIsotropic3D.h
*
- * @brief C++ ElasticIsotropic3D object
+ * @brief C++ MaxwellIsotropic3D object
*
- * 3-D, isotropic, linear elastic material. The physical properties
- * are specified using density, shear-wave speed, and
- * compressional-wave speed. The physical properties are stored
- * internally using density, lambda, and mu, which are directly
+ * 3-D, isotropic, linear Maxwell viscoelastic material. The
+ * physical properties are specified using density, shear-wave speed,
+ * viscosity, and compressional-wave speed. The physical properties are
+ * stored internally using density, lambda, mu, which are directly
* related to the elasticity constants used in the finite-element
- * integration.
+ * integration. The viscosity is stored using Maxwell Time (viscosity/mu).
*/
-#if !defined(pylith_materials_elasticisotropic3d_hh)
-#define pylith_materials_elasticisotropic3d_hh
+#if !defined(pylith_materials_maxwellisotropic3d_hh)
+#define pylith_materials_maxwellisotropic3d_hh
#include "ElasticMaterial.hh"
/// Namespace for pylith package
namespace pylith {
namespace materials {
- class ElasticIsotropic3D;
- class TestElasticIsotropic3D; // unit testing
+ class MaxwellIsotropic3D;
+ class TestMaxwellIsotropic3D; // unit testing
} // materials
} // pylith
-/// 3-D, isotropic, linear elastic material.
-class pylith::materials::ElasticIsotropic3D : public ElasticMaterial
-{ // class ElasticIsotropic3D
- friend class TestElasticIsotropic3D; // unit testing
+/// 3-D, isotropic, linear Maxwell viscoelastic material.
+class pylith::materials::MaxwellIsotropic3D : public ElasticMaterial
+{ // class MaxwellIsotropic3D
+ friend class TestMaxwellIsotropic3D; // unit testing
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
/// Default constructor
- ElasticIsotropic3D(void);
+ MaxwellIsotropic3D(void);
/// Destructor
- ~ElasticIsotropic3D(void);
+ ~MaxwellIsotropic3D(void);
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -76,7 +76,7 @@
*
* @returns Number of parameters
*/
- int _numParameters(void) const;
+ int _numParamValues(int_array* numValues) const;
/** Compute parameters from values in spatial database.
*
@@ -86,7 +86,7 @@
* @param paramVals Array of parameters
* @param dbValues Array of database values
*/
- void _dbToParameters(double_array* paramVals,
+ void _dbToParameters(std::vector<double_array>* paramVals,
const double_array& dbValues) const;
/** Get number of entries in stress/strain tensors.
@@ -115,7 +115,7 @@
* @param parameters Parameters at location
*/
void _calcDensity(double_array* const density,
- const double_array& parameters);
+ const std::vector<double_array>& parameters);
/** Compute stress tensor from parameters.
*
@@ -124,7 +124,7 @@
* @param totalStrain Total strain at locations.
*/
void _calcStress(double_array* const stress,
- const double_array& parameters,
+ const std::vector<double_array>& parameters,
const double_array& totalStrain);
/** Compute derivatives of elasticity matrix from parameters.
@@ -134,23 +134,23 @@
* @param totalStrain Total strain at locations.
*/
void _calcElasticConsts(double_array* const elasticConsts,
- const double_array& parameters,
+ const std::vector<double_array>& parameters,
const double_array& totalStrain);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
/// Not implemented
- ElasticIsotropic3D(const ElasticIsotropic3D& m);
+ MaxwellIsotropic3D(const MaxwellIsotropic3D& m);
/// Not implemented
- const ElasticIsotropic3D& operator=(const ElasticIsotropic3D& m);
+ const MaxwellIsotropic3D& operator=(const MaxwellIsotropic3D& m);
-}; // class ElasticIsotropic3D
+}; // class MaxwellIsotropic3D
-#include "ElasticIsotropic3D.icc" // inline methods
+#include "MaxwellIsotropic3D.icc" // inline methods
-#endif // pylith_materials_elasticisotropic3d_hh
+#endif // pylith_materials_maxwellisotropic3d_hh
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2007-05-30 18:56:06 UTC (rev 7001)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2007-05-30 19:12:45 UTC (rev 7002)
@@ -10,8 +10,8 @@
// ----------------------------------------------------------------------
//
-#if !defined(pylith_materials_elasticisotropic3d_hh)
-#error "ElasticIsotropic3D.icc can only be included from ElasticIsotropic3D.hh"
+#if !defined(pylith_materials_maxwellisotropic3d_hh)
+#error "MaxwellIsotropic3D.icc can only be included from MaxwellIsotropic3D.hh"
#endif
More information about the cig-commits
mailing list