[cig-commits] r14644 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Apr 8 21:54:57 PDT 2009
Author: willic3
Date: 2009-04-08 21:54:56 -0700 (Wed, 08 Apr 2009)
New Revision: 14644
Added:
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc
Modified:
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
Log:
Added initial versions of Maxwell plane strain materials.
They won't work yet because we need a new strain computation function.
Added: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2009-04-09 04:54:56 UTC (rev 14644)
@@ -0,0 +1,602 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "MaxwellPlaneStrain.hh" // implementation of object methods
+
+#include "ViscoelasticMaxwell.hh" // USES computeVisStrain
+
+#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 _MaxwellPlaneStrain{
+
+ /// Number of entries in stress/strain tensors.
+ const int tensorSize = 4;
+
+ /// Number of entries in derivative of elasticity matrix.
+ const int numElasticConsts = 6;
+
+ /// Number of physical properties.
+ const int numProperties = 6;
+
+ /// Physical properties.
+ const Material::PropMetaData properties[] = {
+ { "density", 1, OTHER_FIELD },
+ { "mu", 1, OTHER_FIELD },
+ { "lambda", 1, OTHER_FIELD },
+ { "maxwell_time", 1, OTHER_FIELD },
+ { "total_strain", 4, OTHER_FIELD },
+ { "viscous_strain", 4, OTHER_FIELD },
+ };
+ /// Indices (order) of properties.
+ const int pidDensity = 0;
+ const int pidMu = pidDensity + 1;
+ const int pidLambda = pidMu + 1;
+ const int pidMaxwellTime = pidLambda + 1;
+ const int pidStrainT = pidMaxwellTime + 1;
+ const int pidVisStrain = pidStrainT + tensorSize;
+
+ /// Values expected in spatial database
+ const int numDBValues = 4;
+ const char* namesDBValues[] = {"density", "vs", "vp" , "viscosity"};
+
+ /// Indices (order) of database values
+ const int didDensity = 0;
+ const int didVs = 1;
+ const int didVp = 2;
+ const int didViscosity = 3;
+
+ /// Initial state values expected in spatial database
+ const int numInitialStateDBValues = tensorSize;
+ const char* namesInitialStateDBValues[] = { "stress_xx", "stress_yy",
+ "stress_zz", "stress_xy" };
+
+ } // _MaxwellPlaneStrain
+ } // materials
+} // pylith
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::materials::MaxwellPlaneStrain::MaxwellPlaneStrain(void) :
+ ElasticMaterial(_MaxwellPlaneStrain::tensorSize,
+ _MaxwellPlaneStrain::numElasticConsts,
+ _MaxwellPlaneStrain::namesDBValues,
+ _MaxwellPlaneStrain::namesInitialStateDBValues,
+ _MaxwellPlaneStrain::numDBValues,
+ _MaxwellPlaneStrain::properties,
+ _MaxwellPlaneStrain::numProperties),
+ _calcElasticConstsFn(&pylith::materials::MaxwellPlaneStrain::_calcElasticConstsElastic),
+ _calcStressFn(&pylith::materials::MaxwellPlaneStrain::_calcStressElastic),
+ _updatePropertiesFn(&pylith::materials::MaxwellPlaneStrain::_updatePropertiesElastic)
+{ // constructor
+ _dimension = 2;
+ _visStrain.resize(_MaxwellPlaneStrain::tensorSize);
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::materials::MaxwellPlaneStrain::~MaxwellPlaneStrain(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Compute properties from values in spatial database.
+void
+pylith::materials::MaxwellPlaneStrain::_dbToProperties(
+ double* const propValues,
+ const double_array& dbValues) const
+{ // _dbToProperties
+ assert(0 != propValues);
+ const int numDBValues = dbValues.size();
+ assert(_MaxwellPlaneStrain::numDBValues == numDBValues);
+
+ const double density = dbValues[_MaxwellPlaneStrain::didDensity];
+ const double vs = dbValues[_MaxwellPlaneStrain::didVs];
+ const double vp = dbValues[_MaxwellPlaneStrain::didVp];
+ const double viscosity = dbValues[_MaxwellPlaneStrain::didViscosity];
+
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || viscosity <= 0.0) {
+ std::ostringstream msg;
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
+ << "density: " << density << "\n"
+ << "vp: " << vp << "\n"
+ << "vs: " << vs << "\n"
+ << "viscosity: " << viscosity << "\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);
+
+ const double maxwelltime = viscosity / mu;
+ assert(maxwelltime > 0.0);
+
+ propValues[_MaxwellPlaneStrain::pidDensity] = density;
+ propValues[_MaxwellPlaneStrain::pidMu] = mu;
+ propValues[_MaxwellPlaneStrain::pidLambda] = lambda;
+ propValues[_MaxwellPlaneStrain::pidMaxwellTime] = maxwelltime;
+
+ PetscLogFlops(7);
+} // _dbToProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::MaxwellPlaneStrain::_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();
+ values[_MaxwellPlaneStrain::pidDensity] =
+ _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidDensity],
+ densityScale);
+ values[_MaxwellPlaneStrain::pidMu] =
+ _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidMu],
+ pressureScale);
+ values[_MaxwellPlaneStrain::pidLambda] =
+ _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidLambda],
+ pressureScale);
+ values[_MaxwellPlaneStrain::pidMaxwellTime] =
+ _normalizer->nondimensionalize(values[_MaxwellPlaneStrain::pidMaxwellTime],
+ timeScale);
+
+ PetscLogFlops(4);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::MaxwellPlaneStrain::_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();
+ values[_MaxwellPlaneStrain::pidDensity] =
+ _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidDensity],
+ densityScale);
+ values[_MaxwellPlaneStrain::pidMu] =
+ _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidMu],
+ pressureScale);
+ values[_MaxwellPlaneStrain::pidLambda] =
+ _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidLambda],
+ pressureScale);
+ values[_MaxwellPlaneStrain::pidMaxwellTime] =
+ _normalizer->dimensionalize(values[_MaxwellPlaneStrain::pidMaxwellTime],
+ timeScale);
+
+ PetscLogFlops(4);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::MaxwellPlaneStrain::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _MaxwellPlaneStrain::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::MaxwellPlaneStrain::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _MaxwellPlaneStrain::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
+// Compute density at location from properties.
+void
+pylith::materials::MaxwellPlaneStrain::_calcDensity(double* const density,
+ const double* properties,
+ const int numProperties)
+{ // _calcDensity
+ assert(0 != density);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+
+ density[0] = properties[_MaxwellPlaneStrain::pidDensity];
+} // _calcDensity
+
+// ----------------------------------------------------------------------
+// Compute viscous strain for current time step.
+// material.
+void
+pylith::materials::MaxwellPlaneStrain::_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(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+
+ const int tensorSize = _MaxwellPlaneStrain::tensorSize;
+ const double maxwelltime = properties[_MaxwellPlaneStrain::pidMaxwellTime];
+
+ const double e11 = totalStrain[0];
+ const double e22 = totalStrain[1];
+ const double e33 = totalStrain[2];
+ const double e12 = totalStrain[3];
+
+ const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0 };
+
+ const double meanStrainT =
+ (properties[_MaxwellPlaneStrain::pidStrainT+0] +
+ properties[_MaxwellPlaneStrain::pidStrainT+1] +
+ properties[_MaxwellPlaneStrain::pidStrainT+2])/3.0;
+
+ // Time integration.
+ double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
+ const double expFac = exp(-_dt/maxwelltime);
+
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ devStrainT = properties[_MaxwellPlaneStrain::pidStrainT+iComp] -
+ diag[iComp] * meanStrainT;
+ _visStrain[iComp] = expFac *
+ properties[_MaxwellPlaneStrain::pidVisStrain + iComp] +
+ dq * (devStrainTpdt - devStrainT);
+ } // for
+
+ PetscLogFlops(8 + 7 * tensorSize);
+} // _computeStateVars
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties as an elastic
+// material.
+void
+pylith::materials::MaxwellPlaneStrain::_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(_MaxwellPlaneStrain::tensorSize == stressSize);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+
+ const double mu = properties[_MaxwellPlaneStrain::pidMu];
+ const double lambda = properties[_MaxwellPlaneStrain::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 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 + initialState[2];
+ stress[3] = mu2 * e12 + initialState[3];
+
+ PetscLogFlops(13);
+} // _calcStressElastic
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties as a viscoelastic
+// material.
+void
+pylith::materials::MaxwellPlaneStrain::_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(_MaxwellPlaneStrain::tensorSize == stressSize);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+
+ const int tensorSize = _MaxwellPlaneStrain::tensorSize;
+
+ const double mu = properties[_MaxwellPlaneStrain::pidMu];
+ const double lambda = properties[_MaxwellPlaneStrain::pidLambda];
+ const double maxwelltime = properties[_MaxwellPlaneStrain::pidMaxwellTime];
+
+ 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 };
+
+ // Get viscous strains
+ if (computeStateVars) {
+ pylith::materials::MaxwellPlaneStrain::_computeStateVars(properties,
+ numProperties,
+ totalStrain,
+ strainSize,
+ initialState,
+ initialStateSize);
+ } else {
+ memcpy(&_visStrain[0], &properties[_MaxwellPlaneStrain::pidVisStrain],
+ tensorSize * sizeof(double));
+ } // else
+
+ // Compute new stresses
+ double devStressTpdt = 0.0;
+
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ devStressTpdt = mu2 * _visStrain[iComp];
+
+ 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::MaxwellPlaneStrain::_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(_MaxwellPlaneStrain::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+
+ const double mu = properties[_MaxwellPlaneStrain::pidMu];
+ const double lambda = properties[_MaxwellPlaneStrain::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] = lambda2mu; // C2222
+ elasticConsts[ 5] = lambda; // C2233
+ elasticConsts[ 6] = 0; // C2212
+ elasticConsts[ 7] = lambda2mu; // C3333
+ elasticConsts[ 8] = 0; // C3312
+ elasticConsts[ 9] = mu2; // C1212
+
+ PetscLogFlops(2);
+} // _calcElasticConstsElastic
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from properties
+// as an elastic material.
+void
+pylith::materials::MaxwellPlaneStrain::_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(_MaxwellPlaneStrain::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+
+ const double mu = properties[_MaxwellPlaneStrain::pidMu];
+ const double lambda = properties[_MaxwellPlaneStrain::pidLambda];
+ const double maxwelltime = properties[_MaxwellPlaneStrain::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] = elasticConsts[0]; // C2222
+ elasticConsts[ 5] = elasticConsts[1]; // C2233
+ elasticConsts[ 6] = 0; // C2212
+ elasticConsts[ 7] = elasticConsts[0]; // C3333
+ elasticConsts[ 8] = 0; // C3312
+ elasticConsts[ 9] = 6.0 * visFac; // C1212
+
+ PetscLogFlops(10);
+} // _calcElasticConstsViscoelastic
+
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+pylith::materials::MaxwellPlaneStrain::_stableTimeStepImplicit(const double* properties,
+ const int numProperties) const
+{ // _stableTimeStepImplicit
+ assert(0 != properties);
+ assert(_totalPropsQuadPt == numProperties);
+
+ const double maxwellTime =
+ properties[_MaxwellPlaneStrain::pidMaxwellTime];
+ const double dtStable = 0.1*maxwellTime;
+
+ return dtStable;
+} // _stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::MaxwellPlaneStrain::_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(_MaxwellPlaneStrain::tensorSize == strainSize);
+
+ const double maxwelltime = properties[_MaxwellPlaneStrain::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 };
+
+ for (int iComp=0; iComp < _MaxwellPlaneStrain::tensorSize; ++iComp) {
+ properties[_MaxwellPlaneStrain::pidStrainT+iComp] = totalStrain[iComp];
+ properties[_MaxwellPlaneStrain::pidVisStrain+iComp] =
+ totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
+ } // for
+ PetscLogFlops(3 + 2 * _MaxwellPlaneStrain::tensorSize);
+
+ _needNewJacobian = true;
+} // _updatePropertiesElastic
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::MaxwellPlaneStrain::_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(_MaxwellPlaneStrain::tensorSize == strainSize);
+ assert(_MaxwellPlaneStrain::tensorSize == initialStateSize);
+
+ const int tensorSize = _MaxwellPlaneStrain::tensorSize;
+
+ pylith::materials::MaxwellPlaneStrain::_computeStateVars(properties,
+ numProperties,
+ totalStrain,
+ strainSize,
+ initialState,
+ initialStateSize);
+
+ memcpy(&properties[_MaxwellPlaneStrain::pidVisStrain],
+ &_visStrain[0],
+ tensorSize * sizeof(double));
+ memcpy(&properties[_MaxwellPlaneStrain::pidStrainT],
+ &totalStrain[0],
+ tensorSize * sizeof(double));
+
+ _needNewJacobian = false;
+
+} // _updatePropertiesViscoelastic
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.hh 2009-04-09 04:54:56 UTC (rev 14644)
@@ -0,0 +1,397 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/MaxwellPlaneStrain.h
+ *
+ * @brief C++ MaxwellPlaneStrain object
+ *
+ * 2-D, isotropic, linear Maxwell viscoelastic material for plane strain. 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_maxwellplanestrain_hh)
+#define pylith_materials_maxwellplanestrain_hh
+
+#include "ElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class MaxwellPlaneStrain;
+ class TestMaxwellPlaneStrain; // unit testing
+ } // materials
+} // pylith
+
+/// 2-D, isotropic, linear Maxwell viscoelastic material.
+class pylith::materials::MaxwellPlaneStrain : public ElasticMaterial
+{ // class MaxwellPlaneStrain
+ friend class TestMaxwellPlaneStrain; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor
+ MaxwellPlaneStrain(void);
+
+ /// Destructor
+ ~MaxwellPlaneStrain(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::MaxwellPlaneStrain::*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::MaxwellPlaneStrain::*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::MaxwellPlaneStrain::*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
+ MaxwellPlaneStrain(const MaxwellPlaneStrain& m);
+
+ /// Not implemented
+ const MaxwellPlaneStrain& operator=(const MaxwellPlaneStrain& 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 MaxwellPlaneStrain
+
+#include "MaxwellPlaneStrain.icc" // inline methods
+
+#endif // pylith_materials_maxwellplanestrain_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.icc 2009-04-09 04:54:56 UTC (rev 14644)
@@ -0,0 +1,112 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_materials_maxwellplanestrain_hh)
+#error "MaxwellPlaneStrain.icc can only be included from MaxwellPlaneStrain.hh"
+#endif
+
+#include <assert.h> // USES assert()
+#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+
+// Set current time step.
+inline
+void
+pylith::materials::MaxwellPlaneStrain::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::MaxwellPlaneStrain::useElasticBehavior(const bool flag) {
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcElasticConstsElastic;
+ _updatePropertiesFn =
+ &pylith::materials::MaxwellPlaneStrain::_updatePropertiesElastic;
+ } else {
+ _calcStressFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcStressViscoelastic;
+ _calcElasticConstsFn =
+ &pylith::materials::MaxwellPlaneStrain::_calcElasticConstsViscoelastic;
+ _updatePropertiesFn =
+ &pylith::materials::MaxwellPlaneStrain::_updatePropertiesViscoelastic;
+ } // if/else
+} // useElasticBehavior
+
+// Get flag indicating whether material implements an empty
+inline
+bool
+pylith::materials::MaxwellPlaneStrain::usesUpdateProperties(void) const {
+ return true;
+} // usesUpdateProperties
+
+// Compute stress tensor from parameters.
+inline
+void
+pylith::materials::MaxwellPlaneStrain::_calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* 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::MaxwellPlaneStrain::_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::MaxwellPlaneStrain::_updateProperties(double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialState,
+ const int initialStateSize) {
+ assert(0 != _updatePropertiesFn);
+ CALL_MEMBER_FN(*this, _updatePropertiesFn)(parameters, numParams,
+ totalStrain, strainSize,
+ initialState, initialStateSize);
+} // _updateProperties
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-04-09 03:05:05 UTC (rev 14643)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-04-09 04:54:56 UTC (rev 14644)
@@ -37,7 +37,7 @@
const int numElasticConsts = 21;
/// Number of physical properties.
- const int numProperties = 7;
+ const int numProperties = 8;
/// Physical properties.
const Material::PropMetaData properties[] = {
@@ -46,8 +46,9 @@
{ "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 },
+ { "total_strain", tensorSize, OTHER_FIELD },
+ { "viscous_strain", tensorSize, OTHER_FIELD },
+ { "stress_t", tensorSize, OTHER_FIELD },
};
/// Indices (order) of properties.
const int pidDensity = 0;
@@ -57,6 +58,7 @@
const int pidPowerLawExp = pidViscosityCoeff + 1;
const int pidStrainT = pidPowerLawExp + 1;
const int pidVisStrain = pidStrainT + tensorSize;
+ const int pidStressT = pidVisStrain + tensorSize;
/// Values expected in spatial database
const int numDBValues = 5;
@@ -96,6 +98,7 @@
{ // constructor
_dimension = 3;
_visStrain.resize(_PowerLaw3D::tensorSize);
+ _stressT.resize(_PowerLaw3D::tensorSize);
} // constructor
// ----------------------------------------------------------------------
@@ -171,8 +174,7 @@
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 viscosityCoeffScale = (pressureScale^(1.0/powerLawExp))/timeScale;
const double powerLawExpScale = 1.0;
values[_PowerLaw3D::pidDensity] =
_normalizer->nondimensionalize(values[_PowerLaw3D::pidDensity],
@@ -208,8 +210,7 @@
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 viscosityCoeffScale = (pressureScale^(1.0/powerLawExp))/timeScale;
const double powerLawExpScale = 1.0;
values[_PowerLaw3D::pidDensity] =
_normalizer->dimensionalize(values[_PowerLaw3D::pidDensity],
@@ -278,7 +279,8 @@
// ----------------------------------------------------------------------
// Compute viscous strain for current time step.
-// material.
+//********** May not need this. Check formulation. I probably still need
+// updateState. *************
void
pylith::materials::PowerLaw3D::_computeStateVars(
const double* properties,
@@ -380,19 +382,109 @@
} // _calcStressElastic
// ----------------------------------------------------------------------
+// Effective stress function that computes effective stress function only
+// (no derivative).
+double
+pylith::materials::PowerLaw3D::_effStressFunc(double x,
+ void *params)
+{ // _effStressFunc
+ struct effStressParams *p = (struct effStressParams *) params;
+ double ae = p->ae;
+ double b = p->b;
+ double c = p->c;
+ double d = p->d;
+ double alfa = p->alfa;
+ double dt = p->dt;
+ double effStressT = p->effStressT;
+ double powerLawExp = p->powerLawExp;
+ double viscosityCoeff = p->viscosityCoeff;
+ double factor1 = 1.0-alfa;
+ double effStressTau = factor1 * effStressT + alfa * x;
+ double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
+ (powerLawExp - 1.0))/viscosityCoeff;
+ double a = ae + alfa * dt * gammaTau;
+ double y = a * a * x * x - b +
+ c * gammaTau - d * d * gammaTau * gammaTau;
+ PetscLogFlops(21);
+ return y;
+} // _effStressFunc
+
+// ----------------------------------------------------------------------
+// Effective stress function that computes effective stress function
+// derivative only (no function value).
+double
+pylith::materials::PowerLaw3D::_effStressDFunc(double x,
+ void *params)
+{ // _effStressDFunc
+ struct effStressParams *p = (struct effStressParams *) params;
+ double ae = p->ae;
+ double c = p->c;
+ double d = p->d;
+ double alfa = p->alfa;
+ double dt = p->dt;
+ double effStressT = p->effStressT;
+ double powerLawExp = p->powerLawExp;
+ double viscosityCoeff = p->viscosityCoeff;
+ double factor1 = 1.0-alfa;
+ double effStressTau = factor1 * effStressT + alfa * x;
+ double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
+ (powerLawExp - 1.0))/viscosityCoeff;
+ double a = ae + alfa * dt * gammaTau;
+ double dGammaTau = 0.5 * alfa * (powerLawExp - 1.0) *
+ ((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
+ (viscosityCoeff * viscosityCoeff);
+ double dy = 2.0 * a * a * x + dGammaTau *
+ (2.0 * a * alfa * dt * x * x + c - 2.0 * d * d * gammaTau);
+ PetscLogFlops(36);
+ return dy;
+} // _effStressDFunc
+
+// ----------------------------------------------------------------------
+// Effective stress function that computes effective stress function
+// and derivative.
+void
+pylith::materials::PowerLaw3D::_effStressFuncDFunc(double x,
+ void *params,
+ double *y,
+ double *dy)
+{ // _effStressFunc
+ struct effStressParams *p = (struct effStressParams *) params;
+ double ae = p->ae;
+ double b = p->b;
+ double c = p->c;
+ double d = p->d;
+ double alfa = p->alfa;
+ double dt = p->dt;
+ double effStressT = p->effStressT;
+ double powerLawExp = p->powerLawExp;
+ double viscosityCoeff = p->viscosityCoeff;
+ double factor1 = 1.0-alfa;
+ double effStressTau = factor1 * effStressT + alfa * x;
+ double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
+ (powerLawExp - 1.0))/viscosityCoeff;
+ double dGammaTau = 0.5 * alfa * (powerLawExp - 1.0) *
+ ((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
+ (viscosityCoeff * viscosityCoeff);
+ double a = ae + alfa * dt * gammaTau;
+ double *y = a * a * x * x - b + c * gammaTau - d * d * gammaTau * gammaTau;
+ double *dy = 2.0 * a * a * x + dGammaTau *
+ (2.0 * a * alfa * dt * x * x + c - 2.0 * d * d * gammaTau);
+ PetscLogFlops(46);
+} // _effStressFunc
+// ----------------------------------------------------------------------
// 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)
+ 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);
More information about the CIG-COMMITS
mailing list