[cig-commits] r13808 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Jan 7 20:28:29 PST 2009
Author: willic3
Date: 2009-01-07 20:28:29 -0800 (Wed, 07 Jan 2009)
New Revision: 13808
Added:
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc
Log:
Incomplete files for power-law material.
Added: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-01-08 04:28:29 UTC (rev 13808)
@@ -0,0 +1,652 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "PowerLaw3D.hh" // implementation of object methods
+
+#include "pylith/utils/array.hh" // USES double_array
+
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include "petsc.h" // USES PetscLogFlops
+
+#include <cassert> // USES assert()
+#include <cstring> // USES memcpy()
+#include <sstream> // USES std::ostringstream
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+namespace pylith {
+ namespace materials {
+ namespace _PowerLaw3D{
+
+ /// Number of entries in stress/strain tensors.
+ const int tensorSize = 6;
+
+ /// Number of entries in derivative of elasticity matrix.
+ const int numElasticConsts = 21;
+
+ /// Number of physical properties.
+ const int numProperties = 7;
+
+ /// Physical properties.
+ const Material::PropMetaData properties[] = {
+ { "density", 1, OTHER_FIELD },
+ { "mu", 1, OTHER_FIELD },
+ { "lambda", 1, OTHER_FIELD },
+ { "viscosity_coeff", 1, OTHER_FIELD },
+ { "power_law_exponent", 1, OTHER_FIELD },
+ { "total_strain", 6, OTHER_FIELD },
+ { "viscous_strain", 6, OTHER_FIELD },
+ };
+ /// Indices (order) of properties.
+ const int pidDensity = 0;
+ const int pidMu = pidDensity + 1;
+ const int pidLambda = pidMu + 1;
+ const int pidViscosityCoeff = pidLambda + 1;
+ const int pidPowerLawExp = pidViscosityCoeff + 1;
+ const int pidStrainT = pidPowerLawExp + 1;
+ const int pidVisStrain = pidStrainT + tensorSize;
+
+ /// Values expected in spatial database
+ const int numDBValues = 5;
+ const char* namesDBValues[] = {"density", "vs", "vp" , "viscosity_coeff",
+ "power_law_exponent"};
+
+ /// Indices (order) of database values
+ const int didDensity = 0;
+ const int didVs = 1;
+ const int didVp = 2;
+ const int didViscosityCoeff = 3;
+ const int didPowerLawExp = 4;
+
+ /// 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" };
+
+ } // _PowerLaw3D
+ } // materials
+} // pylith
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::materials::PowerLaw3D::PowerLaw3D(void) :
+ ElasticMaterial(_PowerLaw3D::tensorSize,
+ _PowerLaw3D::numElasticConsts,
+ _PowerLaw3D::namesDBValues,
+ _PowerLaw3D::namesInitialStateDBValues,
+ _PowerLaw3D::numDBValues,
+ _PowerLaw3D::properties,
+ _PowerLaw3D::numProperties),
+ _calcElasticConstsFn(&pylith::materials::PowerLaw3D::_calcElasticConstsElastic),
+ _calcStressFn(&pylith::materials::PowerLaw3D::_calcStressElastic),
+ _updatePropertiesFn(&pylith::materials::PowerLaw3D::_updatePropertiesElastic)
+{ // constructor
+ _dimension = 3;
+ _visStrain.resize(_PowerLaw3D::tensorSize);
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::materials::PowerLaw3D::~PowerLaw3D(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Compute properties from values in spatial database.
+void
+pylith::materials::PowerLaw3D::_dbToProperties(
+ double* const propValues,
+ const double_array& dbValues) const
+{ // _dbToProperties
+ assert(0 != propValues);
+ const int numDBValues = dbValues.size();
+ assert(_PowerLaw3D::numDBValues == numDBValues);
+
+ const double density = dbValues[_PowerLaw3D::didDensity];
+ const double vs = dbValues[_PowerLaw3D::didVs];
+ const double vp = dbValues[_PowerLaw3D::didVp];
+ const double viscosityCoeff = dbValues[_PowerLaw3D::didViscosityCoeff];
+ const double powerLawExp = dbValues[_PowerLaw3D::didPowerLawExp];
+
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || viscosityCoeff <= 0.0
+ || powerLawExp < 1.0) {
+ std::ostringstream msg;
+ msg << "Spatial database returned illegal value for physical "
+ << "properties.\n"
+ << "density: " << density << "\n"
+ << "vp: " << vp << "\n"
+ << "vs: " << vs << "\n"
+ << "viscosityCoeff: " << viscosityCoeff << "\n"
+ << "powerLawExp: " << powerLawExp << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ const double mu = density * vs*vs;
+ const double lambda = density * vp*vp - 2.0*mu;
+
+ if (lambda <= 0.0) {
+ std::ostringstream msg;
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
+ << "density: " << density << "\n"
+ << "vp: " << vp << "\n"
+ << "vs: " << vs << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+ assert(mu > 0);
+
+ propValues[_PowerLaw3D::pidDensity] = density;
+ propValues[_PowerLaw3D::pidMu] = mu;
+ propValues[_PowerLaw3D::pidLambda] = lambda;
+ propValues[_PowerLaw3D::pidViscosityCoeff] = viscosityCoeff;
+ propValues[_PowerLaw3D::pidPowerLawExp] = powerLawExp;
+
+ PetscLogFlops(6);
+} // _dbToProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::PowerLaw3D::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _totalPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ // **** NOTE: Make sure scaling is correct for viscosity coefficient.
+ const double powerLawExp = _PowerLaw3D::pidPowerLawExp;
+ const double viscosityCoeffScale = timeScale *
+ pressureScale^(1.0/powerLawExp);
+ const double powerLawExpScale = 1.0;
+ values[_PowerLaw3D::pidDensity] =
+ _normalizer->nondimensionalize(values[_PowerLaw3D::pidDensity],
+ densityScale);
+ values[_PowerLaw3D::pidMu] =
+ _normalizer->nondimensionalize(values[_PowerLaw3D::pidMu],
+ pressureScale);
+ values[_PowerLaw3D::pidLambda] =
+ _normalizer->nondimensionalize(values[_PowerLaw3D::pidLambda],
+ pressureScale);
+ values[_PowerLaw3D::pidViscosityCoeff] =
+ _normalizer->nondimensionalize(values[_PowerLaw3D::pidViscosityCoeff],
+ viscosityCoeffScale);
+ values[_PowerLaw3D::pidPowerLawExp] =
+ _normalizer->nondimensionalize(values[_PowerLaw3D::pidPowerLawExp],
+ powerLawExpScale);
+
+ PetscLogFlops(8);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::PowerLaw3D::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _totalPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ // **** NOTE: Make sure scaling is correct for viscosity coefficient.
+ const double powerLawExp = _PowerLaw3D::pidPowerLawExp;
+ const double viscosityCoeffScale = timeScale *
+ pressureScale^(1.0/powerLawExp);
+ const double powerLawExpScale = 1.0;
+ values[_PowerLaw3D::pidDensity] =
+ _normalizer->dimensionalize(values[_PowerLaw3D::pidDensity],
+ densityScale);
+ values[_PowerLaw3D::pidMu] =
+ _normalizer->dimensionalize(values[_PowerLaw3D::pidMu],
+ pressureScale);
+ values[_PowerLaw3D::pidLambda] =
+ _normalizer->dimensionalize(values[_PowerLaw3D::pidLambda],
+ pressureScale);
+ values[_PowerLaw3D::pidViscosityCoeff] =
+ _normalizer->dimensionalize(values[_PowerLaw3D::pidViscosityCoeff],
+ viscosityCoeffScale);
+ values[_PowerLaw3D::pidPowerLawExp] =
+ _normalizer->dimensionalize(values[_PowerLaw3D::pidPowerLawExp],
+ powerLawExpScale);
+
+ PetscLogFlops(8);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::PowerLaw3D::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _PowerLaw3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::PowerLaw3D::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _PowerLaw3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
+// Compute density at location from properties.
+void
+pylith::materials::PowerLaw3D::_calcDensity(double* const density,
+ const double* properties,
+ const int numProperties)
+{ // _calcDensity
+ assert(0 != density);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+
+ density[0] = properties[_PowerLaw3D::pidDensity];
+} // _calcDensity
+
+// ----------------------------------------------------------------------
+// Compute viscous strain for current time step.
+// material.
+void
+pylith::materials::PowerLaw3D::_computeStateVars(
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
+{ // _computeStateVars
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(_PowerLaw3D::tensorSize == initialStateSize);
+
+ const int tensorSize = _PowerLaw3D::tensorSize;
+ const double maxwelltime = properties[_PowerLaw3D::pidMaxwellTime];
+
+ 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 meanStrainTpdt = (e11 + e22 + e33)/3.0;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
+ const double meanStrainT =
+ (properties[_PowerLaw3D::pidStrainT+0] +
+ properties[_PowerLaw3D::pidStrainT+1] +
+ properties[_PowerLaw3D::pidStrainT+2])/3.0;
+
+ // Time integration.
+ double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
+ const double expFac = exp(-_dt/maxwelltime);
+
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ devStrainT = properties[_PowerLaw3D::pidStrainT+iComp] -
+ diag[iComp] * meanStrainT;
+ _visStrain[iComp] = expFac *
+ properties[_PowerLaw3D::pidVisStrain + iComp] +
+ dq * (devStrainTpdt - devStrainT);
+ } // for
+
+ PetscLogFlops(8 + 7 * tensorSize);
+} // _computeStateVars
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties as an elastic
+// material.
+void
+pylith::materials::PowerLaw3D::_calcStressElastic(
+ double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
+ const bool computeStateVars)
+{ // _calcStressElastic
+ assert(0 != stress);
+ assert(_PowerLaw3D::tensorSize == stressSize);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(_PowerLaw3D::tensorSize == initialStateSize);
+
+ const double mu = properties[_PowerLaw3D::pidMu];
+ const double lambda = properties[_PowerLaw3D::pidLambda];
+ 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 traceStrainTpdt = e11 + e22 + e33;
+ const double s123 = lambda * traceStrainTpdt;
+
+ stress[0] = s123 + mu2*e11 + initialState[0];
+ stress[1] = s123 + mu2*e22 + initialState[1];
+ stress[2] = s123 + mu2*e33 + initialState[2];
+ stress[3] = mu2 * e12 + initialState[3];
+ stress[4] = mu2 * e23 + initialState[4];
+ stress[5] = mu2 * e13 + initialState[5];
+
+ PetscLogFlops(19);
+} // _calcStressElastic
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties as a viscoelastic
+// material.
+void
+pylith::materials::PowerLaw3D::_calcStressViscoelastic(
+ double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize,
+ const bool computeStateVars)
+{ // _calcStressViscoelastic
+ assert(0 != stress);
+ assert(_PowerLaw3D::tensorSize == stressSize);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(_PowerLaw3D::tensorSize == initialStateSize);
+
+ const int tensorSize = _PowerLaw3D::tensorSize;
+
+ const double mu = properties[_PowerLaw3D::pidMu];
+ const double lambda = properties[_PowerLaw3D::pidLambda];
+ const double viscosityCoeff = properties[_PowerLaw3D::pidViscosityCoeff];
+ const double powerLawExp = properties[_PowerLaw3D::pidPowerLawExp];
+
+ const double mu2 = 2.0 * mu;
+ const double bulkModulus = lambda + mu2/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 meanStressTpdt = bulkModulus * traceStrainTpdt;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
+ // Get viscous strains
+ if (computeStateVars) {
+ pylith::materials::PowerLaw3D::_computeStateVars(properties,
+ numProperties,
+ totalStrain,
+ strainSize,
+ initialState,
+ initialStateSize);
+ } else {
+ memcpy(&_visStrain[0], &properties[_PowerLaw3D::pidVisStrain],
+ tensorSize * sizeof(double));
+ } // else
+
+ // Compute new stresses
+ double devStressTpdt = 0.0;
+
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ devStressTpdt = mu2 * _visStrain[iComp];
+
+ // Later I will want to put in initial stresses.
+ stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
+ initialState[iComp];
+ } // for
+
+ PetscLogFlops(7 + 4 * tensorSize);
+} // _calcStressViscoelastic
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from properties.
+void
+pylith::materials::PowerLaw3D::_calcElasticConstsElastic(
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
+{ // _calcElasticConstsElastic
+ assert(0 != elasticConsts);
+ assert(_PowerLaw3D::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(_PowerLaw3D::tensorSize == initialStateSize);
+
+ const double mu = properties[_PowerLaw3D::pidMu];
+ const double lambda = properties[_PowerLaw3D::pidLambda];
+
+ const double mu2 = 2.0 * mu;
+ const double lambda2mu = lambda + mu2;
+
+ 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] = mu2; // C1212
+ elasticConsts[16] = 0; // C1223
+ elasticConsts[17] = 0; // C1213
+ elasticConsts[18] = mu2; // C2323
+ elasticConsts[19] = 0; // C2313
+ elasticConsts[20] = mu2; // C1313
+
+ PetscLogFlops(4);
+} // _calcElasticConstsElastic
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from properties
+// as an elastic material.
+void
+pylith::materials::PowerLaw3D::_calcElasticConstsViscoelastic(
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
+{ // _calcElasticConstsViscoelastic
+ assert(0 != elasticConsts);
+ assert(_PowerLaw3D::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(_PowerLaw3D::tensorSize == initialStateSize);
+
+ const double mu = properties[_PowerLaw3D::pidMu];
+ const double lambda = properties[_PowerLaw3D::pidLambda];
+ const double maxwelltime = properties[_PowerLaw3D::pidMaxwellTime];
+
+ const double mu2 = 2.0 * mu;
+ const double bulkModulus = lambda + mu2/3.0;
+
+ double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
+
+ const double visFac = mu*dq/3.0;
+ elasticConsts[ 0] = bulkModulus + 4.0*visFac; // C1111
+ elasticConsts[ 1] = bulkModulus - 2.0*visFac; // C1122
+ elasticConsts[ 2] = elasticConsts[1]; // C1133
+ 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] = 6.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
+
+ PetscLogFlops(10);
+} // _calcElasticConstsViscoelastic
+
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+pylith::materials::PowerLaw3D::_stableTimeStepImplicit(const double* properties,
+ const int numProperties) const
+{ // _stableTimeStepImplicit
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+
+ const double maxwellTime =
+ properties[_PowerLaw3D::pidMaxwellTime];
+ const double dtStable = 0.1*maxwellTime;
+
+ return dtStable;
+} // _stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::PowerLaw3D::_updatePropertiesElastic(
+ double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
+{ // _updatePropertiesElastic
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+
+ const double maxwelltime = properties[_PowerLaw3D::pidMaxwellTime];
+
+ 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 diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
+ for (int iComp=0; iComp < _PowerLaw3D::tensorSize; ++iComp) {
+ properties[_PowerLaw3D::pidStrainT+iComp] = totalStrain[iComp];
+ properties[_PowerLaw3D::pidVisStrain+iComp] =
+ totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ } // for
+ PetscLogFlops(3 + 2 * _PowerLaw3D::tensorSize);
+
+ _needNewJacobian = true;
+} // _updatePropertiesElastic
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::PowerLaw3D::_updatePropertiesViscoelastic(
+ double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize)
+{ // _updatePropertiesViscoelastic
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(_PowerLaw3D::tensorSize == initialStateSize);
+
+ const int tensorSize = _PowerLaw3D::tensorSize;
+
+ pylith::materials::PowerLaw3D::_computeStateVars(properties,
+ numProperties,
+ totalStrain,
+ strainSize,
+ initialState,
+ initialStateSize);
+
+ memcpy(&properties[_PowerLaw3D::pidVisStrain],
+ &_visStrain[0],
+ tensorSize * sizeof(double));
+ memcpy(&properties[_PowerLaw3D::pidStrainT],
+ &totalStrain[0],
+ tensorSize * sizeof(double));
+
+ _needNewJacobian = false;
+
+} // _updatePropertiesViscoelastic
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-01-08 04:28:29 UTC (rev 13808)
@@ -0,0 +1,397 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/MaxwellIsotropic3D.h
+ *
+ * @brief C++ MaxwellIsotropic3D object
+ *
+ * 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. The viscosity is stored using Maxwell Time (viscosity/mu).
+ */
+
+#if !defined(pylith_materials_maxwellisotropic3d_hh)
+#define pylith_materials_maxwellisotropic3d_hh
+
+#include "ElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class MaxwellIsotropic3D;
+ class TestMaxwellIsotropic3D; // unit testing
+ } // materials
+} // pylith
+
+/// 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
+ MaxwellIsotropic3D(void);
+
+ /// Destructor
+ ~MaxwellIsotropic3D(void);
+
+ /** Set current time step.
+ *
+ * @param dt Current time step.
+ */
+ void timeStep(const double dt);
+
+ /** Set whether elastic or inelastic constitutive relations are used.
+ *
+ * @param flag True to use elastic, false to use inelastic.
+ */
+ void useElasticBehavior(const bool flag);
+
+ /** Get flag indicating whether material implements an empty
+ * _updateProperties() method.
+ *
+ * @returns False if _updateProperties() is empty, true otherwise.
+ */
+ bool usesUpdateProperties(void) const;
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Compute properties from values in spatial database.
+ *
+ * Order of values in arrays matches order used in dbValues() and
+ * parameterNames().
+ *
+ * @param propValues Array of property values.
+ * @param dbValues Array of database values.
+ */
+ void _dbToProperties(double* const propValues,
+ const double_array& dbValues) const;
+
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _nondimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ 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.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ */
+ void _calcDensity(double* const density,
+ const double* properties,
+ const int numProperties);
+
+ /** 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.
+ *
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _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);
+
+ /** Compute derivatives of elasticity matrix from properties.
+ *
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _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);
+
+ /** Get stable time step for implicit time integration.
+ *
+ * @returns Time step
+ */
+ double _stableTimeStepImplicit(const double* properties,
+ const int numProperties) const;
+
+ /** Update properties (for next time step).
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _updateProperties(double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
+
+ // PRIVATE TYPEDEFS ///////////////////////////////////////////////////
+private :
+
+ /// Member prototype for _calcStress()
+ typedef void (pylith::materials::MaxwellIsotropic3D::*calcStress_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const bool);
+
+ /// Member prototype for _calcElasticConsts()
+ typedef void (pylith::materials::MaxwellIsotropic3D::*calcElasticConsts_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int);
+
+ /// Member prototype for _updateProperties()
+ typedef void (pylith::materials::MaxwellIsotropic3D::*updateProperties_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int);
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+/** Compute viscous strains (state variables) for the current time step.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _computeStateVars(const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
+
+ /** Compute stress tensor from properties as an elastic material.
+ *
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param properties Properties at locations.
+ * @param numProperties Number of properties.
+ * @param 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,
+ const int stressSize,
+ const double* properties,
+ 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.
+ *
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param properties Properties at locations.
+ * @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,
+ const int stressSize,
+ const double* properties,
+ 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
+ * elastic material.
+ *
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _calcElasticConstsElastic(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
+
+ /** Compute derivatives of elasticity matrix from properties as a
+ * viscoelastic material.
+ *
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _calcElasticConstsViscoelastic(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
+
+ /** Update state variables after solve as an elastic material.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _updatePropertiesElastic(double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
+
+ /** Update state variables after solve as a viscoelastic material.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _updatePropertiesViscoelastic(double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize);
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ MaxwellIsotropic3D(const MaxwellIsotropic3D& m);
+
+ /// Not implemented
+ const MaxwellIsotropic3D& operator=(const MaxwellIsotropic3D& m);
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ /// Viscous strain array.
+ double_array _visStrain;
+
+ /// Method to use for _calcElasticConsts().
+ calcElasticConsts_fn_type _calcElasticConstsFn;
+
+ /// Method to use for _calcStress().
+ calcStress_fn_type _calcStressFn;
+
+ /// Method to use for _updateProperties().
+ updateProperties_fn_type _updatePropertiesFn;
+
+}; // class MaxwellIsotropic3D
+
+#include "MaxwellIsotropic3D.icc" // inline methods
+
+#endif // pylith_materials_maxwellisotropic3d_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc 2009-01-08 04:28:29 UTC (rev 13808)
@@ -0,0 +1,112 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_materials_maxwellisotropic3d_hh)
+#error "MaxwellIsotropic3D.icc can only be included from MaxwellIsotropic3D.hh"
+#endif
+
+#include <assert.h> // USES assert()
+#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+
+// Set current time step.
+inline
+void
+pylith::materials::MaxwellIsotropic3D::timeStep(const double dt) {
+ // Jacobian needs to be reformed if the time step size changes.
+ if (_dt > 0.0 && dt != _dt)
+ _needNewJacobian = true;
+ _dt = dt;
+} // timeStep
+
+// Set whether elastic or inelastic constitutive relations are used.
+inline
+void
+pylith::materials::MaxwellIsotropic3D::useElasticBehavior(const bool flag) {
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsElastic;
+ _updatePropertiesFn =
+ &pylith::materials::MaxwellIsotropic3D::_updatePropertiesElastic;
+ } else {
+ _calcStressFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcStressViscoelastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsViscoelastic;
+ _updatePropertiesFn =
+ &pylith::materials::MaxwellIsotropic3D::_updatePropertiesViscoelastic;
+ } // if/else
+} // useElasticBehavior
+
+// Get flag indicating whether material implements an empty
+inline
+bool
+pylith::materials::MaxwellIsotropic3D::usesUpdateProperties(void) const {
+ return true;
+} // usesUpdateProperties
+
+// Compute stress tensor from parameters.
+inline
+void
+pylith::materials::MaxwellIsotropic3D::_calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* 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
+
+// Compute derivatives of elasticity matrix from parameters.
+inline
+void
+pylith::materials::MaxwellIsotropic3D::_calcElasticConsts(
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) {
+ assert(0 != _calcElasticConstsFn);
+ CALL_MEMBER_FN(*this, _calcElasticConstsFn)(elasticConsts, numElasticConsts,
+ parameters, numParams,
+ totalStrain, strainSize,
+ initialState, initialStateSize);
+} // _calcElasticConsts
+
+// Update state variables after solve.
+inline
+void
+pylith::materials::MaxwellIsotropic3D::_updateProperties(double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) {
+ assert(0 != _updatePropertiesFn);
+ CALL_MEMBER_FN(*this, _updatePropertiesFn)(parameters, numParams,
+ totalStrain, strainSize,
+ initialState, initialStateSize);
+} // _updateProperties
+
+// End of file
More information about the CIG-COMMITS
mailing list