[cig-commits] r18489 - in short/3D/PyLith/trunk: libsrc/pylith libsrc/pylith/materials modulesrc/materials pylith pylith/materials unittests/libtests/materials unittests/libtests/materials/data
willic3 at geodynamics.org
willic3 at geodynamics.org
Sun May 29 16:59:03 PDT 2011
Author: willic3
Date: 2011-05-29 16:59:03 -0700 (Sun, 29 May 2011)
New Revision: 18489
Added:
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.hh
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.icc
short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellPlaneStrain.i
short/3D/PyLith/trunk/pylith/materials/GenMaxwellPlaneStrain.py
short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.hh
Modified:
short/3D/PyLith/trunk/libsrc/pylith/Makefile.am
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/pylith/materials/Makefile.am
short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh
short/3D/PyLith/trunk/modulesrc/materials/Makefile.am
short/3D/PyLith/trunk/modulesrc/materials/materials.i
short/3D/PyLith/trunk/pylith/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh
Log:
Added GenMaxwellPlaneStrain and made some changes to MaxwellPlaneStrain.
Primarily, added initial stress in the out-of-plane direction as a state
variable.
Modified: short/3D/PyLith/trunk/libsrc/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/Makefile.am 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/libsrc/pylith/Makefile.am 2011-05-29 23:59:03 UTC (rev 18489)
@@ -115,6 +115,7 @@
materials/ElasticIsotropic3D.cc \
materials/ViscoelasticMaxwell.cc \
materials/GenMaxwellIsotropic3D.cc \
+ materials/GenMaxwellPlaneStrain.cc \
materials/GenMaxwellQpQsIsotropic3D.cc \
materials/MaxwellIsotropic3D.cc \
materials/MaxwellPlaneStrain.cc \
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.hh 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.hh 2011-05-29 23:59:03 UTC (rev 18489)
@@ -44,7 +44,7 @@
* the value is less than one, the remainder of the total shear
* modulus is associated with a spring in parallel with the Maxwell
* models. The physical properties are stored internally using
- * density, lambdaTot, muTot, which are directly related to the
+ * density, lambdaEff, muEff, which are directly related to the
* elasticity constants used in the finite-element integration. The
* viscosity for each model is stored using Maxwell Time
* (viscosity/mu), and the shear ratio is also stored for each Maxwell
@@ -424,12 +424,16 @@
/** Update state variables after solve as a viscoelastic material.
*
+ * @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 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 values.
+ * @param initialStrainSize Size of initial strain array.
*/
void _updateStateVarsViscoelastic(double* const stateVars,
const int numStateVars,
Added: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,951 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "GenMaxwellPlaneStrain.hh" // implementation of object methods
+
+#include "ViscoelasticMaxwell.hh" // USES computeVisStrain
+#include "Metadata.hh" // USES Metadata
+
+#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+
+#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 _GenMaxwellPlaneStrain{
+
+ /// Number of Maxwell models in parallel.
+ const int numMaxwellModels = 3;
+
+ // Dimension of material.
+ const int dimension = 2;
+
+ /// Number of entries in stress/strain tensors.
+ const int tensorSize = 3;
+
+ /// Number of entries in derivative of elasticity matrix.
+ const int numElasticConsts = 9;
+
+ /// Number of physical properties.
+ const int numProperties = 5;
+
+ /// Physical properties.
+ const Metadata::ParamDescription properties[numProperties] = {
+ { "density", 1, pylith::topology::FieldBase::SCALAR },
+ { "mu", 1, pylith::topology::FieldBase::SCALAR },
+ { "lambda", 1, pylith::topology::FieldBase::SCALAR },
+ { "shear_ratio", numMaxwellModels, pylith::topology::FieldBase::OTHER },
+ { "maxwell_time", numMaxwellModels, pylith::topology::FieldBase::OTHER },
+ };
+ // Values expected in properties spatial database. :KLUDGE: Not
+ // generalized over number of models.
+ const int numDBProperties = 3 + 2 * numMaxwellModels;
+ const char* dbProperties[numDBProperties] = {
+ "density", "vs", "vp",
+ "shear-ratio-1",
+ "shear-ratio-2",
+ "shear-ratio-3",
+ "viscosity-1",
+ "viscosity-2",
+ "viscosity-3",
+ };
+
+ /// Number of state variables.
+ const int numStateVars = 2 + numMaxwellModels;
+
+ /// State variables. :KLUDGE: Not generalized over number of models.
+ const Metadata::ParamDescription stateVars[numStateVars] = {
+ { "stress_zz_initial", 1, pylith::topology::FieldBase::SCALAR },
+ { "total_strain", tensorSize, pylith::topology::FieldBase::TENSOR },
+ { "viscous_strain_1", 4, pylith::topology::FieldBase::OTHER },
+ { "viscous_strain_2", 4, pylith::topology::FieldBase::OTHER },
+ { "viscous_strain_3", 4, pylith::topology::FieldBase::OTHER },
+ };
+
+ // Values expected in state variables spatial database
+ const int numDBStateVars = 1 + tensorSize + numMaxwellModels * 4;
+ const char* dbStateVars[numDBStateVars] = {"stress-zz-initial",
+ "total-strain-xx",
+ "total-strain-yy",
+ "total-strain-xy",
+ "viscous-strain-1-xx",
+ "viscous-strain-1-yy",
+ "viscous-strain-1-zz",
+ "viscous-strain-1-xy",
+ "viscous-strain-2-xx",
+ "viscous-strain-2-yy",
+ "viscous-strain-2-zz",
+ "viscous-strain-2-xy",
+ "viscous-strain-3-xx",
+ "viscous-strain-3-yy",
+ "viscous-strain-3-zz",
+ "viscous-strain-3-xy",
+ };
+
+ } // _GenMaxwellPlaneStrain
+ } // materials
+} // pylith
+
+// Indices of physical properties
+const int pylith::materials::GenMaxwellPlaneStrain::p_density = 0;
+
+const int pylith::materials::GenMaxwellPlaneStrain::p_muEff =
+ pylith::materials::GenMaxwellPlaneStrain::p_density + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::p_lambdaEff =
+ pylith::materials::GenMaxwellPlaneStrain::p_muEff + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::p_shearRatio =
+ pylith::materials::GenMaxwellPlaneStrain::p_lambdaEff + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::p_maxwellTime =
+ pylith::materials::GenMaxwellPlaneStrain::p_shearRatio +
+ pylith::materials::_GenMaxwellPlaneStrain::numMaxwellModels;
+
+// Indices of database values (order must match dbProperties)
+const int pylith::materials::GenMaxwellPlaneStrain::db_density = 0;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_vs =
+ pylith::materials::GenMaxwellPlaneStrain::db_density + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_vp =
+ pylith::materials::GenMaxwellPlaneStrain::db_vs + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_shearRatio =
+ pylith::materials::GenMaxwellPlaneStrain::db_vp + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_viscosity =
+ pylith::materials::GenMaxwellPlaneStrain::db_shearRatio +
+ pylith::materials::_GenMaxwellPlaneStrain::numMaxwellModels;
+
+// Indices of state variables
+const int pylith::materials::GenMaxwellPlaneStrain::s_stressZZInitial = 0;
+
+const int pylith::materials::GenMaxwellPlaneStrain::s_totalStrain =
+ pylith::materials::GenMaxwellPlaneStrain::s_stressZZInitial + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::s_viscousStrain1 =
+ pylith::materials::GenMaxwellPlaneStrain::s_totalStrain +
+ pylith::materials::_GenMaxwellPlaneStrain::tensorSize;
+
+const int pylith::materials::GenMaxwellPlaneStrain::s_viscousStrain2 =
+ pylith::materials::GenMaxwellPlaneStrain::s_viscousStrain1 + 4;
+
+const int pylith::materials::GenMaxwellPlaneStrain::s_viscousStrain3 =
+ pylith::materials::GenMaxwellPlaneStrain::s_viscousStrain2 + 4;
+
+// Indices of database values (order must match dbStateVars)
+const int pylith::materials::GenMaxwellPlaneStrain::db_stressZZInitial = 0;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_totalStrain =
+ pylith::materials::GenMaxwellPlaneStrain::db_stressZZInitial + 1;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_viscousStrain1 =
+ pylith::materials::GenMaxwellPlaneStrain::db_totalStrain +
+ pylith::materials::_GenMaxwellPlaneStrain::tensorSize;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_viscousStrain2 =
+ pylith::materials::GenMaxwellPlaneStrain::db_viscousStrain1 + 4;
+
+const int pylith::materials::GenMaxwellPlaneStrain::db_viscousStrain3 =
+ pylith::materials::GenMaxwellPlaneStrain::db_viscousStrain2 + 4;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::materials::GenMaxwellPlaneStrain::GenMaxwellPlaneStrain(void) :
+ ElasticMaterial(_GenMaxwellPlaneStrain::dimension,
+ _GenMaxwellPlaneStrain::tensorSize,
+ _GenMaxwellPlaneStrain::numElasticConsts,
+ Metadata(_GenMaxwellPlaneStrain::properties,
+ _GenMaxwellPlaneStrain::numProperties,
+ _GenMaxwellPlaneStrain::dbProperties,
+ _GenMaxwellPlaneStrain::numDBProperties,
+ _GenMaxwellPlaneStrain::stateVars,
+ _GenMaxwellPlaneStrain::numStateVars,
+ _GenMaxwellPlaneStrain::dbStateVars,
+ _GenMaxwellPlaneStrain::numDBStateVars)),
+ _calcElasticConstsFn(0),
+ _calcStressFn(0),
+ _updateStateVarsFn(0)
+{ // constructor
+ useElasticBehavior(true);
+ _viscousStrain.resize(_GenMaxwellPlaneStrain::numMaxwellModels * 4);
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::materials::GenMaxwellPlaneStrain::~GenMaxwellPlaneStrain(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set whether elastic or inelastic constitutive relations are used.
+void
+pylith::materials::GenMaxwellPlaneStrain::useElasticBehavior(const bool flag)
+{ // useElasticBehavior
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::GenMaxwellPlaneStrain::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::GenMaxwellPlaneStrain::_calcElasticConstsElastic;
+ _updateStateVarsFn =
+ &pylith::materials::GenMaxwellPlaneStrain::_updateStateVarsElastic;
+
+ } else {
+ _calcStressFn =
+ &pylith::materials::GenMaxwellPlaneStrain::_calcStressViscoelastic;
+ _calcElasticConstsFn =
+ &pylith::materials::GenMaxwellPlaneStrain::_calcElasticConstsViscoelastic;
+ _updateStateVarsFn =
+ &pylith::materials::GenMaxwellPlaneStrain::_updateStateVarsViscoelastic;
+ } // if/else
+} // useElasticBehavior
+
+// ----------------------------------------------------------------------
+// Compute parameters from values in spatial database.
+void
+pylith::materials::GenMaxwellPlaneStrain::_dbToProperties(
+ double* const propValues,
+ const double_array& dbValues)
+{ // _dbToProperties
+ assert(0 != propValues);
+ const int numDBValues = dbValues.size();
+ assert(_GenMaxwellPlaneStrain::numDBProperties == numDBValues);
+
+ const double density = dbValues[db_density];
+ const double vs = dbValues[db_vs];
+ const double vp = dbValues[db_vp];
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
+ std::ostringstream msg;
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
+ << "density: " << density << "\n"
+ << "vp: " << vp << "\n"
+ << "vs: " << vs << "\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[p_density] = density;
+ propValues[p_muEff] = mu;
+ propValues[p_lambdaEff] = lambda;
+
+ double visFrac = 0.0;
+ for (int imodel = 0; imodel < numMaxwellModels; ++imodel)
+ visFrac += dbValues[db_shearRatio + imodel];
+ if (visFrac > 1.0) {
+ std::ostringstream msg;
+ msg << "Shear modulus ratios sum to a value greater than 1.0 for\n"
+ << "Generalized Maxwell model.\n"
+ << "Ratio 1: " << propValues[db_shearRatio ] << "\n"
+ << "Ratio 2: " << propValues[db_shearRatio+1] << "\n"
+ << "Ratio 3: " << propValues[db_shearRatio+2] << "\n"
+ << "Total: " << visFrac << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ // Loop over number of Maxwell models.
+ for (int imodel = 0; imodel < numMaxwellModels; ++imodel) {
+ double muRatio = dbValues[db_shearRatio + imodel];
+ double viscosity = dbValues[db_viscosity + imodel];
+ double muFac = muRatio * mu;
+ double maxwellTime = pylith::PYLITH_MAXDOUBLE;
+ if (muFac > 0.0)
+ maxwellTime = viscosity / muFac;
+ if (muRatio < 0.0 || viscosity < 0.0 || muFac < 0.0 || maxwellTime < 0.0) {
+ std::ostringstream msg;
+ msg << "Found negative value(s) for physical properties.\n"
+ << "muRatio: " << muRatio << "\n"
+ << "viscosity: " << viscosity << "\n"
+ << "muFac: " << muFac << "\n"
+ << "maxwellTime: " << maxwellTime << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+ propValues[p_shearRatio + imodel] = muRatio;
+ propValues[p_maxwellTime + imodel] = maxwellTime;
+ } // for
+
+ PetscLogFlops(6 + 2 * numMaxwellModels);
+} // _dbToProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::GenMaxwellPlaneStrain::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _numPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ values[p_density] =
+ _normalizer->nondimensionalize(values[p_density], densityScale);
+ values[p_muEff] =
+ _normalizer->nondimensionalize(values[p_muEff], pressureScale);
+ values[p_lambdaEff] =
+ _normalizer->nondimensionalize(values[p_lambdaEff], pressureScale);
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+ _normalizer->nondimensionalize(&values[p_maxwellTime],
+ numMaxwellModels, timeScale);
+
+ PetscLogFlops(3 + 1 * numMaxwellModels);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::GenMaxwellPlaneStrain::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _numPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ values[p_density] =
+ _normalizer->dimensionalize(values[p_density], densityScale);
+ values[p_muEff] =
+ _normalizer->dimensionalize(values[p_muEff], pressureScale);
+ values[p_lambdaEff] =
+ _normalizer->dimensionalize(values[p_lambdaEff], pressureScale);
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+ _normalizer->dimensionalize(&values[p_maxwellTime],
+ numMaxwellModels, timeScale);
+
+ PetscLogFlops(3 + 1 * numMaxwellModels);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Compute initial state variables from values in spatial database.
+void
+pylith::materials::GenMaxwellPlaneStrain::_dbToStateVars(
+ double* const stateValues,
+ const double_array& dbValues)
+{ // _dbToStateVars
+ assert(0 != stateValues);
+ const int numDBValues = dbValues.size();
+ assert(_GenMaxwellPlaneStrain::numDBStateVars == numDBValues);
+
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+ const int totalSize = 1 + _tensorSize + 4 * numMaxwellModels;
+ assert(totalSize == _numVarsQuadPt);
+ assert(totalSize == numDBValues);
+ for (int i=0; i < totalSize; ++i)
+ stateValues[i] = dbValues[i];
+
+ PetscLogFlops(0);
+} // _dbToStateVars
+
+// ----------------------------------------------------------------------
+// Nondimensionalize state variables.
+void
+pylith::materials::GenMaxwellPlaneStrain::_nondimStateVars(
+ double* const values,
+ const int nvalues) const
+{ // _nondimStateVars
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _numVarsQuadPt);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(&values[s_stressZZInitial], 1, pressureScale);
+
+ PetscLogFlops(1);
+} // _nondimStateVars
+
+// ----------------------------------------------------------------------
+// Dimensionalize state variables.
+void
+pylith::materials::GenMaxwellPlaneStrain::_dimStateVars(
+ double* const values,
+ const int nvalues) const
+{ // _dimStateVars
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _numVarsQuadPt);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(&values[s_stressZZInitial], 1, pressureScale);
+
+ PetscLogFlops(1);
+} // _dimStateVars
+
+// ----------------------------------------------------------------------
+// Compute density at location from properties.
+void
+pylith::materials::GenMaxwellPlaneStrain::_calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars)
+{ // _calcDensity
+ assert(0 != density);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+
+ density[0] = properties[p_density];
+} // _calcDensity
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties as an elastic
+// material.
+void
+pylith::materials::GenMaxwellPlaneStrain::_calcStressElastic(
+ 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)
+{ // _calcStressElastic
+ assert(0 != stress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == stressSize);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != totalStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
+ const double mu2 = 2.0 * mu;
+
+ // :TODO: Need to consider initial state variables????
+ const double e11 = totalStrain[0] - initialStrain[0];
+ const double e22 = totalStrain[1] - initialStrain[1];
+ const double e12 = totalStrain[2] - initialStrain[2];
+
+ const double s12 = lambda * (e11 + e22);
+
+ stress[0] = s12 + mu2 * e11 + initialStress[0];
+ stress[1] = s12 + mu2 * e22 + initialStress[1];
+ stress[2] = mu2 * e12 + initialStress[2];
+
+ PetscLogFlops(14);
+} // _calcStressElastic
+
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties as a viscoelastic
+// material.
+void
+pylith::materials::GenMaxwellPlaneStrain::_calcStressViscoelastic(
+ 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)
+{ // _calcStressViscoelastic
+ assert(0 != stress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == stressSize);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != totalStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+ const int tensorSize = _GenMaxwellPlaneStrain::tensorSize;
+
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
+ const double muRatio[numMaxwellModels] = {
+ properties[p_shearRatio ],
+ properties[p_shearRatio+1],
+ properties[p_shearRatio+2]
+ };
+ const double stressZZInitial = stateVars[s_stressZZInitial];
+
+ const double mu2 = 2.0 * mu;
+ const double bulkModulus = lambda + mu2/3.0;
+
+ // Initial stress and strain values
+ const double meanStrainInitial = (initialStrain[0] + initialStrain[1]) / 3.0;
+ const double meanStressInitial = (initialStress[0] + initialStress[1] +
+ stressZZInitial) / 3.0;
+ const double devStrainInitial[] = {initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ initialStrain[2]};
+ const double devStressInitial[] = {initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ initialStress[2]};
+
+ // Mean stress and strain for t + dt
+ const double meanStrainTpdt = (totalStrain[0] + totalStrain[1]) / 3.0;
+ const double meanStressTpdt = 3.0 * bulkModulus *
+ (meanStrainTpdt - meanStrainInitial) + meanStressInitial;
+
+ const double diag[] = { 1.0, 1.0, 0.0 };
+
+ double visFrac = 0.0;
+ for (int imodel=0; imodel < numMaxwellModels; ++imodel)
+ visFrac += muRatio[imodel];
+ assert(visFrac <= 1.0);
+ const double elasFrac = 1.0 - visFrac;
+
+ PetscLogFlops(18 + numMaxwellModels);
+
+ // Get viscous strains
+ if (computeStateVars) {
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+ } else {
+ int index = 0;
+ for (int iComp=0; iComp < 4; ++iComp)
+ _viscousStrain[index++] = stateVars[s_viscousStrain1+iComp];
+ for (int iComp=0; iComp < 4; ++iComp)
+ _viscousStrain[index++] = stateVars[s_viscousStrain2+iComp];
+ for (int iComp=0; iComp < 4; ++iComp)
+ _viscousStrain[index++] = stateVars[s_viscousStrain3+iComp];
+ } // else
+
+ // Compute new stresses
+ double devStrainTpdt = 0.0;
+ double devStressTpdt = 0.0;
+ const int visIndex[] = {0, 1, 3};
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt -
+ devStrainInitial[iComp];
+ devStressTpdt = elasFrac * devStrainTpdt;
+ for (int model=0; model < numMaxwellModels; ++model) {
+ devStressTpdt += muRatio[model] *
+ _viscousStrain[4 * model + visIndex[iComp]];
+ } // for
+
+ devStressTpdt = mu2 * devStressTpdt + devStressInitial[iComp];
+ stress[iComp] = diag[iComp] * (meanStressTpdt + meanStressInitial) +
+ devStressTpdt;
+ } // for
+
+ PetscLogFlops((9 + 2 * numMaxwellModels) * tensorSize);
+} // _calcStressViscoelastic
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from properties.
+void
+pylith::materials::GenMaxwellPlaneStrain::_calcElasticConstsElastic(
+ 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)
+{ // _calcElasticConstsElastic
+ assert(0 != elasticConsts);
+ assert(_GenMaxwellPlaneStrain::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != totalStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
+
+ const double mu2 = 2.0 * mu;
+ const double lambda2mu = lambda + mu2;
+
+ elasticConsts[ 0] = lambda2mu; // C1111
+ elasticConsts[ 1] = lambda; // C1122
+ elasticConsts[ 2] = 0; // C1112
+ elasticConsts[ 3] = lambda; // C2211
+ elasticConsts[ 4] = lambda2mu; // C2222
+ elasticConsts[ 5] = 0; // C2212
+ elasticConsts[ 6] = 0; // C1211
+ elasticConsts[ 7] = 0; // C1222
+ elasticConsts[ 8] = mu2; // C1212
+
+ PetscLogFlops(2);
+} // _calcElasticConstsElastic
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from properties
+// as a viscoelastic material.
+void
+pylith::materials::GenMaxwellPlaneStrain::_calcElasticConstsViscoelastic(
+ 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)
+{ // _calcElasticConstsViscoelastic
+ assert(0 != elasticConsts);
+ assert(_GenMaxwellPlaneStrain::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != totalStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+
+ const double mu = properties[p_muEff];
+ const double lambda = properties[p_lambdaEff];
+ const double mu2 = 2.0 * mu;
+ const double bulkModulus = lambda + mu2 / 3.0;
+
+ // Compute viscous contribution.
+ double visFac = 0.0;
+ double visFrac = 0.0;
+ double shearRatio = 0.0;
+ for (int imodel = 0; imodel < numMaxwellModels; ++imodel) {
+ shearRatio = properties[p_shearRatio + imodel];
+ double maxwellTime = pylith::PYLITH_MAXDOUBLE;
+ visFrac += shearRatio;
+ if (shearRatio != 0.0) {
+ maxwellTime = properties[p_maxwellTime + imodel];
+ visFac +=
+ shearRatio * ViscoelasticMaxwell::viscousStrainParam(_dt, maxwellTime);
+ } // if
+ } // for
+ double elasFrac = 1.0 - visFrac;
+ double shearFac = elasFrac + visFac;
+
+ elasticConsts[ 0] = bulkModulus + 4.0 * mu / 3.0 * shearFac; // C1111
+ elasticConsts[ 1] = bulkModulus - 2.0 * mu / 3.0 * shearFac; // C1122
+ elasticConsts[ 2] = 0; // C1112
+ elasticConsts[ 3] = elasticConsts[1]; // C2211
+ elasticConsts[ 4] = elasticConsts[0]; // C2222
+ elasticConsts[ 5] = 0; // C2212
+ elasticConsts[ 6] = 0; // C1211
+ elasticConsts[ 7] = 0; // C1222
+ elasticConsts[ 8] = 2.0 * mu * shearFac; // C1212
+
+ PetscLogFlops(15 + 3 * numMaxwellModels);
+} // _calcElasticConstsViscoelastic
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::GenMaxwellPlaneStrain::_updateStateVarsElastic(
+ double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVarsElastic
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ const int tensorSize = _tensorSize;
+
+ const double strainTpdt[] = {totalStrain[0] - initialStrain[0],
+ totalStrain[1] - initialStrain[1],
+ 0.0,
+ totalStrain[2] - initialStrain[2]};
+ const double meanStrainTpdt = (strainTpdt[0] + strainTpdt[1])/3.0;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0};
+
+ // Update total strain
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_totalStrain+iComp] = totalStrain[iComp];
+
+ // Initialize all viscous strains to deviatoric elastic strains.
+ double devStrain = 0.0;
+ double shearRatio = 0.0;
+ for (int iComp=0; iComp < 4; ++iComp) {
+ devStrain = strainTpdt[iComp] - diag[iComp] * meanStrainTpdt;
+ // Maxwell model 1
+ stateVars[s_viscousStrain1+iComp] = devStrain;
+ // Maxwell model 2
+ stateVars[s_viscousStrain2+iComp] = devStrain;
+ // Maxwell model 3
+ stateVars[s_viscousStrain3+iComp] = devStrain;
+ } // for
+ PetscLogFlops(2 + 2 * 4);
+
+ _needNewJacobian = true;
+} // _updateStateVarsElastic
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::GenMaxwellPlaneStrain::_updateStateVarsViscoelastic(
+ double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVarsViscoelastic
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ _computeStateVars(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+
+ const int tensorSize = _tensorSize;
+
+ // Total strain
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_totalStrain+iComp] = totalStrain[iComp];
+
+ // Viscous strain
+ int index = 0;
+ for (int iComp=0; iComp < 4; ++iComp)
+ stateVars[s_viscousStrain1+iComp] = _viscousStrain[index++];
+ for (int iComp=0; iComp < 4; ++iComp)
+ stateVars[s_viscousStrain2+iComp] = _viscousStrain[index++];
+ for (int iComp=0; iComp < 4; ++iComp)
+ stateVars[s_viscousStrain3+iComp] = _viscousStrain[index++];
+
+ _needNewJacobian = false;
+} // _updateStateVarsViscoelastic
+
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+pylith::materials::GenMaxwellPlaneStrain::_stableTimeStepImplicit(
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const
+{ // _stableTimeStepImplicit
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+
+ double dtStable = pylith::PYLITH_MAXDOUBLE;
+
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+ for (int i=0; i < numMaxwellModels; ++i) {
+ const double maxwellTime = properties[p_maxwellTime+i];
+ const double dt = 0.2*maxwellTime;
+ if (dt < dtStable)
+ dtStable = dt;
+ } // for
+
+ return dtStable;
+} // _stableTimeStepImplicit
+
+#include <iostream>
+// ----------------------------------------------------------------------
+// Compute viscous strain for current time step.
+void
+pylith::materials::GenMaxwellPlaneStrain::_computeStateVars(
+ const double* stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _computeStateVars
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_GenMaxwellPlaneStrain::tensorSize == initialStrainSize);
+
+ const int tensorSize = _tensorSize;
+ const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
+
+ const double muRatio[numMaxwellModels] = {
+ properties[p_shearRatio ],
+ properties[p_shearRatio+1],
+ properties[p_shearRatio+2]
+ };
+ const double maxwellTime[numMaxwellModels] = {
+ properties[p_maxwellTime ],
+ properties[p_maxwellTime+1],
+ properties[p_maxwellTime+2]
+ };
+
+ const double strainTpdt[] = {totalStrain[0],
+ totalStrain[1],
+ 0.0,
+ totalStrain[2]};
+ const double strainT[] = {stateVars[s_totalStrain+0],
+ stateVars[s_totalStrain+1],
+ 0.0,
+ stateVars[s_totalStrain+2]};
+
+ const double meanStrainTpdt = (strainTpdt[0] + strainTpdt[1])/3.0;
+ const double meanStrainT = (strainT[0] + strainT[1]) / 3.0;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0 };
+
+ PetscLogFlops(4);
+
+ // Compute Prony series terms
+ double_array dq(numMaxwellModels);
+ dq = 0.0;
+ for (int i=0; i < numMaxwellModels; ++i)
+ if (muRatio[i] != 0.0)
+ dq[i] = ViscoelasticMaxwell::viscousStrainParam(_dt, maxwellTime[i]);
+
+ // Compute new viscous strains
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+ double deltaStrain = 0.0;
+
+ for (int iComp=0; iComp < 4; ++iComp) {
+ devStrainTpdt = strainTpdt[iComp] - diag[iComp] * meanStrainTpdt;
+ devStrainT = strainT[iComp] - diag[iComp] * meanStrainT;
+ deltaStrain = devStrainTpdt - devStrainT;
+
+ // Maxwell model 1
+ int imodel = 0;
+ if (0.0 != muRatio[imodel]) {
+ _viscousStrain[imodel * 4 + iComp] = exp(-_dt / maxwellTime[imodel]) *
+ stateVars[s_viscousStrain1 + iComp] + dq[imodel] * deltaStrain;
+ PetscLogFlops(6);
+ } // if
+
+ // Maxwell model 2
+ imodel = 1;
+ if (0.0 != muRatio[imodel]) {
+ _viscousStrain[imodel * 4 + iComp] = exp(-_dt / maxwellTime[imodel]) *
+ stateVars[s_viscousStrain2 + iComp] + dq[imodel] * deltaStrain;
+ PetscLogFlops(6);
+ } // if
+
+ // Maxwell model 3
+ imodel = 2;
+ if (0.0 != muRatio[imodel]) {
+ _viscousStrain[imodel * 4 + iComp] = exp(-_dt / maxwellTime[imodel]) *
+ stateVars[s_viscousStrain3 + iComp] + dq[imodel] * deltaStrain;
+ PetscLogFlops(6);
+ } // if
+
+ } // for
+
+ PetscLogFlops(numMaxwellModels * (4 * (5 + 18)));
+} // _computeStateVars
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.hh 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,540 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/GenMaxwellPlaneStrain.hh
+ *
+ * @brief 2-D, isotropic, generalized linear Maxwell viscoelastic material for
+ * plane strain.
+ */
+
+#if !defined(pylith_materials_genmaxwellplanestrain_hh)
+#define pylith_materials_genmaxwellplanestrain_hh
+
+// Include directives ---------------------------------------------------
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
+
+// GenMaxwellPlaneStrain ---------------------------------------------------
+/** @brief 2-D, isotropic, generalized linear Maxwell viscoelastic material for
+ * plane strain.
+ *
+ * This consists of several Maxwell models in parallel. At present,
+ * the number of models is fixed at 3, but this will be changed in the
+ * future. The physical properties are specified using density,
+ * shear-wave speed, and compressional-wave speed. A viscosity and a
+ * shear ratio are also given for each Maxwell model. The shear ratio
+ * specifies how much of the total shear modulus is associated with
+ * that model. The shear ratios must sum to a value less than one. If
+ * the value is less than one, the remainder of the total shear
+ * modulus is associated with a spring in parallel with the Maxwell
+ * models. The physical properties are stored internally using
+ * density, lambdaEff, muEff, which are directly related to the
+ * elasticity constants used in the finite-element integration. The
+ * viscosity for each model is stored using Maxwell Time
+ * (viscosity/mu), and the shear ratio is also stored for each Maxwell
+ * model.
+ */
+class pylith::materials::GenMaxwellPlaneStrain : public ElasticMaterial
+{ // class GenMaxwellPlaneStrain
+ friend class TestGenMaxwellPlaneStrain; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor
+ GenMaxwellPlaneStrain(void);
+
+ /// Destructor
+ ~GenMaxwellPlaneStrain(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);
+
+ // 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);
+
+ /** 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;
+
+ /** Compute initial state variables from values in spatial database.
+ *
+ * @param stateValues Array of state variable values.
+ * @param dbValues Array of database values.
+ */
+ void _dbToStateVars(double* const stateValues,
+ const double_array& dbValues);
+
+ /** Nondimensionalize state variables..
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _nondimStateVars(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize state variables.
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _dimStateVars(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.
+ * @param stateVars Number of state variables.
+ * @param numStateVars Number of state variables.
+ */
+ void _calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars);
+
+ /** 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 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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ 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 stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @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.
+ */
+ 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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** 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.
+ */
+ void _updateStateVars(double* const stateVars,
+ const int numStateVars,
+ const double* 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.
+ * @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 double* stateVars,
+ const int numStateVars) const;
+
+ // PRIVATE TYPEDEFS ///////////////////////////////////////////////////
+private :
+
+ /// Member prototype for _calcStress()
+ typedef void (pylith::materials::GenMaxwellPlaneStrain::*calcStress_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const bool);
+
+ /// Member prototype for _calcElasticConsts()
+ typedef void (pylith::materials::GenMaxwellPlaneStrain::*calcElasticConsts_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int);
+
+ /// Member prototype for _updateStateVars()
+ typedef void (pylith::materials::GenMaxwellPlaneStrain::*updateStateVars_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int);
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** 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 stateVars State variables at locations.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at locations.
+ * @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 computeStateVars Flag indicating to compute updated state vars.
+ */
+ void _calcStressElastic(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);
+
+ /** Compute stress tensor from properties as a 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 stateVars State variables at locations.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at locations.
+ * @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 computeStateVars Flag indicating to compute updated state vars.
+ */
+ void _calcStressViscoelastic(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);
+
+ /** 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 stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @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.
+ */
+ void _calcElasticConstsElastic(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);
+
+ /** 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 stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @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.
+ */
+ void _calcElasticConstsViscoelastic(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);
+
+ /** Update state variables after solve as an elastic material.
+ *
+ * @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.
+ */
+ void _updateStateVarsElastic(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** Update state variables after solve as a viscoelastic material.
+ *
+ * @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.
+ */
+ void _updateStateVarsViscoelastic(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** Compute viscous strains (state variables) for the current 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.
+ */
+ void _computeStateVars(const double* stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ /// Viscous strain array.
+ double_array _viscousStrain;
+
+ /// Method to use for _calcElasticConsts().
+ calcElasticConsts_fn_type _calcElasticConstsFn;
+
+ /// Method to use for _calcStress().
+ calcStress_fn_type _calcStressFn;
+
+ /// Method to use for _updateStateVars().
+ updateStateVars_fn_type _updateStateVarsFn;
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ static const int p_density;
+ static const int p_muEff;
+ static const int p_lambdaEff;
+ static const int p_shearRatio;
+ static const int p_maxwellTime;
+ static const int db_density;
+ static const int db_vs;
+ static const int db_vp;
+ static const int db_shearRatio;
+ static const int db_viscosity;
+
+ static const int s_stressZZInitial;
+ static const int s_totalStrain;
+ static const int s_viscousStrain1;
+ static const int s_viscousStrain2;
+ static const int s_viscousStrain3;
+ static const int db_stressZZInitial;
+ static const int db_totalStrain;
+ static const int db_viscousStrain1;
+ static const int db_viscousStrain2;
+ static const int db_viscousStrain3;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ GenMaxwellPlaneStrain(const GenMaxwellPlaneStrain&);
+
+ /// Not implemented
+ const GenMaxwellPlaneStrain& operator=(const GenMaxwellPlaneStrain&);
+
+}; // class GenMaxwellPlaneStrain
+
+#include "GenMaxwellPlaneStrain.icc" // inline methods
+
+#endif // pylith_materials_genmaxwellplanestrain_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.icc 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,109 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_materials_genmaxwellplanestrain_hh)
+#error "GenMaxwellPlaneStrain.icc can only be included from GenMaxwellPlaneStrain.hh"
+#endif
+
+#include <assert.h> // USES assert()
+#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+
+// Set current time step.
+inline
+void
+pylith::materials::GenMaxwellPlaneStrain::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
+
+// Compute stress tensor from parameters.
+inline
+void
+pylith::materials::GenMaxwellPlaneStrain::_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) {
+ assert(0 != _calcStressFn);
+ CALL_MEMBER_FN(*this, _calcStressFn)(stress, stressSize,
+ properties, numProperties,
+ stateVars, numStateVars,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize,
+ computeStateVars);
+} // _calcStress
+
+// Compute derivatives of elasticity matrix from parameters.
+inline
+void
+pylith::materials::GenMaxwellPlaneStrain::_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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize) {
+ assert(0 != _calcElasticConstsFn);
+ CALL_MEMBER_FN(*this, _calcElasticConstsFn)(elasticConsts, numElasticConsts,
+ properties, numProperties,
+ stateVars, numStateVars,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+} // _calcElasticConsts
+
+// Update state variables after solve.
+inline
+void
+pylith::materials::GenMaxwellPlaneStrain::_updateStateVars(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize) {
+ assert(0 != _updateStateVarsFn);
+ CALL_MEMBER_FN(*this, _updateStateVarsFn)(stateVars, numStateVars,
+ properties, numProperties,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+} // _updateStateVars
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Makefile.am 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Makefile.am 2011-05-29 23:59:03 UTC (rev 18489)
@@ -30,6 +30,8 @@
ElasticPlaneStress.hh \
GenMaxwellIsotropic3D.hh \
GenMaxwellIsotropic3D.icc \
+ GenMaxwellPlaneStrain.hh \
+ GenMaxwellPlaneStrain.icc \
GenMaxwellQpQsIsotropic3D.hh \
GenMaxwellQpQsIsotropic3D.icc \
MaxwellIsotropic3D.hh \
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh 2011-05-29 23:59:03 UTC (rev 18489)
@@ -44,6 +44,7 @@
class MaxwellIsotropic3D;
class MaxwellPlaneStrain;
class GenMaxwellIsotropic3D;
+ class GenMaxwellPlaneStrain;
class GenMaxwellQpQsIsotropic3D;
class PowerLaw3D;
class DruckerPrager3D;
Added: short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellPlaneStrain.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellPlaneStrain.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/materials/GenMaxwellPlaneStrain.i 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,209 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file modulesrc/materials/GenMaxwellPlaneStrain.i
+ *
+ * @brief Python interface to C++ GenMaxwellPlaneStrain object.
+ */
+
+namespace pylith {
+ namespace materials {
+
+ class pylith::materials::GenMaxwellPlaneStrain : public ElasticMaterial
+ { // class GenMaxwellPlaneStrain
+
+ // PUBLIC METHODS /////////////////////////////////////////////////
+ public :
+
+ /// Default constructor
+ GenMaxwellPlaneStrain(void);
+
+ /// Destructor
+ ~GenMaxwellPlaneStrain(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);
+
+ // 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 pylith::double_array& dbValues);
+
+ /** 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;
+
+ /** Compute initial state variables from values in spatial database.
+ *
+ * @param stateValues Array of state variable values.
+ * @param dbValues Array of database values.
+ */
+ void _dbToStateVars(double* const stateValues,
+ const pylith::double_array& dbValues);
+
+ /** Compute density from properties.
+ *
+ * @param density Array for density.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ */
+ void _calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars);
+
+ /** 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 initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ 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 stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @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* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** 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 values.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain values.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _updateStateVars(double* const stateVars,
+ const int numStateVars,
+ const double* 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.
+ * @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 double* stateVars,
+ const int numStateVars) const;
+
+ }; // class GenMaxwellPlaneStrain
+
+ } // materials
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2011-05-29 23:59:03 UTC (rev 18489)
@@ -35,6 +35,7 @@
MaxwellIsotropic3D.i \
MaxwellPlaneStrain.i \
GenMaxwellIsotropic3D.i \
+ GenMaxwellPlaneStrain.i \
GenMaxwellQpQsIsotropic3D.i \
PowerLaw3D.i \
DruckerPrager3D.i
Modified: short/3D/PyLith/trunk/modulesrc/materials/materials.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/materials.i 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/modulesrc/materials/materials.i 2011-05-29 23:59:03 UTC (rev 18489)
@@ -35,6 +35,7 @@
#include "pylith/materials/MaxwellIsotropic3D.hh"
#include "pylith/materials/MaxwellPlaneStrain.hh"
#include "pylith/materials/GenMaxwellIsotropic3D.hh"
+#include "pylith/materials/GenMaxwellPlaneStrain.hh"
#include "pylith/materials/GenMaxwellQpQsIsotropic3D.hh"
#include "pylith/materials/PowerLaw3D.hh"
#include "pylith/materials/DruckerPrager3D.hh"
@@ -77,6 +78,7 @@
%include "MaxwellIsotropic3D.i"
%include "MaxwellPlaneStrain.i"
%include "GenMaxwellIsotropic3D.i"
+%include "GenMaxwellPlaneStrain.i"
%include "GenMaxwellQpQsIsotropic3D.i"
%include "PowerLaw3D.i"
%include "DruckerPrager3D.i"
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2011-05-29 23:59:03 UTC (rev 18489)
@@ -76,6 +76,7 @@
materials/ElasticStrain1D.py \
materials/ElasticStress1D.py \
materials/GenMaxwellIsotropic3D.py \
+ materials/GenMaxwellPlaneStrain.py \
materials/GenMaxwellQpQsIsotropic3D.py \
materials/Homogeneous.py \
materials/Material.py \
Added: short/3D/PyLith/trunk/pylith/materials/GenMaxwellPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/GenMaxwellPlaneStrain.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/materials/GenMaxwellPlaneStrain.py 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/materials/GenMaxwellPlaneStrain.py
+##
+## @brief Python object implementing 3-D isotropic linear GenMaxwell
+## viscoelastic material.
+##
+## Factory: material.
+
+from ElasticMaterial import ElasticMaterial
+from materials import GenMaxwellPlaneStrain as ModuleGenMaxwellPlaneStrain
+
+# GenMaxwellPlaneStrain class
+class GenMaxwellPlaneStrain(ElasticMaterial, ModuleGenMaxwellPlaneStrain):
+ """
+ Python object implementing plane strain isotropic linear GenMaxwell
+ viscoelastic material.
+
+ Factory: material.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="genmaxwellplanestrain"):
+ """
+ Constructor.
+ """
+ ElasticMaterial.__init__(self, name)
+ self.availableFields = \
+ {'vertex': \
+ {'info': [],
+ 'data': []},
+ 'cell': \
+ {'info': ["mu", "lambda", "density",
+ "shear_ratio", "maxwell_time"],
+ 'data': ["stress_zz_initial",
+ "total_strain", "stress",
+ "viscous_strain_1",
+ "viscous_strain_2",
+ "viscous_strain_3",
+ ]}}
+ self._loggingPrefix = "MaMx2D "
+ return
+
+
+ def _createModuleObj(self):
+ """
+ Call constructor for module object for access to C++ object.
+ """
+ ModuleGenMaxwellPlaneStrain.__init__(self)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def material():
+ """
+ Factory associated with GenMaxwellPlaneStrain.
+ """
+ return GenMaxwellPlaneStrain()
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2011-05-29 23:59:03 UTC (rev 18489)
@@ -38,6 +38,7 @@
TestMaxwellIsotropic3D.cc \
TestMaxwellPlaneStrain.cc \
TestGenMaxwellIsotropic3D.cc \
+ TestGenMaxwellPlaneStrain.cc \
TestGenMaxwellQpQsIsotropic3D.cc \
TestPowerLaw3D.cc \
TestDruckerPrager3D.cc \
@@ -55,6 +56,7 @@
TestElasticPlaneStress.hh \
TestElasticIsotropic3D.hh \
TestGenMaxwellIsotropic3D.hh \
+ TestGenMaxwellPlaneStrain.hh \
TestGenMaxwellQpQsIsotropic3D.hh \
TestMaxwellIsotropic3D.hh \
TestMaxwellPlaneStrain.hh \
@@ -77,6 +79,8 @@
data/MaxwellPlaneStrainTimeDepData.cc \
data/GenMaxwellIsotropic3DElasticData.cc \
data/GenMaxwellIsotropic3DTimeDepData.cc \
+ data/GenMaxwellPlaneStrainElasticData.cc \
+ data/GenMaxwellPlaneStrainTimeDepData.cc \
data/GenMaxwellQpQsIsotropic3DElasticData.cc \
data/GenMaxwellQpQsIsotropic3DTimeDepData.cc \
data/PowerLaw3DElasticData.cc \
@@ -95,6 +99,8 @@
data/ElasticStress1DData.hh \
data/GenMaxwellIsotropic3DElasticData.hh \
data/GenMaxwellIsotropic3DTimeDepData.hh \
+ data/GenMaxwellPlaneStrainElasticData.hh \
+ data/GenMaxwellPlaneStrainTimeDepData.hh \
data/GenMaxwellQpQsIsotropic3DElasticData.hh \
data/GenMaxwellQpQsIsotropic3DTimeDepData.hh \
data/MaxwellIsotropic3DElasticData.hh \
Added: short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.cc 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,196 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestGenMaxwellPlaneStrain.hh" // Implementation of class methods
+
+#include "data/GenMaxwellPlaneStrainElasticData.hh" // USES GenMaxwellPlaneStrainElasticData
+#include "data/GenMaxwellPlaneStrainTimeDepData.hh" // USES GenMaxwellPlaneStrainTimeDepData
+
+#include "pylith/materials/GenMaxwellPlaneStrain.hh" // USES GenMaxwellPlaneStrain
+
+#include <cstring> // USES memcpy()
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestGenMaxwellPlaneStrain );
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::materials::TestGenMaxwellPlaneStrain::setUp(void)
+{ // setUp
+ _material = new GenMaxwellPlaneStrain();
+ _matElastic = new GenMaxwellPlaneStrain();
+
+ _data = new GenMaxwellPlaneStrainElasticData();
+ _dataElastic = new GenMaxwellPlaneStrainElasticData();
+ setupNormalizer();
+} // setUp
+
+// ----------------------------------------------------------------------
+// Test timeStep()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::testTimeStep(void)
+{ // testTimeStep
+ GenMaxwellPlaneStrain material;
+
+ CPPUNIT_ASSERT_EQUAL(false, material._needNewJacobian);
+
+ const double dt1 = 1.0;
+ material.timeStep(dt1);
+ CPPUNIT_ASSERT_EQUAL(dt1, material.Material::timeStep());
+ CPPUNIT_ASSERT_EQUAL(false, material.needNewJacobian());
+
+ const double dt2 = 2.0;
+ material.timeStep(dt2);
+ CPPUNIT_ASSERT_EQUAL(dt2, material.Material::timeStep());
+ CPPUNIT_ASSERT_EQUAL(true, material.needNewJacobian());
+} // testTimeStep
+
+// ----------------------------------------------------------------------
+// Test useElasticBehavior()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::testUseElasticBehavior(void)
+{ // testUseElasticBehavior
+ GenMaxwellPlaneStrain material;
+
+ material.useElasticBehavior(true);
+ // Some compilers/operating systems (cygwin) don't allow comparing
+ // pointers. Use first test to determine if we can compare pointers.
+ if (&pylith::materials::GenMaxwellPlaneStrain::_calcStressElastic ==
+ material._calcStressFn) {
+ CPPUNIT_ASSERT(&pylith::materials::GenMaxwellPlaneStrain::_calcStressElastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(&pylith::materials::GenMaxwellPlaneStrain::_calcElasticConstsElastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(&pylith::materials::GenMaxwellPlaneStrain::_updateStateVarsElastic ==
+ material._updateStateVarsFn);
+
+ material.useElasticBehavior(false);
+ CPPUNIT_ASSERT(&pylith::materials::GenMaxwellPlaneStrain::_calcStressViscoelastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(&pylith::materials::GenMaxwellPlaneStrain::_calcElasticConstsViscoelastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(&pylith::materials::GenMaxwellPlaneStrain::_updateStateVarsViscoelastic ==
+ material._updateStateVarsFn);
+ } // if
+} // testUseElasticBehavior
+
+// ----------------------------------------------------------------------
+// Test usesHasStateVars()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::testHasStateVars(void)
+{ // testHasStateVars
+ GenMaxwellPlaneStrain material;
+ CPPUNIT_ASSERT_EQUAL(true, material.hasStateVars());
+} // testHasStateVars
+
+// ----------------------------------------------------------------------
+// Test calcStressElastic()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::test_calcStressElastic(void)
+{ // testCalcStressElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_calcStress();
+} // testCalcStressElastic
+
+// ----------------------------------------------------------------------
+// Test calcElasticConstsElastic()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::test_calcElasticConstsElastic(void)
+{ // testElasticConstsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_calcElasticConsts();
+} // testElasticConstsElastic
+
+// ----------------------------------------------------------------------
+// Test updateStateVarsElastic()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::test_updateStateVarsElastic(void)
+{ // testUpdateStateVarsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_updateStateVars();
+} // testUpdateStateVarsElastic
+
+// ----------------------------------------------------------------------
+// Test calcStressTimeDep()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::test_calcStressTimeDep(void)
+{ // testCalcStressTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new GenMaxwellPlaneStrainTimeDepData();
+
+ double dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcStress();
+} // testCalcStressTimeDep
+
+// ----------------------------------------------------------------------
+// Test calcElasticConstsTimeDep()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::test_calcElasticConstsTimeDep(void)
+{ // testElasticConstsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new GenMaxwellPlaneStrainTimeDepData();
+
+ double dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcElasticConsts();
+} // testElasticConstsTimeDep
+
+// ----------------------------------------------------------------------
+// Test updateStateVarsTimeDep()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::test_updateStateVarsTimeDep(void)
+{ // testUpdateStateVarsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new GenMaxwellPlaneStrainTimeDepData();
+
+ double dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_updateStateVars();
+
+} // testUpdateStateVarsTimeDep
+
+// ----------------------------------------------------------------------
+// Test _stableTimeStepImplicit()
+void
+pylith::materials::TestGenMaxwellPlaneStrain::test_stableTimeStepImplicit(void)
+{ // test_stableTimeStepImplicit
+ CPPUNIT_ASSERT(0 != _matElastic);
+
+ delete _dataElastic; _dataElastic = new GenMaxwellPlaneStrainTimeDepData();
+
+ TestElasticMaterial::test_stableTimeStepImplicit();
+} // test_stableTimeStepImplicit
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellPlaneStrain.hh 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,114 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/materials/TestGenMaxwellPlaneStrain.hh
+ *
+ * @brief C++ TestGenMaxwellPlaneStrain object
+ *
+ * C++ unit testing for GenMaxwellPlaneStrain.
+ */
+
+#if !defined(pylith_materials_testgenmaxwellplanestrain_hh)
+#define pylith_materials_testgenmaxwellplanestrain_hh
+
+#include "TestElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class TestGenMaxwellPlaneStrain;
+ class GenMaxwellPlaneStrainElasticData;
+ class GenMaxwellPlaneStrainTimeDepData;
+ } // materials
+} // pylith
+
+/// C++ unit testing for GenMaxwellPlaneStrain
+class pylith::materials::TestGenMaxwellPlaneStrain : public TestElasticMaterial
+{ // class TestGenMaxwellPlaneStrain
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestGenMaxwellPlaneStrain );
+
+ CPPUNIT_TEST( testDimension );
+ CPPUNIT_TEST( testTensorSize );
+ CPPUNIT_TEST( testDBToProperties );
+ CPPUNIT_TEST( testNonDimProperties );
+ CPPUNIT_TEST( testDimProperties );
+ CPPUNIT_TEST( testDBToStateVars );
+ CPPUNIT_TEST( testNonDimStateVars );
+ CPPUNIT_TEST( testDimStateVars );
+ CPPUNIT_TEST( test_calcDensity );
+ CPPUNIT_TEST( test_stableTimeStepImplicit );
+
+ // Need to test Maxwell viscoelastic specific behavior.
+ CPPUNIT_TEST( testTimeStep );
+ CPPUNIT_TEST( testUseElasticBehavior );
+ CPPUNIT_TEST( testHasStateVars );
+
+ CPPUNIT_TEST( test_calcStressElastic );
+ CPPUNIT_TEST( test_calcStressTimeDep );
+ CPPUNIT_TEST( test_calcElasticConstsElastic );
+ CPPUNIT_TEST( test_calcElasticConstsTimeDep );
+ CPPUNIT_TEST( test_updateStateVarsElastic );
+ CPPUNIT_TEST( test_updateStateVarsTimeDep );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+ /// Test timeStep()
+ void testTimeStep(void);
+
+ /// Test useElasticBehavior()
+ void testUseElasticBehavior(void);
+
+ /// Test hasStateVars()
+ void testHasStateVars(void);
+
+ /// Test _calcStressElastic()
+ void test_calcStressElastic(void);
+
+ /// Test _calcElasticConstsElastic()
+ void test_calcElasticConstsElastic(void);
+
+ /// Test _updateStateVarsElastic()
+ void test_updateStateVarsElastic(void);
+
+ /// Test _calcStressTimeDep()
+ void test_calcStressTimeDep(void);
+
+ /// Test _calcElasticConstsTimeDep()
+ void test_calcElasticConstsTimeDep(void);
+
+ /// Test _updateStatevarsTimeDep()
+ void test_updateStateVarsTimeDep(void);
+
+ /// Test _stableTimeStepImplicit()
+ void test_stableTimeStepImplicit(void);
+
+}; // class TestGenMaxwellPlaneStrain
+
+#endif // pylith_materials_testgenmaxwellplanestrain_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,252 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/GenMaxwellPlaneStrainElastic.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ GenMaxwellPlaneStrain object with elastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+
+# ----------------------------------------------------------------------
+dimension = 2
+numElasticConsts = 9
+tensorSize = 3
+numMaxwellModels = 3
+
+# GenMaxwellPlaneStrainElastic class
+class GenMaxwellPlaneStrainElastic(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ GenMaxwellPlaneStrain object with elastic behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="genmaxwellplanestrainelastic"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ import pdb
+ pdb.set_trace()
+
+ numLocs = 2
+
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "shear-ratio-1", "shear-ratio-2", "shear-ratio-3",
+ "viscosity-1", "viscosity-2", "viscosity-3"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1, 1, 1, 1],
+ dtype=numpy.int32)
+
+ self.dbStateVarValues = ["stress-zz-initial",
+ "total-strain-xx",
+ "total-strain-yy",
+ "total-strain-xy",
+ "viscous-strain-1-xx",
+ "viscous-strain-1-yy",
+ "viscous-strain-1-zz",
+ "viscous-strain-1-xy",
+ "viscous-strain-2-xx",
+ "viscous-strain-2-yy",
+ "viscous-strain-2-zz",
+ "viscous-strain-2-xy",
+ "viscous-strain-3-xx",
+ "viscous-strain-3-yy",
+ "viscous-strain-3-zz",
+ "viscous-strain-3-xy"
+ ]
+ self.numStateVarValues = numpy.array([1, tensorSize, 4, 4, 4],
+ dtype=numpy.int32)
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ shearRatioA = [0.5, 0.1, 0.2]
+ viscosityA = [1.0e18, 1.0e17, 1.0e19]
+ strainA = [1.1e-4, 1.2e-4, 1.4e-4]
+ initialStressA = [2.1e4, 2.2e4, 2.4e4]
+ initialStrainA = [3.1e-5, 3.2e-5, 3.4e-5]
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+ stressInitialZZA = numpy.array([1.5e4], dtype=numpy.float64)
+
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ shearRatioB = [0.2, 0.2, 0.2]
+ viscosityB = [1.0e18, 1.0e19, 1.0e20]
+ strainB = [4.1e-4, 4.2e-4, 4.4e-4]
+ initialStressB = [5.1e4, 5.2e4, 5.4e4]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.4e-5]
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ stressInitialZZB = numpy.array([4.5e4], dtype=numpy.float64)
+
+ maxwellTimeA = [1.0e30, 1.0e30, 1.0e30]
+ maxwellTimeB = [1.0e30, 1.0e30, 1.0e30]
+ for i in xrange(numMaxwellModels):
+ if shearRatioA[i] != 0.0:
+ maxwellTimeA[i] = viscosityA[i]/(muA*shearRatioA[i])
+ if shearRatioB[i] != 0.0:
+ maxwellTimeB[i] = viscosityB[i]/(muB*shearRatioB[i])
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.timeScale = 1.0
+ self.densityScale = 1.0e+3
+
+ propA = [densityA, vsA, vpA] + shearRatioA + viscosityA
+ propB = [densityB, vsB, vpB] + shearRatioB + viscosityB
+ self.dbProperties = numpy.array([propA, propB], dtype=numpy.float64)
+ propA = [densityA, muA, lambdaA] + shearRatioA + maxwellTimeA
+ propB = [densityB, muB, lambdaB] + shearRatioB + maxwellTimeB
+ self.properties = numpy.array([propA, propB], dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ # At present, only the first (stressInitialZZ) is being used.
+ self.dbStateVars = numpy.zeros( (numLocs,
+ 1 + tensorSize + 4 * numMaxwellModels),
+ dtype=numpy.float64)
+ self.dbStateVars[0, 0] = stressInitialZZA
+ self.dbStateVars[1, 0] = stressInitialZZB
+
+ self.stateVars = numpy.zeros( (numLocs,
+ 1 + tensorSize + 4 * numMaxwellModels),
+ dtype=numpy.float64)
+ self.stateVars[0, 0] = stressInitialZZA
+ self.stateVars[1, 0] = stressInitialZZB
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ time0 = self.timeScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0,
+ shearRatioA[0], shearRatioA[1], shearRatioA[2],
+ maxwellTimeA[0]/time0, maxwellTimeA[1]/time0,
+ maxwellTimeA[2]/time0],
+ [densityB/density0, muB/mu0, lambdaB/mu0,
+ shearRatioB[0], shearRatioB[1], shearRatioB[2],
+ maxwellTimeB[0]/time0, maxwellTimeB[1]/time0,
+ maxwellTimeB[2]/time0] ],
+ dtype=numpy.float64)
+
+ stressInitialZZANondim = stressInitialZZA/mu0
+ stressInitialZZBNondim = stressInitialZZB/mu0
+
+ self.stateVarsNondim = numpy.zeros( (numLocs,
+ 1 + tensorSize + 4 * numMaxwellModels),
+ dtype=numpy.float64)
+
+ self.stateVarsNondim[0, 0] = stressInitialZZANondim
+ self.stateVarsNondim[1, 0] = stressInitialZZBNondim
+
+ self.initialStress = numpy.array([initialStressA,
+ initialStressB],
+ dtype=numpy.float64)
+ self.initialStrain = numpy.array([initialStrainA,
+ initialStrainB],
+ dtype=numpy.float64)
+
+ self.density = numpy.array([densityA,
+ densityB],
+ dtype=numpy.float64)
+
+ self.strain = numpy.array([strainA,
+ strainB],
+ dtype=numpy.float64)
+
+ self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts), \
+ dtype=numpy.float64)
+ self.stateVarsUpdated = numpy.zeros((self.numLocs,
+ 1 + tensorSize + 4 * numMaxwellModels),
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA,
+ initialStressA, initialStrainA,
+ self.stateVars[0,:])
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB,
+ initialStressB, initialStrainB,
+ self.stateVars[1,:])
+ self.dtStableImplicit = 0.2*min(min(maxwellTimeA), min(maxwellTimeB))
+
+ return
+
+
+ def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV,
+ stateVars):
+ """
+ Compute stress and derivative of elasticity matrix.
+ """
+ C1111 = lambdaV + 2.0 * muV
+ C1122 = lambdaV
+ C1112 = 0.0
+ C2211 = lambdaV
+ C2222 = lambdaV + 2.0 * muV
+ C2212 = 0.0
+ C1211 = 0.0
+ C1222 = 0.0
+ C1212 = 2.0 * muV
+ elasticConsts = numpy.array([C1111, C1122, C1112,
+ C2211, C2222, C2212,
+ C1211, C1222, C1212],
+ dtype=numpy.float64)
+
+ initialStress = numpy.reshape(initialStressV, (tensorSize,1))
+ initialStrain = numpy.reshape(initialStrainV, (tensorSize,1))
+ stressZZInitial = numpy.array([stateVars[0]], dtype=numpy.float64)
+ strain = numpy.reshape(strainV, (tensorSize,1)) - initialStrain
+
+ elastic = numpy.array([ [C1111, C1122, C1112],
+ [C2211, C2222, C2212],
+ [C1211, C1222, C1212] ],
+ dtype=numpy.float64)
+ stress = numpy.dot(elastic, strain) + initialStress
+ meanStrain = (strain[0,0] + strain[1,0])/3.0
+ viscousStrain = numpy.array([strain[0,0] - meanStrain,
+ strain[1,0] - meanStrain,
+ -meanStrain,
+ strain[2,0]],
+ dtype=numpy.float64)
+ viscousStrainVec = numpy.ravel(viscousStrain)
+ strainVec = numpy.array(strainV, dtype=numpy.float64)
+
+ stateVarsUpdated = numpy.concatenate((stressZZInitial, strainVec,
+ viscousStrainVec, viscousStrainVec,
+ viscousStrainVec))
+ return (elasticConsts, numpy.ravel(stress), numpy.ravel(stateVarsUpdated))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = GenMaxwellPlaneStrainElastic()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.cc 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,404 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application genmaxwellplanestrainelastic.
+
+#include "GenMaxwellPlaneStrainElasticData.hh"
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_dimension = 2;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numLocs = 2;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numProperties = 9;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numStateVars = 5;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numDBProperties = 9;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numDBStateVars = 16;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numPropsQuadPt = 9;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numVarsQuadPt = 16;
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_lengthScale = 1.00000000e+03;
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_timeScale = 1.00000000e+00;
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_pressureScale = 2.25000000e+10;
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_densityScale = 1.00000000e+03;
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_dtStableImplicit = 8.88888889e+06;
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::GenMaxwellPlaneStrainElasticData::_numStateVarValues[] = {
+1,
+3,
+4,
+4,
+4,
+};
+
+const char* pylith::materials::GenMaxwellPlaneStrainElasticData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"shear-ratio-1",
+"shear-ratio-2",
+"shear-ratio-3",
+"viscosity-1",
+"viscosity-2",
+"viscosity-3",
+};
+
+const char* pylith::materials::GenMaxwellPlaneStrainElasticData::_dbStateVarValues[] = {
+"stress-zz-initial",
+"total-strain-xx",
+"total-strain-yy",
+"total-strain-xy",
+"viscous-strain-1-xx",
+"viscous-strain-1-yy",
+"viscous-strain-1-zz",
+"viscous-strain-1-xy",
+"viscous-strain-2-xx",
+"viscous-strain-2-yy",
+"viscous-strain-2-zz",
+"viscous-strain-2-xy",
+"viscous-strain-3-xx",
+"viscous-strain-3-yy",
+"viscous-strain-3-zz",
+"viscous-strain-3-xy",
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 5.00000000e-01,
+ 1.00000000e-01,
+ 2.00000000e-01,
+ 1.00000000e+18,
+ 1.00000000e+17,
+ 1.00000000e+19,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 1.00000000e+18,
+ 1.00000000e+19,
+ 1.00000000e+20,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_dbStateVars[] = {
+ 1.50000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 5.00000000e-01,
+ 1.00000000e-01,
+ 2.00000000e-01,
+ 8.88888889e+07,
+ 4.44444444e+07,
+ 2.22222222e+09,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 1.73611111e+09,
+ 1.73611111e+10,
+ 1.73611111e+11,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_stateVars[] = {
+ 1.50000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 5.00000000e-01,
+ 1.00000000e-01,
+ 2.00000000e-01,
+ 8.88888889e+07,
+ 4.44444444e+07,
+ 2.22222222e+09,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 1.73611111e+09,
+ 1.73611111e+10,
+ 1.73611111e+11,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_stateVarsNondim[] = {
+ 6.66666667e-07,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.00000000e-06,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_strain[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_stress[] = {
+ 7.33350000e+06,
+ 7.73950000e+06,
+ 4.79400000e+06,
+ 4.09740000e+06,
+ 4.15024000e+06,
+ 2.21976000e+06,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_elasticConsts[] = {
+ 6.75000000e+10,
+ 2.25000000e+10,
+ 0.00000000e+00,
+ 2.25000000e+10,
+ 6.75000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 8.64000000e+09,
+ 2.88000000e+09,
+ 0.00000000e+00,
+ 2.88000000e+09,
+ 8.64000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.76000000e+09,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.40000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.40000000e+04,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_initialStrain[] = {
+ 3.10000000e-05,
+ 3.20000000e-05,
+ 3.40000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.40000000e-05,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainElasticData::_stateVarsUpdated[] = {
+ 1.50000000e+04,
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 2.33333333e-05,
+ 3.23333333e-05,
+ -5.56666667e-05,
+ 1.06000000e-04,
+ 2.33333333e-05,
+ 3.23333333e-05,
+ -5.56666667e-05,
+ 1.06000000e-04,
+ 2.33333333e-05,
+ 3.23333333e-05,
+ -5.56666667e-05,
+ 1.06000000e-04,
+ 4.50000000e+04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+ 1.13333333e-04,
+ 1.22333333e-04,
+ -2.35666667e-04,
+ 3.76000000e-04,
+ 1.13333333e-04,
+ 1.22333333e-04,
+ -2.35666667e-04,
+ 3.76000000e-04,
+ 1.13333333e-04,
+ 1.22333333e-04,
+ -2.35666667e-04,
+ 3.76000000e-04,
+};
+
+pylith::materials::GenMaxwellPlaneStrainElasticData::GenMaxwellPlaneStrainElasticData(void)
+{ // constructor
+ dimension = _dimension;
+ numLocs = _numLocs;
+ numProperties = _numProperties;
+ numStateVars = _numStateVars;
+ numDBProperties = _numDBProperties;
+ numDBStateVars = _numDBStateVars;
+ numPropsQuadPt = _numPropsQuadPt;
+ numVarsQuadPt = _numVarsQuadPt;
+ lengthScale = _lengthScale;
+ timeScale = _timeScale;
+ pressureScale = _pressureScale;
+ densityScale = _densityScale;
+ dtStableImplicit = _dtStableImplicit;
+ numPropertyValues = const_cast<int*>(_numPropertyValues);
+ numStateVarValues = const_cast<int*>(_numStateVarValues);
+ dbPropertyValues = const_cast<char**>(_dbPropertyValues);
+ dbStateVarValues = const_cast<char**>(_dbStateVarValues);
+ dbProperties = const_cast<double*>(_dbProperties);
+ dbStateVars = const_cast<double*>(_dbStateVars);
+ properties = const_cast<double*>(_properties);
+ stateVars = const_cast<double*>(_stateVars);
+ propertiesNondim = const_cast<double*>(_propertiesNondim);
+ stateVarsNondim = const_cast<double*>(_stateVarsNondim);
+ density = const_cast<double*>(_density);
+ strain = const_cast<double*>(_strain);
+ stress = const_cast<double*>(_stress);
+ elasticConsts = const_cast<double*>(_elasticConsts);
+ initialStress = const_cast<double*>(_initialStress);
+ initialStrain = const_cast<double*>(_initialStrain);
+ stateVarsUpdated = const_cast<double*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::GenMaxwellPlaneStrainElasticData::~GenMaxwellPlaneStrainElasticData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainElasticData.hh 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,110 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application genmaxwellplanestrainelastic.
+
+#if !defined(pylith_materials_genmaxwellplanestrainelasticdata_hh)
+#define pylith_materials_genmaxwellplanestrainelasticdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class GenMaxwellPlaneStrainElasticData;
+ } // pylith
+} // materials
+
+class pylith::materials::GenMaxwellPlaneStrainElasticData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ GenMaxwellPlaneStrainElasticData(void);
+
+ /// Destructor
+ ~GenMaxwellPlaneStrainElasticData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numLocs;
+
+ static const int _numProperties;
+
+ static const int _numStateVars;
+
+ static const int _numDBProperties;
+
+ static const int _numDBStateVars;
+
+ static const int _numPropsQuadPt;
+
+ static const int _numVarsQuadPt;
+
+ static const double _lengthScale;
+
+ static const double _timeScale;
+
+ static const double _pressureScale;
+
+ static const double _densityScale;
+
+ static const double _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const double _dbProperties[];
+
+ static const double _dbStateVars[];
+
+ static const double _properties[];
+
+ static const double _stateVars[];
+
+ static const double _propertiesNondim[];
+
+ static const double _stateVarsNondim[];
+
+ static const double _density[];
+
+ static const double _strain[];
+
+ static const double _stress[];
+
+ static const double _elasticConsts[];
+
+ static const double _initialStress[];
+
+ static const double _initialStrain[];
+
+ static const double _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_genmaxwellplanestrainelasticdata_hh
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,442 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDep.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ GenMaxwellPlaneStrain object with viscoelastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+import math
+
+# ----------------------------------------------------------------------
+dimension = 2
+numElasticConsts = 9
+tensorSize = 3
+numMaxwellModels = 3
+
+# GenMaxwellPlaneStrainTimeDep class
+class GenMaxwellPlaneStrainTimeDep(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ GenMaxwellPlaneStrain object using viscoelastic behavior.
+ """
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="genmaxwellplanestraintimedep"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ import pdb
+ pdb.set_trace()
+
+ numLocs = 2
+ self.dt = 2.0e5
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "shear-ratio-1", "shear-ratio-2", "shear-ratio-3",
+ "viscosity-1", "viscosity-2", "viscosity-3"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1, 1, 1, 1],
+ dtype=numpy.int32)
+
+ self.dbStateVarValues = ["stress-zz-initial",
+ "total-strain-xx",
+ "total-strain-yy",
+ "total-strain-xy",
+ "viscous-strain-1-xx",
+ "viscous-strain-1-yy",
+ "viscous-strain-1-zz",
+ "viscous-strain-1-xy",
+ "viscous-strain-2-xx",
+ "viscous-strain-2-yy",
+ "viscous-strain-2-zz",
+ "viscous-strain-2-xy",
+ "viscous-strain-3-xx",
+ "viscous-strain-3-yy",
+ "viscous-strain-3-zz",
+ "viscous-strain-3-xy"
+ ]
+ self.numStateVarValues = numpy.array([1, tensorSize, 4, 4, 4],
+ dtype=numpy.int32)
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ shearRatioA = [0.5, 0.1, 0.2]
+ viscosityA = [1.0e18, 1.0e17, 1.0e19]
+ strainA = [1.1e-4, 1.2e-4, 1.4e-4]
+ initialStressA = [2.1e4, 2.2e4, 2.4e4]
+ initialStrainA = [3.6e-5, 3.5e-5, 3.3e-5]
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+ meanStrainA = (strainA[0] + strainA[1])/3.0
+ stressInitialZZA = numpy.array([2.0e4], dtype=numpy.float64)
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ shearRatioB = [0.2, 0.2, 0.2]
+ viscosityB = [1.0e18, 1.0e19, 1.0e20]
+ strainB = [4.1e-4, 4.2e-4, 4.4e-4]
+ initialStressB = [5.1e4, 5.2e4, 5.4e4]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.6e-5]
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ meanStrainB = (strainB[0] + strainB[1])/3.0
+ stressInitialZZB = numpy.array([5.0e4], dtype=numpy.float64)
+
+ # Compute Maxwell time and viscous strain.
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0], dtype=numpy.float64)
+ strainTA = [strainA[0], strainA[1], 0.0, strainA[2]]
+ strainTB = [strainB[0], strainB[1], 0.0, strainB[2]]
+
+ maxwellTimeA = [0.0, 0.0, 0.0]
+ maxwellTimeB = [0.0, 0.0, 0.0]
+ visStrainA = numpy.zeros( (numMaxwellModels, 4), dtype=numpy.float64)
+ visStrainB = numpy.zeros( (numMaxwellModels, 4), dtype=numpy.float64)
+ for imodel in xrange(numMaxwellModels):
+ if shearRatioA[imodel] != 0.0:
+ maxwellTimeA[imodel] = viscosityA[imodel]/(muA*shearRatioA[imodel])
+ visStrainA[imodel,:] = numpy.array(strainTA, dtype=numpy.float64) - \
+ diag * meanStrainA
+ if shearRatioB[imodel] != 0.0:
+ maxwellTimeB[imodel] = viscosityB[imodel]/(muB*shearRatioB[imodel])
+ visStrainB[imodel,:] = numpy.array(strainTB, dtype=numpy.float64) - \
+ diag * meanStrainB
+
+ dbPropA = [densityA, vsA, vpA] + shearRatioA + viscosityA
+ dbPropB = [densityB, vsB, vpB] + shearRatioB + viscosityB
+ self.dbProperties = numpy.array([dbPropA, dbPropB], dtype=numpy.float64)
+ propA = [densityA, muA, lambdaA] + shearRatioA + maxwellTimeA
+ propB = [densityB, muB, lambdaB] + shearRatioB + maxwellTimeB
+ self.properties = numpy.array([propA, propB], dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ # At present, only the first (stressInitialZZ) is being used.
+ self.dbStateVars = numpy.zeros((numLocs,
+ 1 + tensorSize + 4 * numMaxwellModels),
+ dtype=numpy.float64)
+ self.dbStateVars[0, 0] = stressInitialZZA
+ self.dbStateVars[1, 0] = stressInitialZZB
+
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.timeScale = 1.0
+ self.densityScale = 1.0e+3
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ time0 = self.timeScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0,
+ shearRatioA[0], shearRatioA[1], shearRatioA[2],
+ maxwellTimeA[0]/time0, maxwellTimeA[1]/time0,
+ maxwellTimeA[2]/time0],
+ [densityB/density0, muB/mu0, lambdaB/mu0,
+ shearRatioB[0], shearRatioB[1], shearRatioB[2],
+ maxwellTimeB[0]/time0, maxwellTimeB[1]/time0,
+ maxwellTimeB[2]/time0] ],
+ dtype=numpy.float64)
+
+ self.initialStress = numpy.array([initialStressA,
+ initialStressB],
+ dtype=numpy.float64)
+ self.initialStrain = numpy.array([initialStrainA,
+ initialStrainB],
+ dtype=numpy.float64)
+
+ self.density = numpy.array([densityA,
+ densityB],
+ dtype=numpy.float64)
+
+ # Simplest approach for now is to assume this is the first step
+ # after the elastic solution. In that case, both the total strain
+ # from the last step (total_strain) and the total viscous strain
+ # (viscous_strain) are defined by the assigned elastic strain.
+ # Revised approach. For a better test, I am setting the total strain
+ # for the current time step to be equal to the strain from the previous
+ # time step plus a constant amount.
+ totalStrainA = [strainA[0] + 1.0e-5,
+ strainA[1] + 1.0e-5,
+ strainA[2] + 1.0e-5]
+ totalStrainB = [strainB[0] + 1.0e-5,
+ strainB[1] + 1.0e-5,
+ strainB[2] + 1.0e-5]
+
+ visStrainVecA = numpy.ravel(visStrainA)
+ visStrainVecB = numpy.ravel(visStrainB)
+ strainVecA = numpy.array(strainA, dtype=numpy.float64)
+ strainVecB = numpy.array(strainB, dtype=numpy.float64)
+ stateVarsA = numpy.concatenate((stressInitialZZA, strainVecA,
+ visStrainVecA))
+ stateVarsB = numpy.concatenate((stressInitialZZB, strainVecB,
+ visStrainVecB))
+
+ self.stateVars = numpy.array([stateVarsA, stateVarsB], dtype=numpy.float64)
+
+ stressInitialZZANondim = stressInitialZZA/mu0
+ stressInitialZZBNondim = stressInitialZZB/mu0
+ self.stateVarsNondim = numpy.zeros( (numLocs,
+ 1 + tensorSize + 4 * numMaxwellModels),
+ dtype=numpy.float64)
+ self.stateVarsNondim[:] = self.stateVars[:] # no scaling
+ # Scale stressInitialZZ
+ self.stateVarsNondim[0, 0] = stressInitialZZANondim
+ self.stateVarsNondim[1, 0] = stressInitialZZBNondim
+
+ self.strain = numpy.array([totalStrainA, totalStrainB], dtype=numpy.float64)
+ self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.stateVarsUpdated = numpy.zeros((numLocs,
+ 1 + tensorSize + 4 * numMaxwellModels),
+ dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros((numLocs, numElasticConsts), \
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA,
+ shearRatioA, maxwellTimeA,
+ totalStrainA, visStrainA,
+ initialStressA, initialStrainA,
+ stateVarsA)
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB,
+ shearRatioB, maxwellTimeB,
+ totalStrainB, visStrainB,
+ initialStressB, initialStrainB,
+ stateVarsB)
+ self.dtStableImplicit = 0.2*min(min(maxwellTimeA), min(maxwellTimeB))
+
+ return
+
+
+ def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
+ muV, lambdaV, shearRatioV, maxwellTimeV, dqV,
+ totalStrainT, viscousStrainT,
+ initialStress, initialStrain,
+ stateVars):
+ """
+ Function to compute a particular stress component as a function of a
+ strain component.
+ """
+ strainTest = numpy.array(strainTpdt, dtype=numpy.float64)
+ strainTest[strainComp] = strainVal
+ stressTpdt, viscousStrainTpdt = self._computeStress(strainTest,
+ muV,
+ lambdaV,
+ shearRatioV,
+ maxwellTimeV,
+ dqV,
+ totalStrainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain,
+ stateVars)
+ return stressTpdt[stressComp]
+
+
+ def _computeStress(self, strainTpdt, muV, lambdaV, shearRatioV, maxwellTimeV,
+ dqV, strainT, viscousStrainT,
+ initialStress, initialStrain, stateVars):
+ """
+ Function to compute stresses and viscous strains for a given strain.
+ """
+
+ bulkModulus = lambdaV + 2.0 * muV / 3.0
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0], dtype=numpy.float64)
+
+ # Initial stresses and strains
+ meanStrainInitial = (initialStrain[0] + initialStrain[1]) / 3.0
+ meanStressInitial = \
+ (initialStress[0] + initialStress[1] + stateVars[0]) / 3.0
+
+ initialStressV = numpy.array([initialStress[0],
+ initialStress[1],
+ stateVars[0],
+ initialStress[2]], dtype=numpy.float64)
+ initialStrainV = numpy.array([initialStrain[0],
+ initialStrain[1],
+ 0.0,
+ initialStrain[2]], dtype=numpy.float64)
+
+ devStrainInitial = initialStrainV - diag * meanStrainInitial
+ devStressInitial = initialStressV - diag * meanStressInitial
+
+ # Strains from previous time step
+ strainTV = numpy.array([strainT[0], strainT[1], 0.0, strainT[2]],
+ dtype=numpy.float64)
+ meanStrainT = (strainT[0] + strainT[1]) / 3.0
+ devStrainT = strainTV - diag * meanStrainT - devStrainInitial
+
+ # Strains and mean stress for current time step
+ strainTpdtV = numpy.array([strainTpdt[0], strainTpdt[1], 0.0,
+ strainTpdt[2]], dtype=numpy.float64)
+ meanStrainTpdt = (strainTpdt[0] + strainTpdt[1]) / 3.0
+ devStrainTpdt = strainTpdtV - diag * meanStrainTpdt - devStrainInitial
+ meanStressTpdt = 3.0 * bulkModulus * \
+ (meanStrainTpdt - meanStrainInitial) + meanStressInitial
+
+ # Compute viscous factors
+ visFrac = 0.0
+ visFac = 0.0
+ expFac = numpy.zeros((numMaxwellModels), dtype=numpy.float64)
+ for model in range(numMaxwellModels):
+ visFrac += shearRatioV[model]
+ visFac += shearRatioV[model] * dqV[model]
+ expFac[model] = math.exp(-self.dt/maxwellTimeV[model])
+
+
+ elasFrac = 1.0 - visFrac
+ stressTpdt = numpy.zeros( (4), dtype=numpy.float64)
+ viscousStrainTpdt = numpy.zeros( (3, 4), dtype=numpy.float64)
+ elasFac = 2.0 * muV
+
+ # Compute viscous strain and stress
+ for iComp in range(4):
+ deltaStrain = devStrainTpdt[iComp] - devStrainT[iComp]
+ devStressTpdt = elasFrac * devStrainTpdt[iComp]
+ for model in range(numMaxwellModels):
+ if shearRatioV[model] != 0.0:
+ viscousStrainTpdt[model, iComp] = expFac[model] * \
+ viscousStrainT[model, iComp] + \
+ dqV[model] * deltaStrain
+ devStressTpdt += shearRatioV[model] * viscousStrainTpdt[model, iComp]
+ devStressTpdt = elasFac * devStressTpdt + devStressInitial[iComp]
+ stressTpdt[iComp] = diag[iComp] * (meanStressTpdt + meanStressInitial) + \
+ devStressTpdt
+
+ stressTpdtV = numpy.array([stressTpdt[0], stressTpdt[1], stressTpdt[3]],
+ dtype=numpy.float64)
+
+ return stressTpdtV, viscousStrainTpdt
+
+
+ def _computeViscousFactor(self, maxwellTime):
+ """
+ Compute viscous strain factor for a given Maxwell time.
+ """
+
+ timeFrac = 1.0e-5
+ numTerms = 5
+ dq = 0.0
+ if maxwellTime < timeFrac*self.dt:
+ fSign = 1.0
+ factorial = 1.0
+ fraction = 1.0
+ dq = 1.0
+ for iTerm in range(2, numTerms + 1):
+ factorial *= iTerm
+ fSign *= -1.0
+ fraction *= self.dt/maxwellTime
+ dq += fSign*fraction/factorial
+ else:
+ dq = maxwellTime*(1.0-math.exp(-self.dt/maxwellTime))/self.dt
+
+ return dq
+
+
+ def _calcStress(self, strainV, muV, lambdaV, shearRatioV, maxwellTimeV,
+ totalStrainV, viscousStrainV, initialStressV, initialStrainV,
+ stateVars):
+ """
+ Compute stress and derivative of elasticity matrix.
+ This assumes behavior is always viscoelastic.
+ """
+ import scipy.misc
+
+ # Define some numpy arrays
+ strainTpdt = numpy.array(totalStrainV, dtype=numpy.float64)
+ strainT = numpy.array(strainV, dtype=numpy.float64)
+ viscousStrainT = numpy.array(viscousStrainV, dtype=numpy.float64)
+ initialStress = numpy.array(initialStressV, dtype=numpy.float64)
+ initialStrain = numpy.array(initialStrainV, dtype=numpy.float64)
+ stressZZInitial = numpy.array([stateVars[0]], dtype=numpy.float64)
+
+ # Compute viscous factors for each Maxwell model
+ dqV = numpy.zeros(numMaxwellModels, dtype=numpy.float64)
+ for model in range(numMaxwellModels):
+ dqV[model] = self._computeViscousFactor(maxwellTimeV[model])
+
+ # Compute current stress and viscous strain
+ stressTpdt, viscousStrainTpdt = self._computeStress(strainTpdt,
+ muV,
+ lambdaV,
+ shearRatioV,
+ maxwellTimeV,
+ dqV,
+ strainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain,
+ stateVars)
+ # Form updated state variables
+ strainTpdtVec = numpy.ravel(strainTpdt)
+ viscousStrainTpdtVec = numpy.ravel(viscousStrainTpdt)
+ stateVarsUpdated = numpy.concatenate((stressZZInitial, strainTpdtVec,
+ viscousStrainTpdtVec))
+
+ # Compute components of tangent constitutive matrix using numerical
+ # derivatives.
+ derivDx = 1.0e-12
+ derivOrder = 3
+ elasticConstsList = []
+
+ for stressComp in range(tensorSize):
+ for strainComp in range(tensorSize):
+ dStressDStrain = stressComp + strainComp
+ dStressDStrain = scipy.misc.derivative(self._calcStressComponent,
+ strainTpdt[strainComp],
+ dx=derivDx,
+ args=(strainComp,
+ stressComp,
+ strainTpdt,
+ muV,
+ lambdaV,
+ shearRatioV,
+ maxwellTimeV,
+ dqV,
+ strainT,
+ viscousStrainT,
+ initialStress,
+ initialStrain,
+ stateVars),
+ order=derivOrder)
+ elasticConstsList.append(dStressDStrain)
+
+ elasticConsts = numpy.array(elasticConstsList, dtype=numpy.float64)
+
+ return (elasticConsts, numpy.ravel(stressTpdt),
+ numpy.ravel(stateVarsUpdated))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = GenMaxwellPlaneStrainTimeDep()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.cc 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,404 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application genmaxwellplanestraintimedep.
+
+#include "GenMaxwellPlaneStrainTimeDepData.hh"
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_dimension = 2;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numLocs = 2;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numProperties = 9;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numStateVars = 5;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numDBProperties = 9;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numDBStateVars = 16;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numPropsQuadPt = 9;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numVarsQuadPt = 16;
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_lengthScale = 1.00000000e+03;
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_timeScale = 1.00000000e+00;
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_pressureScale = 2.25000000e+10;
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_densityScale = 1.00000000e+03;
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_dtStableImplicit = 8.88888889e+06;
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::GenMaxwellPlaneStrainTimeDepData::_numStateVarValues[] = {
+1,
+3,
+4,
+4,
+4,
+};
+
+const char* pylith::materials::GenMaxwellPlaneStrainTimeDepData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"shear-ratio-1",
+"shear-ratio-2",
+"shear-ratio-3",
+"viscosity-1",
+"viscosity-2",
+"viscosity-3",
+};
+
+const char* pylith::materials::GenMaxwellPlaneStrainTimeDepData::_dbStateVarValues[] = {
+"stress-zz-initial",
+"total-strain-xx",
+"total-strain-yy",
+"total-strain-xy",
+"viscous-strain-1-xx",
+"viscous-strain-1-yy",
+"viscous-strain-1-zz",
+"viscous-strain-1-xy",
+"viscous-strain-2-xx",
+"viscous-strain-2-yy",
+"viscous-strain-2-zz",
+"viscous-strain-2-xy",
+"viscous-strain-3-xx",
+"viscous-strain-3-yy",
+"viscous-strain-3-zz",
+"viscous-strain-3-xy",
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 5.00000000e-01,
+ 1.00000000e-01,
+ 2.00000000e-01,
+ 1.00000000e+18,
+ 1.00000000e+17,
+ 1.00000000e+19,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 1.00000000e+18,
+ 1.00000000e+19,
+ 1.00000000e+20,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_dbStateVars[] = {
+ 2.00000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.00000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 5.00000000e-01,
+ 1.00000000e-01,
+ 2.00000000e-01,
+ 8.88888889e+07,
+ 4.44444444e+07,
+ 2.22222222e+09,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 1.73611111e+09,
+ 1.73611111e+10,
+ 1.73611111e+11,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_stateVars[] = {
+ 2.00000000e+04,
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 3.33333333e-05,
+ 4.33333333e-05,
+ -7.66666667e-05,
+ 1.40000000e-04,
+ 3.33333333e-05,
+ 4.33333333e-05,
+ -7.66666667e-05,
+ 1.40000000e-04,
+ 3.33333333e-05,
+ 4.33333333e-05,
+ -7.66666667e-05,
+ 1.40000000e-04,
+ 5.00000000e+04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+ 1.33333333e-04,
+ 1.43333333e-04,
+ -2.76666667e-04,
+ 4.40000000e-04,
+ 1.33333333e-04,
+ 1.43333333e-04,
+ -2.76666667e-04,
+ 4.40000000e-04,
+ 1.33333333e-04,
+ 1.43333333e-04,
+ -2.76666667e-04,
+ 4.40000000e-04,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 5.00000000e-01,
+ 1.00000000e-01,
+ 2.00000000e-01,
+ 8.88888889e+07,
+ 4.44444444e+07,
+ 2.22222222e+09,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 2.00000000e-01,
+ 1.73611111e+09,
+ 1.73611111e+10,
+ 1.73611111e+11,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_stateVarsNondim[] = {
+ 8.88888889e-07,
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 3.33333333e-05,
+ 4.33333333e-05,
+ -7.66666667e-05,
+ 1.40000000e-04,
+ 3.33333333e-05,
+ 4.33333333e-05,
+ -7.66666667e-05,
+ 1.40000000e-04,
+ 3.33333333e-05,
+ 4.33333333e-05,
+ -7.66666667e-05,
+ 1.40000000e-04,
+ 2.22222222e-06,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+ 1.33333333e-04,
+ 1.43333333e-04,
+ -2.76666667e-04,
+ 4.40000000e-04,
+ 1.33333333e-04,
+ 1.43333333e-04,
+ -2.76666667e-04,
+ 4.40000000e-04,
+ 1.33333333e-04,
+ 1.43333333e-04,
+ -2.76666667e-04,
+ 4.40000000e-04,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_strain[] = {
+ 1.20000000e-04,
+ 1.30000000e-04,
+ 1.50000000e-04,
+ 4.20000000e-04,
+ 4.30000000e-04,
+ 4.50000000e-04,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_stress[] = {
+ 8.29099455e+06,
+ 8.75027873e+06,
+ 6.46662036e+06,
+ 4.33270011e+06,
+ 4.38899464e+06,
+ 2.49387045e+06,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_elasticConsts[] = {
+ 6.74761273e+10,
+ 2.25119358e+10,
+ 0.00000000e+00,
+ 2.25119358e+10,
+ 6.74761273e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.49641915e+10,
+ 8.63995077e+09,
+ 2.88002472e+09,
+ 0.00000000e+00,
+ 2.88002472e+09,
+ 8.63995124e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.75992628e+09,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.40000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.40000000e+04,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_initialStrain[] = {
+ 3.60000000e-05,
+ 3.50000000e-05,
+ 3.30000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.60000000e-05,
+};
+
+const double pylith::materials::GenMaxwellPlaneStrainTimeDepData::_stateVarsUpdated[] = {
+ 2.00000000e+04,
+ 1.20000000e-04,
+ 1.30000000e-04,
+ 1.50000000e-04,
+ 3.65880038e-05,
+ 4.65655291e-05,
+ -8.31535329e-05,
+ 1.49674113e-04,
+ 3.65095149e-05,
+ 4.64646160e-05,
+ -8.29741309e-05,
+ 1.49348949e-04,
+ 3.66635168e-05,
+ 4.66626168e-05,
+ -8.33261337e-05,
+ 1.49986951e-04,
+ 5.00000000e+04,
+ 4.20000000e-04,
+ 4.30000000e-04,
+ 4.50000000e-04,
+ 1.36651116e-04,
+ 1.46649964e-04,
+ -2.83301079e-04,
+ 4.49948739e-04,
+ 1.36665111e-04,
+ 1.46664996e-04,
+ -2.83330108e-04,
+ 4.49994874e-04,
+ 1.36666511e-04,
+ 1.46666500e-04,
+ -2.83333011e-04,
+ 4.49999487e-04,
+};
+
+pylith::materials::GenMaxwellPlaneStrainTimeDepData::GenMaxwellPlaneStrainTimeDepData(void)
+{ // constructor
+ dimension = _dimension;
+ numLocs = _numLocs;
+ numProperties = _numProperties;
+ numStateVars = _numStateVars;
+ numDBProperties = _numDBProperties;
+ numDBStateVars = _numDBStateVars;
+ numPropsQuadPt = _numPropsQuadPt;
+ numVarsQuadPt = _numVarsQuadPt;
+ lengthScale = _lengthScale;
+ timeScale = _timeScale;
+ pressureScale = _pressureScale;
+ densityScale = _densityScale;
+ dtStableImplicit = _dtStableImplicit;
+ numPropertyValues = const_cast<int*>(_numPropertyValues);
+ numStateVarValues = const_cast<int*>(_numStateVarValues);
+ dbPropertyValues = const_cast<char**>(_dbPropertyValues);
+ dbStateVarValues = const_cast<char**>(_dbStateVarValues);
+ dbProperties = const_cast<double*>(_dbProperties);
+ dbStateVars = const_cast<double*>(_dbStateVars);
+ properties = const_cast<double*>(_properties);
+ stateVars = const_cast<double*>(_stateVars);
+ propertiesNondim = const_cast<double*>(_propertiesNondim);
+ stateVarsNondim = const_cast<double*>(_stateVarsNondim);
+ density = const_cast<double*>(_density);
+ strain = const_cast<double*>(_strain);
+ stress = const_cast<double*>(_stress);
+ elasticConsts = const_cast<double*>(_elasticConsts);
+ initialStress = const_cast<double*>(_initialStress);
+ initialStrain = const_cast<double*>(_initialStrain);
+ stateVarsUpdated = const_cast<double*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::GenMaxwellPlaneStrainTimeDepData::~GenMaxwellPlaneStrainTimeDepData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellPlaneStrainTimeDepData.hh 2011-05-29 23:59:03 UTC (rev 18489)
@@ -0,0 +1,110 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application genmaxwellplanestraintimedep.
+
+#if !defined(pylith_materials_genmaxwellplanestraintimedepdata_hh)
+#define pylith_materials_genmaxwellplanestraintimedepdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class GenMaxwellPlaneStrainTimeDepData;
+ } // pylith
+} // materials
+
+class pylith::materials::GenMaxwellPlaneStrainTimeDepData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ GenMaxwellPlaneStrainTimeDepData(void);
+
+ /// Destructor
+ ~GenMaxwellPlaneStrainTimeDepData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numLocs;
+
+ static const int _numProperties;
+
+ static const int _numStateVars;
+
+ static const int _numDBProperties;
+
+ static const int _numDBStateVars;
+
+ static const int _numPropsQuadPt;
+
+ static const int _numVarsQuadPt;
+
+ static const double _lengthScale;
+
+ static const double _timeScale;
+
+ static const double _pressureScale;
+
+ static const double _densityScale;
+
+ static const double _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const double _dbProperties[];
+
+ static const double _dbStateVars[];
+
+ static const double _properties[];
+
+ static const double _stateVars[];
+
+ static const double _propertiesNondim[];
+
+ static const double _stateVarsNondim[];
+
+ static const double _density[];
+
+ static const double _strain[];
+
+ static const double _stress[];
+
+ static const double _elasticConsts[];
+
+ static const double _initialStress[];
+
+ static const double _initialStrain[];
+
+ static const double _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_genmaxwellplanestraintimedepdata_hh
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElastic.py 2011-05-29 23:59:03 UTC (rev 18489)
@@ -77,6 +77,7 @@
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
maxwellTimeA = viscosityA / muA
+ stressInitialZZA = numpy.array([1.5e4], dtype=numpy.float64)
densityB = 2000.0
vsB = 1200.0
@@ -88,6 +89,7 @@
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
maxwellTimeB = viscosityB / muB
+ stressInitialZZB = numpy.array([4.5e4], dtype=numpy.float64)
self.lengthScale = 1.0e+3
self.pressureScale = muA
@@ -102,21 +104,36 @@
dtype=numpy.float64)
# TEMPORARY, need to determine how to use initial state variables
- self.dbStateVars = numpy.zeros( (numLocs, tensorSize+5),
+ # At present, only the first (stressInitialZZ) is being used.
+ self.dbStateVars = numpy.zeros( (numLocs, 1 + tensorSize + 4),
dtype=numpy.float64)
- self.stateVars = numpy.zeros( (numLocs, tensorSize+5),
+ self.dbStateVars[0, 0] = stressInitialZZA
+ self.dbStateVars[1, 0] = stressInitialZZB
+
+ self.stateVars = numpy.zeros( (numLocs, 1 + tensorSize + 4),
dtype=numpy.float64)
+ self.stateVars[0, 0] = stressInitialZZA
+ self.stateVars[1, 0] = stressInitialZZB
mu0 = self.pressureScale
density0 = self.densityScale
time0 = self.timeScale
self.propertiesNondim = \
- numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, maxwellTimeA/time0],
- [densityB/density0, muB/mu0, lambdaB/mu0, maxwellTimeB/time0] ],
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0,
+ maxwellTimeA/time0],
+ [densityB/density0, muB/mu0, lambdaB/mu0,
+ maxwellTimeB/time0] ],
dtype=numpy.float64)
- self.stateVarsNondim = self.stateVars # no scaling
+ stressInitialZZANondim = stressInitialZZA/mu0
+ stressInitialZZBNondim = stressInitialZZB/mu0
+ self.stateVarsNondim = numpy.zeros( (numLocs, 1 + tensorSize + 4),
+ dtype=numpy.float64)
+
+ self.stateVarsNondim[0, 0] = stressInitialZZANondim
+ self.stateVarsNondim[1, 0] = stressInitialZZBNondim
+
self.initialStress = numpy.array([initialStressA,
initialStressB],
dtype=numpy.float64)
@@ -135,32 +152,35 @@
self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts), \
dtype=numpy.float64)
- self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + 5), \
+ self.stateVarsUpdated = numpy.zeros( (numLocs, 1 + tensorSize + 4), \
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
- self._calcStress(strainA, muA, lambdaA, \
- initialStressA, initialStrainA)
+ self._calcStress(strainA, muA, lambdaA,
+ initialStressA, initialStrainA,
+ self.stateVars[0,:])
(self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
- self._calcStress(strainB, muB, lambdaB, \
- initialStressB, initialStrainB)
+ self._calcStress(strainB, muB, lambdaB,
+ initialStressB, initialStrainB,
+ self.stateVars[1,:])
self.dtStableImplicit = 0.2*min(maxwellTimeA, maxwellTimeB)
return
- def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV):
+ def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV,
+ stateVars):
"""
Compute stress and derivative of elasticity matrix.
"""
- C1111 = lambdaV + 2.0*muV
+ C1111 = lambdaV + 2.0 * muV
C1122 = lambdaV
C1112 = 0.0
C2211 = lambdaV
- C2222 = lambdaV + 2.0*muV
+ C2222 = lambdaV + 2.0 * muV
C2212 = 0.0
C1211 = 0.0
C1222 = 0.0
- C1212 = 2.0*muV
+ C1212 = 2.0 * muV
elasticConsts = numpy.array([C1111, C1122, C1112,
C2211, C2222, C2212,
C1211, C1222, C1212], dtype=numpy.float64)
@@ -168,6 +188,8 @@
strain = numpy.reshape(strainV, (tensorSize,1))
initialStress = numpy.reshape(initialStressV, (tensorSize,1))
initialStrain = numpy.reshape(initialStrainV, (tensorSize,1))
+ stressZZInitial = numpy.array([stateVars[0]], dtype=numpy.float64)
+
elastic = numpy.array([ [C1111, C1122, C1112],
[C2211, C2222, C2212],
[C1211, C1222, C1212] ], dtype=numpy.float64)
@@ -180,7 +202,6 @@
dtype=numpy.float64)
viscousStrainVec = numpy.ravel(viscousStrain)
strainVec = numpy.ravel(strain)
- stressZZInitial = numpy.zeros(1, dtype=numpy.float64)
stateVarsUpdated = numpy.concatenate((stressZZInitial, strainVec,
viscousStrainVec))
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainElasticData.cc 2011-05-29 23:59:03 UTC (rev 18489)
@@ -90,6 +90,7 @@
};
const double pylith::materials::MaxwellPlaneStrainElasticData::_dbStateVars[] = {
+ 1.50000000e+04,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
@@ -97,6 +98,7 @@
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
+ 4.50000000e+04,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
@@ -104,8 +106,6 @@
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 0.00000000e+00,
- 0.00000000e+00,
};
const double pylith::materials::MaxwellPlaneStrainElasticData::_properties[] = {
@@ -120,6 +120,7 @@
};
const double pylith::materials::MaxwellPlaneStrainElasticData::_stateVars[] = {
+ 1.50000000e+04,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
@@ -127,6 +128,7 @@
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
+ 4.50000000e+04,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
@@ -134,8 +136,6 @@
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 0.00000000e+00,
- 0.00000000e+00,
};
const double pylith::materials::MaxwellPlaneStrainElasticData::_propertiesNondim[] = {
@@ -150,6 +150,7 @@
};
const double pylith::materials::MaxwellPlaneStrainElasticData::_stateVarsNondim[] = {
+ 6.66666667e-07,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
@@ -157,6 +158,7 @@
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
+ 2.00000000e-06,
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
@@ -164,8 +166,6 @@
0.00000000e+00,
0.00000000e+00,
0.00000000e+00,
- 0.00000000e+00,
- 0.00000000e+00,
};
const double pylith::materials::MaxwellPlaneStrainElasticData::_density[] = {
@@ -231,7 +231,7 @@
};
const double pylith::materials::MaxwellPlaneStrainElasticData::_stateVarsUpdated[] = {
- 0.00000000e+00,
+ 1.50000000e+04,
1.10000000e-04,
1.20000000e-04,
1.40000000e-04,
@@ -239,7 +239,7 @@
4.33333333e-05,
-7.66666667e-05,
1.40000000e-04,
- 0.00000000e+00,
+ 4.50000000e+04,
4.10000000e-04,
4.20000000e-04,
4.40000000e-04,
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellPlaneStrainTimeDep.py 2011-05-29 23:59:03 UTC (rev 18489)
@@ -24,6 +24,7 @@
from ElasticMaterialApp import ElasticMaterialApp
import numpy
+import math
# ----------------------------------------------------------------------
dimension = 2
@@ -66,7 +67,6 @@
"viscous-strain-zz",
"viscous-strain-xy"
]
- self.stateVarValues = ["stress-zz-initial","total-strain", "viscous-strain"]
self.numStateVarValues = numpy.array([1, 3, 4], dtype=numpy.int32)
self.dt = 2.0e5
@@ -97,8 +97,7 @@
meanStrainB = (strainB[0] + strainB[1])/3.0
stressInitialZZB = numpy.array([5.0e4], dtype=numpy.float64)
- diag = numpy.array([1.0, 1.0, 1.0, 0.0],
- dtype=numpy.float64)
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0], dtype=numpy.float64)
strainTA = [strainA[0], strainA[1], 0.0, strainA[2]]
strainTB = [strainB[0], strainB[1], 0.0, strainB[2]]
@@ -117,7 +116,7 @@
# TEMPORARY, need to determine how to use initial state variables
# At present, only the first (stressInitialZZ) is being used.
- self.dbStateVars = numpy.zeros( (numLocs, tensorSize+5),
+ self.dbStateVars = numpy.zeros( (numLocs, 1 + tensorSize + 4),
dtype=numpy.float64)
self.dbStateVars[0, 0] = stressInitialZZA
self.dbStateVars[1, 0] = stressInitialZZB
@@ -156,32 +155,32 @@
totalStrainB = [strainB[0] + 1.0e-5,
strainB[1] + 1.0e-5,
strainB[2] + 1.0e-5]
+
viscousStrainA = numpy.array(strainTA) - diag * meanStrainA
viscousStrainB = numpy.array(strainTB) - diag * meanStrainB
viscousStrainVecA = numpy.ravel(viscousStrainA)
viscousStrainVecB = numpy.ravel(viscousStrainB)
strainVecA = numpy.array(strainA, dtype=numpy.float64)
strainVecB = numpy.array(strainB, dtype=numpy.float64)
+
stateVarsA = numpy.concatenate((stressInitialZZA, strainVecA,
viscousStrainVecA))
stateVarsB = numpy.concatenate((stressInitialZZB, strainVecB,
viscousStrainVecB))
- self.stateVars = numpy.array([ stateVarsA,
- stateVarsB ],
- dtype=numpy.float64)
+ self.stateVars = numpy.array([stateVarsA, stateVarsB], dtype=numpy.float64)
+
stressInitialZZANondim = stressInitialZZA/mu0
stressInitialZZBNondim = stressInitialZZB/mu0
- self.stateVarsNondim = numpy.zeros( (numLocs, tensorSize + 5),
+ self.stateVarsNondim = numpy.zeros( (numLocs, 1 + tensorSize + 4),
dtype=numpy.float64)
self.stateVarsNondim[:] = self.stateVars[:] # no scaling
# Scale stressInitialZZ
self.stateVarsNondim[0, 0] = stressInitialZZANondim
self.stateVarsNondim[1, 0] = stressInitialZZBNondim
- self.strain = numpy.array([totalStrainA, totalStrainB],
- dtype=numpy.float64)
+ self.strain = numpy.array([totalStrainA, totalStrainB], dtype=numpy.float64)
self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
- self.stateVarsUpdated = numpy.zeros( (numLocs, tensorSize + 5),
+ self.stateVarsUpdated = numpy.zeros( (numLocs, 1 + tensorSize + 4),
dtype=numpy.float64)
self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts),
dtype=numpy.float64)
@@ -199,9 +198,6 @@
initialStressB, initialStrainB,
stateVarsB)
- self.stateVarsUpdated[0, 0] = stressInitialZZA
- self.stateVarsUpdated[1, 0] = stressInitialZZB
-
self.dtStableImplicit = 0.2*min(maxwellTimeA, maxwellTimeB)
return
@@ -235,14 +231,12 @@
"""
Function to compute stresses and viscous strains for a given strain.
"""
- import math
bulkModulus = lambdaV + 2.0 * muV / 3.0
diag = [1.0, 1.0, 0.0]
# Initial stresses and strains
- meanStrainInitial = \
- (initialStrain[0] + initialStrain[1]) / 3.0
+ meanStrainInitial = (initialStrain[0] + initialStrain[1]) / 3.0
meanStressInitial = \
(initialStress[0] + initialStress[1] + stateVars[0]) / 3.0
@@ -301,7 +295,6 @@
"""
Compute viscous strain factor for a given Maxwell time.
"""
- import math
timeFrac = 1.0e-5
numTerms = 5
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh 2011-05-29 15:46:39 UTC (rev 18488)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh 2011-05-29 23:59:03 UTC (rev 18489)
@@ -73,6 +73,11 @@
--data.object=MaxwellPlaneStrainElasticData \
--data.parent=ElasticMaterialData
+ python GenMaxwellPlaneStrainElastic.py \
+ --data.namespace=pylith,materials \
+ --data.object=GenMaxwellPlaneStrainElasticData \
+ --data.parent=ElasticMaterialData
+
# 1-D ----------------------------------------------------------------
python ElasticStrain1D.py \
@@ -105,6 +110,11 @@
--data.object=MaxwellIsotropic3DTimeDepData \
--data.parent=ElasticMaterialData
+ python GenMaxwellPlaneStrainTimeDep.py \
+ --data.namespace=pylith,materials \
+ --data.object=GenMaxwellPlaneStrainTimeDepData \
+ --data.parent=ElasticMaterialData
+
python MaxwellPlaneStrainTimeDep.py \
--data.namespace=pylith,materials \
--data.object=MaxwellPlaneStrainTimeDepData \
More information about the CIG-COMMITS
mailing list