[cig-commits] r14110 - in short/3D/PyLith/branches/pylith-swig/libsrc: . bc materials
brad at geodynamics.org
brad at geodynamics.org
Fri Feb 20 15:29:29 PST 2009
Author: brad
Date: 2009-02-20 15:29:29 -0800 (Fri, 20 Feb 2009)
New Revision: 14110
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.hh
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.icc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.hh
short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.icc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.hh
Log:
Updated ElasticMaterial and ElasticIsotropic3D.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-02-20 23:29:29 UTC (rev 14110)
@@ -52,6 +52,8 @@
feassemble/Quadrature3D.cc \
materials/Metadata.cc \
materials/Material.cc \
+ materials/ElasticMaterial.cc \
+ materials/ElasticIsotropic3D.cc \
meshio/BinaryIO.cc \
meshio/GMVFile.cc \
meshio/GMVFileAscii.cc \
@@ -94,7 +96,6 @@
# feassemble/Quadrature3D.cc \
# materials/ElasticStress1D.cc \
# materials/ElasticStrain1D.cc \
-# materials/ElasticIsotropic3D.cc \
# materials/ElasticMaterial.cc \
# materials/ElasticPlaneStrain.cc \
# materials/ElasticPlaneStress.cc \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc 2009-02-20 23:29:29 UTC (rev 14110)
@@ -205,8 +205,7 @@
&quadPtsGlobal[iSpace], spaceDim, cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find traction values at \n"
- << "(";
+ msg << "Could not find traction values at (";
for (int i=0; i < spaceDim; ++i)
msg << " " << quadPtsGlobal[i+iSpace];
msg << ") for traction boundary condition " << _label << "\n"
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.cc 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.cc 2009-02-20 23:29:29 UTC (rev 14110)
@@ -14,6 +14,8 @@
#include "ElasticIsotropic3D.hh" // implementation of object methods
+#include "Metadata.hh" // USES Metadata
+
#include "pylith/utils/array.hh" // USES double_array
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
@@ -30,58 +32,64 @@
namespace materials {
namespace _ElasticIsotropic3D {
- /// Number of entries in stress tensor.
+ // Dimension of material.
+ const int dimension = 3;
+
+ // Number of entries in stress tensor.
const int tensorSize = 6;
- /// Number of elastic constants (for general 3-D elastic material)
+ // Number of elastic constants (for general 3-D elastic material)
const int numElasticConsts = 21;
- /// Number of physical properties.
+ // Number of physical properties.
const int numProperties = 3;
- /// Physical properties.
- const Material::PropMetaData properties[] = {
- { "density", 1, OTHER_FIELD },
- { "mu", 1, OTHER_FIELD },
- { "lambda", 1, OTHER_FIELD },
+ // Physical properties.
+ const Metadata::ParamDescription properties[] = {
+ { "density", 1, pylith::topology::FieldBase::SCALAR },
+ { "mu", 1, pylith::topology::FieldBase::SCALAR },
+ { "lambda", 1, pylith::topology::FieldBase::SCALAR },
};
- /// Indices of physical properties
- const int pidDensity = 0;
- const int pidMu = pidDensity + 1;
- const int pidLambda = pidMu + 1;
- /// Values expected in spatial database
- const int numDBValues = 3;
- const char* namesDBValues[] = { "density", "vs", "vp" };
+ // Values expected in spatial database
+ const int numDBProperties = 3;
+ const char* dbProperties[] = { "density", "vs", "vp" };
- /// Indices of database values
- 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", "stress_yy",
- "stress_zz", "stress_xy",
- "stress_yz", "stress_xz" };
-
} // _ElasticIsotropic3D
} // materials
} // pylith
+// Indices of physical properties
+const int pylith::materials::ElasticIsotropic3D::p_density = 0;
+const int pylith::materials::ElasticIsotropic3D::p_mu =
+ pylith::materials::ElasticIsotropic3D::p_density + 1;
+
+const int pylith::materials::ElasticIsotropic3D::p_lambda =
+ pylith::materials::ElasticIsotropic3D::p_mu + 1;
+
+// Indices of database values (order must match dbProperties)
+const int pylith::materials::ElasticIsotropic3D::db_density = 0;
+
+const int pylith::materials::ElasticIsotropic3D::db_vs =
+ pylith::materials::ElasticIsotropic3D::db_density + 1;
+
+const int pylith::materials::ElasticIsotropic3D::db_vp =
+ pylith::materials::ElasticIsotropic3D::db_vs + 1;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::ElasticIsotropic3D::ElasticIsotropic3D(void) :
- ElasticMaterial(_ElasticIsotropic3D::tensorSize,
+ ElasticMaterial(_ElasticIsotropic3D::dimension,
+ _ElasticIsotropic3D::tensorSize,
_ElasticIsotropic3D::numElasticConsts,
- _ElasticIsotropic3D::namesDBValues,
- _ElasticIsotropic3D::namesInitialStateDBValues,
- _ElasticIsotropic3D::numDBValues,
- _ElasticIsotropic3D::properties,
- _ElasticIsotropic3D::numProperties)
+ Metadata(_ElasticIsotropic3D::properties,
+ _ElasticIsotropic3D::numProperties,
+ _ElasticIsotropic3D::dbProperties,
+ _ElasticIsotropic3D::numDBProperties,
+ 0, 0,
+ 0, 0))
{ // constructor
- _dimension = 3;
} // constructor
// ----------------------------------------------------------------------
@@ -94,16 +102,16 @@
// Compute properties from values in spatial database.
void
pylith::materials::ElasticIsotropic3D::_dbToProperties(
- double* const propValues,
- const double_array& dbValues) const
+ double* const propValues,
+ const double_array& dbValues) const
{ // _dbToProperties
assert(0 != propValues);
const int numDBValues = dbValues.size();
- assert(_ElasticIsotropic3D::numDBValues == numDBValues);
+ assert(_ElasticIsotropic3D::numDBProperties == numDBValues);
- const double density = dbValues[_ElasticIsotropic3D::didDensity];
- const double vs = dbValues[_ElasticIsotropic3D::didVs];
- const double vp = dbValues[_ElasticIsotropic3D::didVp];
+ const double density = dbValues[db_density];
+ const double vs = dbValues[db_vs];
+ const double vp = dbValues[db_vp];
if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
std::ostringstream msg;
@@ -127,9 +135,9 @@
throw std::runtime_error(msg.str());
} // if
- propValues[_ElasticIsotropic3D::pidDensity] = density;
- propValues[_ElasticIsotropic3D::pidMu] = mu;
- propValues[_ElasticIsotropic3D::pidLambda] = lambda;
+ propValues[p_density] = density;
+ propValues[p_mu] = mu;
+ propValues[p_lambda] = lambda;
PetscLogFlops(6);
} // _dbToProperties
@@ -146,16 +154,14 @@
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
- values[_ElasticIsotropic3D::pidDensity] =
- _normalizer->nondimensionalize(values[_ElasticIsotropic3D::pidDensity],
- densityScale);
- values[_ElasticIsotropic3D::pidMu] =
- _normalizer->nondimensionalize(values[_ElasticIsotropic3D::pidMu],
- pressureScale);
- values[_ElasticIsotropic3D::pidLambda] =
- _normalizer->nondimensionalize(values[_ElasticIsotropic3D::pidLambda],
- pressureScale);
+ values[p_density] =
+ _normalizer->nondimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->nondimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->nondimensionalize(values[p_lambda], pressureScale);
+
PetscLogFlops(3);
} // _nondimProperties
@@ -171,137 +177,120 @@
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
- values[_ElasticIsotropic3D::pidDensity] =
- _normalizer->dimensionalize(values[_ElasticIsotropic3D::pidDensity],
- densityScale);
- values[_ElasticIsotropic3D::pidMu] =
- _normalizer->dimensionalize(values[_ElasticIsotropic3D::pidMu],
- pressureScale);
- values[_ElasticIsotropic3D::pidLambda] =
- _normalizer->dimensionalize(values[_ElasticIsotropic3D::pidLambda],
- pressureScale);
+ values[p_density] =
+ _normalizer->dimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->dimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->dimensionalize(values[p_lambda], pressureScale);
+
PetscLogFlops(3);
} // _dimProperties
// ----------------------------------------------------------------------
-// Nondimensionalize initial state.
-void
-pylith::materials::ElasticIsotropic3D::_nondimInitState(double* const values,
- const int nvalues) const
-{ // _nondimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _ElasticIsotropic3D::numInitialStateDBValues);
-
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->nondimensionalize(values, nvalues, pressureScale);
-
- PetscLogFlops(nvalues);
-} // _nondimInitState
-
-// ----------------------------------------------------------------------
-// Dimensionalize initial state.
-void
-pylith::materials::ElasticIsotropic3D::_dimInitState(double* const values,
- const int nvalues) const
-{ // _dimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _ElasticIsotropic3D::numInitialStateDBValues);
-
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->dimensionalize(values, nvalues, pressureScale);
-
- PetscLogFlops(nvalues);
-} // _dimInitState
-
-// ----------------------------------------------------------------------
// Compute density at location from properties.
void
-pylith::materials::ElasticIsotropic3D::_calcDensity(
- double* const density,
- const double* properties,
- const int numProperties)
+pylith::materials::ElasticIsotropic3D::_calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars)
{ // _calcDensity
assert(0 != density);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 == numStateVars);
- density[0] = properties[_ElasticIsotropic3D::pidDensity];
+ density[0] = properties[p_density];
} // _calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from properties.
void
-pylith::materials::ElasticIsotropic3D::_calcStress(
- double* const stress,
- const int stressSize,
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize,
- const bool computeStateVars)
+pylith::materials::ElasticIsotropic3D::_calcStress(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars)
{ // _calcStress
assert(0 != stress);
assert(_ElasticIsotropic3D::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 == numStateVars);
assert(0 != totalStrain);
assert(_ElasticIsotropic3D::tensorSize == strainSize);
- assert(_ElasticIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_ElasticIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_ElasticIsotropic3D::tensorSize == initialStrainSize);
- const double density = properties[_ElasticIsotropic3D::pidDensity];
- const double mu = properties[_ElasticIsotropic3D::pidMu];
- const double lambda = properties[_ElasticIsotropic3D::pidLambda];
+ const double density = properties[p_density];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0*mu;
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
- const double e12 = totalStrain[3];
- const double e23 = totalStrain[4];
- const double e13 = totalStrain[5];
+ const double e11 = totalStrain[0] + initialStrain[0];
+ const double e22 = totalStrain[1] + initialStrain[1];
+ const double e33 = totalStrain[2] + initialStrain[2];
+ const double e12 = totalStrain[3] + initialStrain[3];
+ const double e23 = totalStrain[4] + initialStrain[4];
+ const double e13 = totalStrain[5] + initialStrain[5];
const double s123 = lambda * (e11 + e22 + e33);
- stress[0] = s123 + mu2*e11 + initialState[0];
- stress[1] = s123 + mu2*e22 + initialState[1];
- stress[2] = s123 + mu2*e33 + initialState[2];
- stress[3] = mu2 * e12 + initialState[3];
- stress[4] = mu2 * e23 + initialState[4];
- stress[5] = mu2 * e13 + initialState[5];
+ stress[0] = s123 + mu2*e11 + initialStress[0];
+ stress[1] = s123 + mu2*e22 + initialStress[1];
+ stress[2] = s123 + mu2*e33 + initialStress[2];
+ stress[3] = mu2 * e12 + initialStress[3];
+ stress[4] = mu2 * e23 + initialStress[4];
+ stress[5] = mu2 * e13 + initialStress[5];
- PetscLogFlops(19);
+ PetscLogFlops(25);
} // _calcStress
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from properties.
void
pylith::materials::ElasticIsotropic3D::_calcElasticConsts(
- double* const elasticConsts,
- const int numElasticConsts,
- const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConsts
assert(0 != elasticConsts);
assert(_ElasticIsotropic3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 == numStateVars);
assert(0 != totalStrain);
assert(_ElasticIsotropic3D::tensorSize == strainSize);
- assert(_ElasticIsotropic3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_ElasticIsotropic3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_ElasticIsotropic3D::tensorSize == initialStrainSize);
- const double density = properties[_ElasticIsotropic3D::pidDensity];
- const double mu = properties[_ElasticIsotropic3D::pidMu];
- const double lambda = properties[_ElasticIsotropic3D::pidLambda];
+ const double density = properties[p_density];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
@@ -334,8 +323,11 @@
// ----------------------------------------------------------------------
// Get stable time step for implicit time integration.
double
-pylith::materials::ElasticIsotropic3D::_stableTimeStepImplicit(const double* properties,
- const int numProperties) const
+pylith::materials::ElasticIsotropic3D::_stableTimeStepImplicit(
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const
{ // _stableTimeStepImplicit
return pylith::PYLITH_MAXDOUBLE;
} // _stableTimeStepImplicit
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.hh 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticIsotropic3D.hh 2009-02-20 23:29:29 UTC (rev 14110)
@@ -25,17 +25,10 @@
#if !defined(pylith_materials_elasticisotropic3d_hh)
#define pylith_materials_elasticisotropic3d_hh
-#include "ElasticMaterial.hh"
+// Include directives ---------------------------------------------------
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
-/// Namespace for pylith package
-namespace pylith {
- namespace materials {
- class ElasticIsotropic3D;
- class TestElasticIsotropic3D; // unit testing
- } // materials
-} // pylith
-
-/// 3-D, isotropic, linear elastic material.
+// ElasticIsotropic3D ---------------------------------------------------
class pylith::materials::ElasticIsotropic3D : public ElasticMaterial
{ // class ElasticIsotropic3D
friend class TestElasticIsotropic3D; // unit testing
@@ -79,22 +72,6 @@
void _dimProperties(double* const values,
const int nvalues) const;
- /** Nondimensionalize initial state.
- *
- * @param values Array of initial state values.
- * @param nvalues Number of values.
- */
- void _nondimInitState(double* const values,
- const int nvalues) const;
-
- /** Dimensionalize initial state.
- *
- * @param values Array of initial state values.
- * @param nvalues Number of values.
- */
- void _dimInitState(double* const values,
- const int nvalues) const;
-
/** Compute density from properties.
*
* @param density Array for density.
@@ -103,31 +80,42 @@
*/
void _calcDensity(double* const density,
const double* properties,
- const int numProperties);
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars);
- /** Compute stress tensor from properties. If the state variables
- * are from the previous time step, then the computeStateVars flag
- * should be set to true so that the state variables are updated
- * (but not stored) when computing the stresses.
+ /** Compute stress tensor from properties and state variables. If
+ * the state variables are from the previous time step, then the
+ * computeStateVars flag should be set to true so that the state
+ * variables are updated (but not stored) when computing the
+ * stresses.
*
* @param stress Array for stress tensor.
* @param stressSize Size of stress tensor.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
- * @param computeStateVars Flag indicating to compute updated state vars.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ * @param computeStateVars Flag indicating to compute updated state variables.
*/
void _calcStress(double* const stress,
const int stressSize,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -136,27 +124,52 @@
* @param numElasticConsts Number of elastic constants.
* @param properties Properties at location.
* @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
* @param totalStrain Total strain at location.
* @param strainSize Size of strain tensor.
- * @param initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
*/
void _calcElasticConsts(double* const elasticConsts,
const int numElasticConsts,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize);
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
/** Get stable time step for implicit time integration.
*
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
* @returns Time step
*/
double _stableTimeStepImplicit(const double* properties,
- const int numProperties) const;
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const;
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ static const int p_density;
+ static const int p_mu;
+ static const int p_lambda;
+ static const int db_density;
+ static const int db_vs;
+ static const int db_vp;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc 2009-02-20 23:29:29 UTC (rev 14110)
@@ -29,6 +29,10 @@
#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::ElasticMaterial::ElasticMaterial(const int dimension,
const int tensorSize,
@@ -61,167 +65,57 @@
// Initialize material by getting physical property parameters from
// database.
void
-initialize(const topology::Mesh& mesh,
- feassemble::Quadrature<topology::Mesh>* quadrature)
+pylith::materials::ElasticMaterial::initialize(
+ const topology::Mesh& mesh,
+ feassemble::Quadrature<topology::Mesh>* quadrature)
{ // initialize
Material::initialize(mesh, quadrature);
- delete _initialStress; _initialStress = 0;
- delete _initialStrain; _initialStrain = 0;
+ assert(0 != quadrature);
+ _numQuadPts = quadrature->numQuadPts();
- if (0 != _dbInitialStress)
- _initializeInitialStress(mesh, quadrature);
-
- if (0 != _dbInitialStrain)
- _initializeInitialStrain(mesh, quadrature);
+ _initializeInitialStress(mesh, quadrature);
+ _initializeInitialStrain(mesh, quadrature);
+ _allocateCellArrays();
} // initialize
- // Get quadrature information
- const int numQuadPts = quadrature->numQuadPts();
- const int spaceDim = quadrature->spaceDim();
+// ----------------------------------------------------------------------
+// Retrieve parameters for physical properties and state variables for cell.
+void
+pylith::materials::ElasticMaterial::retrievePropsAndVars(const int cell)
+{ // retrievePropsAndVars
+ assert(0 != _properties);
+ assert(0 != _stateVars);
- // Create arrays for querying
- const int tensorSize = tensorSize();
- double_array quadPtsGlobal(numQuadPts*spaceDim);
- double_array valuesQuery(tensorSize);
- double_array valuesCell(numQuadPts*tensorSize);
+ const ALE::Obj<RealSection>& propertiesSection = _properties->section();
+ assert(!propertiesSection.isNull());
+ propertiesSection->restrictPoint(cell, &_propertiesCell[0],
+ _propertiesCell.size());
+
+ const ALE::Obj<RealSection>& stateVarsSection = _stateVars->section();
+ assert(!stateVarsSection.isNull());
+ stateVarsSection->restrictPoint(cell, &_stateVarsCell[0],
+ _stateVarsCell.size());
+} // retrievePropsAndVars
- // Get cells associated with material
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
-
- // Create field to hold initial stress state.
- delete _initialStress; _initialStress = 0;
- if (0 != _dbInitialStress) {
- _initialStress = new topology::Field<topology::Mesh>(mesh);
- assert(0 != _initialStress);
- fiberDim = numQuadPts * tensorSize;
- _initialStress->newSection(cells, fiberDim);
- _initialStress->allocate();
- _initialStress->zero();
- } // if
- const ALE::Obj<RealSection>& initialStressSection =
- (0 != _initialStress) ? _initialStress->section() : 0;
-
- // Create field to hold initial strain state.
- delete _initialStrain; _initialStrain = 0;
- if (0 != _dbInitialStrain) {
- _initialStrain = new topology::Field<topology::Mesh>(mesh);
- assert(0 != _initialStrain);
- fiberDim = numQuadPts * tensorSize;
- if (0 == _initialStress)
- _initialStrain->newSection(cells, fiberDim);
- else { // reuse chart from inital stress
- const ALE::Obj<RealSection::chart_type>& chart =
- initialStressSection->getChart();
- assert(!chart.isNull());
- _initialStrain->newSection(*chart, fiberDim);
- } // else
- _initialStrain->allocate();
- _initialStrain->zero();
- } // if
- const ALE::Obj<RealSection>& initialStrainSection =
- (0 != _initialStrain) ? _initialStrain->section() : 0;
-
- // Setup databases as necessary.
- if (0 != _dbInitialStress) {
- _dbInitialStress->open();
- _dbInitialStress->queryVals(stressDBValues, tensorSize);
- } // if
- if (0 != _dbInitialStrain) {
- _dbInitialStrain->open();
- _dbInitialStrain->queryVals(strainDBValues, tensorSize);
- } // if
-
- assert(0 != _normalizer);
- const double pressureScale = _normalizer->pressureScale();
- const double lengthScale = _normalizer->pressureScale();
-
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
- c_iter != cellsEnd;
- ++c_iter) {
- // Compute geometry information for current cell
- quadrature->retrieveGeometry(*c_iter);
-
- const double_array& quadPtsNonDim = quadrature->quadPts();
- quadPtsGlobal = quadPtsNonDim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
-
- // Loop over quadrature points in cell and query database
- for (int iQuadPt=0, index=0;
- iQuadPt < numQuadPts;
- ++iQuadPt, index+=spaceDim) {
- int err = _dbInitalStress->query(&valuesQuery[0], tensorSize,
- &quadPtsGlobal[index], spaceDim, cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find initial stress state at \n"
- << "(";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[index+i];
- msg << ") in material " << _label << "\n"
- << "using spatial database '" << _dbInitialStress->label() << "'.";
- throw std::runtime_error(msg.str());
- } // if
- _dbToProperties(&propertiesCell[iQuadPt*_numPropsQuadPt],
- propertiesQuery);
- _nondimProperties(&propertiesCell[iQuadPt*_numPropsQuadPt],
- _numPropsQuadPt);
-
- if (0 != _dbInitialState) {
- err = _dbInitialState->query(&stateVarsQuery[0], numDBStateVars,
- &quadPtsGlobal[index], spaceDim, cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find initial state values at \n" << "(";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[index+i];
- msg << ") in material " << _label << "\n"
- << "using spatial database '" << _dbInitialState->label()
- << "'.";
- throw std::runtime_error(msg.str());
- } // if
- _dbToStateVars(&stateVarsCell[iQuadPt*_numVarsQuadPt],
- stateVarsQuery);
- _nondimStateVars(&stateVarsCell[iQuadPt*_numVarsQuadPt],
- _numVarsQuadPt);
- } // if
-
- } // for
- // Insert cell contribution into fields
- propertiesSection->updatePoint(*c_iter, &propertiesCell[0]);
- if (0 != _dbInitialState)
- stateVarsSection->updatePoint(*c_iter, &stateVarsCell[0]);
- } // for
-
- // Close databases
- _dbProperties->close();
- if (0 != _dbInitialState)
- _dbInitialState->close();
-
-} // initialize
-
// ----------------------------------------------------------------------
// Compute density for cell at quadrature points.
const pylith::double_array&
pylith::materials::ElasticMaterial::calcDensity(void)
{ // calcDensity
const int numQuadPts = _numQuadPts;
- const int totalPropsQuadPt = _totalPropsQuadPt;
- assert(_propertiesCell.size() == numQuadPts*totalPropsQuadPt);
- assert(_density.size() == numQuadPts);
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
+ assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
+ assert(_densityCell.size() == numQuadPts*1);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- _calcDensity(&_density[iQuad],
- &_propertiesCell[iQuad*totalPropsQuadPt], totalPropsQuadPt);
+ _calcDensity(&_densityCell[iQuad],
+ &_propertiesCell[iQuad*numPropsQuadPt], numPropsQuadPt,
+ &_stateVarsCell[iQuad*numVarsQuadPt], numVarsQuadPt);
- return _density;
+ return _densityCell;
} // calcDensity
// ----------------------------------------------------------------------
@@ -231,20 +125,26 @@
const bool computeStateVars)
{ // calcStress
const int numQuadPts = _numQuadPts;
- const int totalPropsQuadPt = _totalPropsQuadPt;
- assert(_propertiesCell.size() == numQuadPts*totalPropsQuadPt);
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int tensorSize = _tensorSize;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
+ assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
+ assert(_stressCell.size() == numQuadPts*_tensorSize);
+ assert(_initialStressCell.size() == numQuadPts*_tensorSize);
+ assert(_initialStrainCell.size() == numQuadPts*_tensorSize);
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,
+ _calcStress(&_stressCell[iQuad*_tensorSize], _tensorSize,
+ &_propertiesCell[iQuad*numPropsQuadPt], numPropsQuadPt,
+ &_stateVarsCell[iQuad*numVarsQuadPt], numVarsQuadPt,
&totalStrain[iQuad*_tensorSize], _tensorSize,
- &_initialStateCell[iQuad*_initialStateSize],
- _initialStateSize, computeStateVars);
+ &_initialStressCell[iQuad*_tensorSize], _tensorSize,
+ &_initialStrainCell[iQuad*_tensorSize], _tensorSize,
+ computeStateVars);
- return _stress;
+ return _stressCell;
} // calcStress
// ----------------------------------------------------------------------
@@ -254,126 +154,371 @@
const double_array& totalStrain)
{ // calcDerivElastic
const int numQuadPts = _numQuadPts;
- const int totalPropsQuadPt = _totalPropsQuadPt;
- assert(_propertiesCell.size() == numQuadPts*totalPropsQuadPt);
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int tensorSize = _tensorSize;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
+ assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
+ assert(_elasticConstsCell.size() == numQuadPts*_numElasticConsts);
+ assert(_initialStressCell.size() == numQuadPts*_tensorSize);
+ assert(_initialStrainCell.size() == numQuadPts*_tensorSize);
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],
+ _calcElasticConsts(&_elasticConstsCell[iQuad*_numElasticConsts],
_numElasticConsts,
- &_propertiesCell[iQuad*totalPropsQuadPt],
- totalPropsQuadPt,
+ &_propertiesCell[iQuad*numPropsQuadPt],
+ numPropsQuadPt,
+ &_stateVarsCell[iQuad*numVarsQuadPt], numVarsQuadPt,
&totalStrain[iQuad*_tensorSize], _tensorSize,
- &_initialStateCell[iQuad*_initialStateSize],
- _initialStateSize);
+ &_initialStressCell[iQuad*_tensorSize], _tensorSize,
+ &_initialStrainCell[iQuad*_tensorSize], _tensorSize);
- return _elasticConsts;
+ return _elasticConstsCell;
} // calcDerivElastic
// ----------------------------------------------------------------------
+// Update state variables (for next time step).
+void
+pylith::materials::ElasticMaterial::updateStateVars(
+ const double_array& totalStrain,
+ const int cell)
+{ // updateStateVars
+ const int numQuadPts = _numQuadPts;
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int tensorSize = _tensorSize;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
+ assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
+ assert(_initialStressCell.size() == numQuadPts*_tensorSize);
+ assert(_initialStrainCell.size() == numQuadPts*_tensorSize);
+ assert(totalStrain.size() == numQuadPts*_tensorSize);
+
+ const ALE::Obj<RealSection>& stateVarsSection = _stateVars->section();
+ assert(!stateVarsSection.isNull());
+ stateVarsSection->restrictPoint(cell, &_stateVarsCell[0],
+ _stateVarsCell.size());
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ _updateStateVars(&_stateVarsCell[iQuad*numVarsQuadPt], numVarsQuadPt,
+ &_propertiesCell[iQuad*numPropsQuadPt],
+ numPropsQuadPt,
+ &totalStrain[iQuad*_tensorSize], _tensorSize,
+ &_initialStressCell[iQuad*_tensorSize], _tensorSize,
+ &_initialStrainCell[iQuad*_tensorSize], _tensorSize);
+
+ stateVarsSection->updatePoint(cell, &_stateVarsCell[0]);
+} // updateStateVars
+
+// ----------------------------------------------------------------------
// Get stable time step for implicit time integration.
double
pylith::materials::ElasticMaterial::stableTimeStepImplicit(void) const
{ // stableTimeStepImplicit
const int numQuadPts = _numQuadPts;
- const int totalPropsQuadPt = _totalPropsQuadPt;
- assert(_propertiesCell.size() == numQuadPts*totalPropsQuadPt);
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int tensorSize = _tensorSize;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
+ assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
+ assert(_elasticConstsCell.size() == numQuadPts*_numElasticConsts);
+ assert(_initialStressCell.size() == numQuadPts*_tensorSize);
+ assert(_initialStrainCell.size() == numQuadPts*_tensorSize);
double dtStable = pylith::PYLITH_MAXDOUBLE;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double dt =
- _stableTimeStepImplicit(&_propertiesCell[iQuad*totalPropsQuadPt],
- totalPropsQuadPt);
+ _stableTimeStepImplicit(&_propertiesCell[iQuad*numPropsQuadPt],
+ numPropsQuadPt,
+ &_stateVarsCell[iQuad*numVarsQuadPt],
+ numVarsQuadPt);
if (dt < dtStable)
dtStable = dt;
} // for
-
+
return dtStable;
} // stableTimeStepImplicit
// ----------------------------------------------------------------------
-// Get cell's property information from material's sections.
+// Allocate cell arrays.
void
-pylith::materials::ElasticMaterial::getPropertiesCell(const Mesh::point_type& cell,
- const int numQuadPts)
-{ // getPropertiesCell
- if (_numQuadPts != numQuadPts) {
- const int totalPropsQuadPt = _totalPropsQuadPt;
+pylith::materials::ElasticMaterial::_allocateCellArrays(void)
+{ // _allocateCellArrays
+ const int numQuadPts = _numQuadPts;
+ const int tensorSize = _tensorSize;
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ const int numElasticConsts = _numElasticConsts;
- _numQuadPts = numQuadPts;
+ _propertiesCell.resize(numQuadPts * numPropsQuadPt);
+ _stateVarsCell.resize(numQuadPts * numVarsQuadPt);
+ _initialStressCell.resize(numQuadPts * tensorSize);
+ _initialStrainCell.resize(numQuadPts * tensorSize);
+ _densityCell.resize(numQuadPts);
+ _stressCell.resize(numQuadPts * tensorSize);
+ _elasticConstsCell.resize(numQuadPts * numElasticConsts);
+} // _allocateCellArrays
- _propertiesCell.resize(numQuadPts*totalPropsQuadPt);
- _density.resize(numQuadPts*1);
- _stress.resize(numQuadPts*_tensorSize);
- _elasticConsts.resize(numQuadPts*_numElasticConsts);
- _initialStateCell.resize(numQuadPts*_tensorSize);
- } // if
-
- _getProperties(cell);
-} // getPropertiesCell
-
// ----------------------------------------------------------------------
-// Update properties (state variables for next time step).
+// Initialize initial stress field.
void
-pylith::materials::ElasticMaterial::updateProperties(
- const double_array& totalStrain,
- const Mesh::point_type& cell)
-{ // updateProperties
- const int numQuadPts = _numQuadPts;
- const int totalPropsQuadPt = _totalPropsQuadPt;
- getPropertiesCell(cell, numQuadPts);
+pylith::materials::ElasticMaterial::_initializeInitialStress(
+ const topology::Mesh& mesh,
+ feassemble::Quadrature<topology::Mesh>* quadrature)
+{ // _initializeInitialStress
+ delete _initialStress; _initialStress = 0;
+ if (0 == _dbInitialStress)
+ return;
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- _updateProperties(&_propertiesCell[iQuad*totalPropsQuadPt],
- totalPropsQuadPt,
- &totalStrain[iQuad*_tensorSize], _tensorSize,
- &_initialStateCell[iQuad*_initialStateSize],
- _initialStateSize);
+ assert(0 != _dbInitialStress);
+ assert(0 != quadrature);
+
+ // Get quadrature information
+ const int numQuadPts = quadrature->numQuadPts();
+ const int spaceDim = quadrature->spaceDim();
+
+ // Create arrays for querying
+ const int tensorSize = _tensorSize;
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+ double_array stressCell(numQuadPts*tensorSize);
+
+ // Get cells associated with material
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+
+ // Create field to hold initial stress state.
+ delete _initialStress; _initialStress = 0;
+ _initialStress = new topology::Field<topology::Mesh>(mesh);
+ assert(0 != _initialStress);
+ const int fiberDim = numQuadPts * tensorSize;
+ assert(fiberDim > 0);
+ _initialStress->newSection(cells, fiberDim);
+ _initialStress->allocate();
+ _initialStress->zero();
+ const ALE::Obj<RealSection>& initialStressSection =
+ _initialStress->section();
+
+ // Setup databases for querying
+ _dbInitialStress->open();
+ switch (dimension())
+ { // switch
+ case 1: {
+ const char* stressDBValues[] = { "stress" };
+ _dbInitialStress->queryVals(stressDBValues, tensorSize);
+ } // case 1
+ case 2 : {
+ const char* stressDBValues[] = {
+ "stress-xx",
+ "stress-yy",
+ "stress-xy",
+ };
+ _dbInitialStress->queryVals(stressDBValues, tensorSize);
+ } // case 3
+ case 3 : {
+ const char* stressDBValues[] = {
+ "stress-xx",
+ "stress-yy",
+ "stress-zz",
+ "stress-xy",
+ "stress-yz",
+ "stress-xz",
+ };
+ _dbInitialStress->queryVals(stressDBValues, tensorSize);
+ } // case 3
+ default :
+ assert(0);
+ } // switch
- _properties->updatePoint(cell, &_propertiesCell[0]);
-} // updateProperties
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->pressureScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const ALE::Obj<RealSection>& stressSection = _initialStress->section();
+ assert(!stressSection.isNull());
+
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ quadrature->retrieveGeometry(*c_iter);
+
+ // Dimensionalize coordinates for querying
+ const double_array& quadPtsNonDim = quadrature->quadPts();
+ quadPtsGlobal = quadPtsNonDim;
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
+
+ // Loop over quadrature points in cell and query database
+ for (int iQuadPt=0, iCoord=0, iStress=0;
+ iQuadPt < numQuadPts;
+ ++iQuadPt, iCoord+=spaceDim, iStress+=tensorSize) {
+ int err = _dbInitialStress->query(&stressCell[iStress], tensorSize,
+ &quadPtsGlobal[iCoord], spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find initial stress at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPtsGlobal[iCoord+i];
+ msg << ") in material " << label() << "\n"
+ << "using spatial database '" << _dbInitialStress->label() << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // for
+
+ // Nondimensionalize stress
+ _normalizer->nondimensionalize(&stressCell[0], stressCell.size(),
+ pressureScale);
+
+ stressSection->updatePoint(*c_iter, &stressCell[0]);
+ } // for
+
+ // Close databases
+ _dbInitialStress->close();
+} // _initializeInitialStress
+
// ----------------------------------------------------------------------
-// Get properties for cell.
+// Initialize initial strain field.
void
-pylith::materials::ElasticMaterial::_getProperties(const Mesh::point_type& cell)
-{ // _getProperties
- assert(!_properties.isNull());
+pylith::materials::ElasticMaterial::_initializeInitialStrain(
+ const topology::Mesh& mesh,
+ feassemble::Quadrature<topology::Mesh>* quadrature)
+{ // _initializeInitialStrain
+ delete _initialStrain; _initialStrain = 0;
+ if (0 == _dbInitialStrain)
+ return;
- const int numQuadPts = _numQuadPts;
- 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 (_initialState.isNull())
- _initialStateCell = 0.0;
- else {
- 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
+ assert(0 != _dbInitialStrain);
+ assert(0 != quadrature);
-} // _getProperties
+ // Get quadrature information
+ const int numQuadPts = quadrature->numQuadPts();
+ const int spaceDim = quadrature->spaceDim();
+ // Create arrays for querying
+ const int tensorSize = _tensorSize;
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+
+ // Get cells associated with material
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+
+ // Create field to hold initial strain state.
+ delete _initialStrain; _initialStrain = 0;
+ _initialStrain = new topology::Field<topology::Mesh>(mesh);
+ assert(0 != _initialStrain);
+ const int fiberDim = numQuadPts * tensorSize;
+ assert(fiberDim > 0);
+ _initialStrain->newSection(cells, fiberDim);
+ _initialStrain->allocate();
+ _initialStrain->zero();
+ const ALE::Obj<RealSection>& initialStrainSection =
+ _initialStrain->section();
+
+ // Setup databases for querying
+ _dbInitialStrain->open();
+ switch (dimension())
+ { // switch
+ case 1: {
+ const char* strainDBValues[] = { "strain" };
+ _dbInitialStrain->queryVals(strainDBValues, tensorSize);
+ } // case 1
+ case 2 : {
+ const char* strainDBValues[] = {
+ "strain-xx",
+ "strain-yy",
+ "strain-xy",
+ };
+ _dbInitialStrain->queryVals(strainDBValues, tensorSize);
+ } // case 3
+ case 3 : {
+ const char* strainDBValues[] = {
+ "strain-xx",
+ "strain-yy",
+ "strain-zz",
+ "strain-xy",
+ "strain-yz",
+ "strain-xz",
+ };
+ _dbInitialStrain->queryVals(strainDBValues, tensorSize);
+ } // case 3
+ default :
+ assert(0);
+ } // switch
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->pressureScale();
+ const double pressureScale = _normalizer->pressureScale();
+
+ const ALE::Obj<RealSection>& strainSection = _initialStrain->section();
+ assert(!strainSection.isNull());
+
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ quadrature->retrieveGeometry(*c_iter);
+
+ // Dimensionalize coordinates for querying
+ const double_array& quadPtsNonDim = quadrature->quadPts();
+ quadPtsGlobal = quadPtsNonDim;
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
+
+ // Loop over quadrature points in cell and query database
+ for (int iQuadPt=0, iCoord=0, iStrain=0;
+ iQuadPt < numQuadPts;
+ ++iQuadPt, iCoord+=spaceDim, iStrain+=tensorSize) {
+ int err = _dbInitialStrain->query(&strainCell[iStrain], tensorSize,
+ &quadPtsGlobal[iCoord], spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find initial strain at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPtsGlobal[iCoord+i];
+ msg << ") in material " << label() << "\n"
+ << "using spatial database '" << _dbInitialStrain->label() << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // for
+
+ // Nondimensionalize strain
+ _normalizer->nondimensionalize(&strainCell[0], strainCell.size(),
+ pressureScale);
+
+ strainSection->updatePoint(*c_iter, &strainCell[0]);
+ } // for
+
+ // Close databases
+ _dbInitialStrain->close();
+} // _initializeInitialStrain
+
// ----------------------------------------------------------------------
-// Update properties (for next time step).
+// Update stateVars (for next time step).
void
-pylith::materials::ElasticMaterial::_updateProperties(double* const properties,
- const int totalPropsQuadPt,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _updateProperties
-} // _updateProperties
+pylith::materials::ElasticMaterial::_updateStateVars(
+ double* const stateVars,
+ const int numStateVars,
+ double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVars
+} // _updateStateVars
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh 2009-02-20 23:29:29 UTC (rev 14110)
@@ -38,10 +38,10 @@
* @param numElasticConsts Number of elastic constants.
* @param metadata Metadata for physical properties and state variables.
*/
- Material(const int dimension,
- const int tensorSize,
- const int numElasticConsts,
- const Metadata& metadata);
+ ElasticMaterial(const int dimension,
+ const int tensorSize,
+ const int numElasticConsts,
+ const Metadata& metadata);
/// Destructor.
virtual
@@ -51,13 +51,13 @@
*
* @param db Spatial database.
*/
- void dbInitialStress(const spatialdata::spatialdb::SpatialDB* db);
+ void dbInitialStress(spatialdata::spatialdb::SpatialDB* db);
/** Set database for initial strain state.
*
* @param db Spatial database.
*/
- void dbInitialStrain(const spatialdata::spatialdb::SpatialDB* db);
+ void dbInitialStrain(spatialdata::spatialdb::SpatialDB* db);
/** Initialize material by getting physical property parameters from
* database.
@@ -69,6 +69,13 @@
void initialize(const topology::Mesh& mesh,
feassemble::Quadrature<topology::Mesh>* quadrature);
+ /** Retrieve parameters for physical properties and state variables
+ * for cell.
+ *
+ * @param cell Finite-element cell
+ */
+ void retrievePropsAndVars(const int cell);
+
/** Compute density for cell at quadrature points.
*
* @returns Array of density values at cell's quadrature points.
@@ -127,6 +134,22 @@
const double_array&
calcDerivElastic(const double_array& totalStrain);
+ /** Update state variables (for next time step).
+ *
+ * @param totalStrain Total strain tensor at quadrature points
+ * [numQuadPts][tensorSize]
+ * @param cell Finite element cell
+ */
+ void updateStateVars(const double_array& totalStrain,
+ const int cell);
+
+ /** Get flag indicating whether material implements an empty
+ * _updateProperties() method.
+ *
+ * @returns False if _updateProperties() is empty, true otherwise.
+ */
+ bool hasStateVars(void) const;
+
/** Get stable time step for implicit time integration.
*
* Default is MAXFLOAT (or 1.0e+30 if MAXFLOAT is not defined in math.h).
@@ -136,15 +159,6 @@
virtual
double stableTimeStepImplicit(void) const;
- /** Update properties (for next time step).
- *
- * @param totalStrain Total strain tensor at quadrature points
- * [numQuadPts][tensorSize]
- * @param cell Finite element cell
- */
- void updateProperties(const double_array& totalStrain,
- const Mesh::point_type& cell);
-
/** Set whether elastic or inelastic constitutive relations are used.
*
* @param flag True to use elastic, false to use inelastic.
@@ -152,22 +166,11 @@
virtual
void useElasticBehavior(const bool flag);
- /** Get flag indicating whether material implements an empty
- * _updateProperties() method.
- *
- * @returns False if _updateProperties() is empty, true otherwise.
- */
- virtual
- bool usesUpdateProperties(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /** Allocate cell arrays.
- *
- * @param numQuadPts Number of quadrature points.
- */
- void _allocateCellArrays(void);
+ /// These methods must be implemented by every elasticity
+ /// constitutive model.
/** Compute density from properties.
*
@@ -246,6 +249,31 @@
const double* initialStrain,
const int initialStrainSize) = 0;
+ /** Update state variables (for next time step).
+ *
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ virtual
+ void _updateStateVars(double* const stateVars,
+ const int numStateVars,
+ double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
/** Get stable time step for implicit time integration.
*
* @param properties Properties at location.
@@ -261,48 +289,39 @@
const double* stateVars,
const int numStateVars) const = 0;
- /** Update state variables (for next time step).
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Allocate cell arrays.
*
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param totalStrain Total strain at location.
- * @param strainSize Size of strain tensor.
- * @param initialStress Initial stress tensor at location.
- * @param initialStressSize Size of initial stress array.
- * @param initialStrain Initial strain tensor at location.
- * @param initialStrainSize Size of initial strain array.
+ * @param numQuadPts Number of quadrature points.
*/
- virtual
- void _updateProperties(double* const stateVars,
- const int numStateVars,
- double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialStress,
- const int initialStressSize,
- const double* initialStrain,
- const int initialStrainSize);
+ void _allocateCellArrays(void);
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
+ /** Initialize initial stress field.
+ *
+ * @param mesh Finite-element mesh.
+ * @param quadrature Quadrature for finite-element integration
+ */
+ void _initializeInitialStress(const topology::Mesh& mesh,
+ feassemble::Quadrature<topology::Mesh>* quadrature);
- /** Get properties for cell.
+ /** Initialize initial strain field.
*
- * @param cell Finite-element cell
+ * @param mesh Finite-element mesh.
+ * @param quadrature Quadrature for finite-element integration
*/
- void _getProperties(const Mesh::point_type& cell);
+ void _initializeInitialStrain(const topology::Mesh& mesh,
+ feassemble::Quadrature<topology::Mesh>* quadrature);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
/// Database for initial stress tensor;
- spatialdata::spatialdb:SpatialDB* _dbInitialStress;
+ spatialdata::spatialdb::SpatialDB* _dbInitialStress;
/// Database for initial strain tensor;
- spatialdata::spatialdb:SpatialDB* _dbInitialStrain;
+ spatialdata::spatialdb::SpatialDB* _dbInitialStrain;
/// Initial stress field.
topology::Field<topology::Mesh>* _initialStress;
@@ -343,24 +362,24 @@
* size = numQuadPts
* index = iQuadPt
*/
- double_array _density;
+ double_array _densityCell;
/** Stress tensor at quadrature points for current cell.
*
* size = numQuadPts * tensorSize
* index = iQuadPt * tensorSize + iStress
*/
- double_array _stress;
+ double_array _stressCell;
/** Elasticity matrix at quadrature points for current cell.
*
* size = numQuadPts * numElasticConsts
* index = iQuadPt * numElasticConsts + iConstant
*/
- double_array _elasticConsts;
+ double_array _elasticConstsCell;
int _numQuadPts; ///< Number of quadrature points
- int _numElasticConsts; ///< Number of elastic constants.
+ const int _numElasticConsts; ///< Number of elastic constants.
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.icc 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.icc 2009-02-20 23:29:29 UTC (rev 14110)
@@ -17,14 +17,14 @@
// Set database for initial stress state.
inline
void
-pylith::materials::ElasticMaterial::dbInitialStress(const spatialdata::spatialdb::SpatialDB* db) {
+pylith::materials::ElasticMaterial::dbInitialStress(spatialdata::spatialdb::SpatialDB* db) {
_dbInitialStress = db;
}
// Set database for initial strain state.
inline
void
-pylith::materials::ElasticMaterial::dbInitialStrain(const spatialdata::spatialdb::SpatialDB* db) {
+pylith::materials::ElasticMaterial::dbInitialStrain(spatialdata::spatialdb::SpatialDB* db) {
_dbInitialStrain = db;
}
@@ -38,8 +38,8 @@
// _updateProperties() method.
inline
bool
-pylith::materials::ElasticMaterial::usesUpdateProperties(void) const {
- return false;
+pylith::materials::ElasticMaterial::hasStateVars(void) const {
+ return _numVarsQuadPt > 0;
} // usesUpdateProperties
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.cc 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.cc 2009-02-20 23:29:29 UTC (rev 14110)
@@ -170,7 +170,7 @@
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
@@ -207,7 +207,7 @@
&quadPtsGlobal[index], spaceDim, cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find initial state values at \n" << "(";
+ msg << "Could not find initial state variables at \n" << "(";
for (int i=0; i < spaceDim; ++i)
msg << " " << quadPtsGlobal[index+i];
msg << ") in material " << _label << "\n"
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.hh 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.hh 2009-02-20 23:29:29 UTC (rev 14110)
@@ -156,6 +156,8 @@
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
+ /// These methods should be implemented by every constitutive model.
+
/** Compute properties from values in spatial database.
*
* @param propValues Array of property values.
@@ -190,7 +192,7 @@
*/
virtual
void _dbToStateVars(double* const stateValues,
- const double_array& dbValues) const = 0;
+ const double_array& dbValues) const;
/** Nondimensionalize state variables.
*
@@ -199,7 +201,7 @@
*/
virtual
void _nondimStateVars(double* const values,
- const int nvalues) const = 0;
+ const int nvalues) const;
/** Dimensionalize state variables.
*
@@ -208,7 +210,7 @@
*/
virtual
void _dimStateVars(double* const values,
- const int nvalues) const = 0;
+ const int nvalues) const;
// PROTECTED MEMBERS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.icc 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/Material.icc 2009-02-20 23:29:29 UTC (rev 14110)
@@ -93,5 +93,26 @@
_needNewJacobian = false;
} // resetNeedNewJacobian
+// Compute initial state variables from values in spatial database.
+inline
+void
+pylith::materials::Material::_dbToStateVars(double* const stateValues,
+ const double_array& dbValues) const
+{}
+// Nondimensionalize state variables.
+inline
+void
+pylith::materials::Material::_nondimStateVars(double* const values,
+ const int nvalues) const
+{}
+
+// Dimensionalize state variables.
+inline
+void
+pylith::materials::Material::_dimStateVars(double* const values,
+ const int nvalues) const
+{}
+
+
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.cc 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.cc 2009-02-20 23:29:29 UTC (rev 14110)
@@ -20,14 +20,14 @@
// ----------------------------------------------------------------------
// Constructor.
-pylith::materials::Metadata::Metadata(const int numProps,
- const ParamDescription* props,
+pylith::materials::Metadata::Metadata(const ParamDescription* props,
+ const int numProps,
+ const char* dbProps[],
const int numDBProps,
- const char* dbProps[],
+ const ParamDescription* vars,
const int numVars,
- const ParamDescription* vars,
- const int numDBVars,
- const char* dbVars[]) :
+ const char* dbVars[],
+ const int numDBVars) :
_numDBProperties(numDBProps),
_dbProperties(dbProps),
_numDBStateVars(numDBVars),
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.hh 2009-02-20 14:00:52 UTC (rev 14109)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/Metadata.hh 2009-02-20 23:29:29 UTC (rev 14110)
@@ -57,23 +57,23 @@
/** Constructor.
*
+ * @param properties Array of property descriptions.
* @param numProperties Number of physical properties for material.
- * @param properties Array of property descriptions.
+ * @param dbProperties Array of names for database values for properties.
* @param numDBProperties Number of database values for physical properties.
- * @param dbProperties Array of names for database values for properties.
+ * @param stateVars Array of state variable descriptions.
* @param numStateVars Number of state variables for material.
- * @param stateVars Array of state variable descriptions.
+ * @param dbStateVars Array of names for database values for state variables.
* @param numDBStateVars Number of database values for state variables.
- * @param dbStateVars Array of names for database values for state variables.
*/
- Metadata(const int numProperties,
- const ParamDescription* properties,
+ Metadata(const ParamDescription* properties,
+ const int numProperties,
+ const char* dbProperties[],
const int numDBProperties,
- const char* dbProperties[],
+ const ParamDescription* stateVars,
const int numStateVars,
- const ParamDescription* stateVars,
- const int numDBStateVars,
- const char* dbStateVars[]);
+ const char* dbStateVars[],
+ const int numDBStateVars);
/// Default destructor
~Metadata(void);
More information about the CIG-COMMITS
mailing list