[cig-commits] r19599 - in short/3D/PyLith/branches/pylith-scecdynrup: . libsrc/pylith libsrc/pylith/materials modulesrc/materials pylith pylith/materials unittests/libtests/materials unittests/libtests/materials/data unittests/pytests/materials
brad at geodynamics.org
brad at geodynamics.org
Tue Feb 7 13:21:03 PST 2012
Author: brad
Date: 2012-02-07 13:21:03 -0800 (Tue, 07 Feb 2012)
New Revision: 19599
Added:
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i
short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPragerPlaneStrain.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElastic.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/TODO
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.icc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/materialsfwd.hh
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/materials.i
short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/generate.sh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/testmaterials.py
Log:
Merge from tunk.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-02-07 21:21:03 UTC (rev 19599)
@@ -2,6 +2,12 @@
CURRENT ISSUES/PRIORITIES (1.7.0)
======================================================================
+* Manual
+
+ - Order of tensor components for Xdmf files
+ - Drucker Prager fit to yield surface
+ - Slip-weakening healing parameter
+
* configure
+ Check compatibility of PyLith options with PETSc
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/Makefile.am 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/Makefile.am 2012-02-07 21:21:03 UTC (rev 19599)
@@ -121,6 +121,7 @@
materials/MaxwellPlaneStrain.cc \
materials/PowerLaw3D.cc \
materials/DruckerPrager3D.cc \
+ materials/DruckerPragerPlaneStrain.cc \
meshio/BinaryIO.cc \
meshio/GMVFile.cc \
meshio/GMVFileAscii.cc \
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc 2012-02-07 21:21:03 UTC (rev 19599)
@@ -216,7 +216,7 @@
const PylithScalar cohesion = dbValues[db_cohesion];
const PylithScalar dilatationAngle = dbValues[db_dilatationAngle];
- const double pi = M_PI;
+ const PylithScalar pi = M_PI;
if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || cohesion < 0.0 ||
frictionAngle < 0.0 || frictionAngle > pi/2 ||
Copied: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc (from rev 19598, short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,1129 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "DruckerPragerPlaneStrain.hh" // implementation of object methods
+
+#include "Metadata.hh" // USES Metadata
+
+#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include "petsc.h" // USES PetscLogFlops
+
+#include <cmath> // USES fabs()
+#include <cassert> // USES assert()
+#include <cstring> // USES memcpy()
+#include <sstream> // USES std::ostringstream
+#include <iostream> // USES std::cout
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+namespace pylith {
+ namespace materials {
+ namespace _DruckerPragerPlaneStrain{
+
+ /// 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 = 6;
+
+ /// Physical properties.
+ const Metadata::ParamDescription properties[] = {
+ { "density", 1, pylith::topology::FieldBase::SCALAR },
+ { "mu", 1, pylith::topology::FieldBase::SCALAR },
+ { "lambda", 1, pylith::topology::FieldBase::SCALAR },
+ { "alpha_yield", 1, pylith::topology::FieldBase::SCALAR },
+ { "beta", 1, pylith::topology::FieldBase::SCALAR },
+ { "alpha_flow", 1, pylith::topology::FieldBase::SCALAR }
+ };
+
+ // Values expected in properties spatial database
+ const int numDBProperties = 6;
+ const char* dbProperties[6] = {
+ "density",
+ "vs",
+ "vp" ,
+ "friction-angle",
+ "cohesion",
+ "dilatation-angle"
+ };
+
+ /// Number of state variables.
+ const int numStateVars = 2;
+
+ /// State variables.
+ const Metadata::ParamDescription stateVars[] = {
+ { "stress_zz_initial", 1, pylith::topology::FieldBase::SCALAR },
+ { "plastic_strain", 4, pylith::topology::FieldBase::OTHER }
+ };
+
+ // Values expected in state variables spatial database.
+ const int numDBStateVars = 5;
+ const char* dbStateVars[5] = {
+ "stress-zz-initial",
+ "plastic-strain-xx",
+ "plastic-strain-yy",
+ "plastic-strain-zz",
+ "plastic-strain-xy"
+ };
+
+ } // _DruckerPragerPlaneStrain
+ } // materials
+} // pylith
+
+// Indices of physical properties.
+const int pylith::materials::DruckerPragerPlaneStrain::p_density = 0;
+
+const int pylith::materials::DruckerPragerPlaneStrain::p_mu =
+ pylith::materials::DruckerPragerPlaneStrain::p_density + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::p_lambda =
+ pylith::materials::DruckerPragerPlaneStrain::p_mu + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::p_alphaYield =
+ pylith::materials::DruckerPragerPlaneStrain::p_lambda + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::p_beta =
+ pylith::materials::DruckerPragerPlaneStrain::p_alphaYield + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::p_alphaFlow =
+ pylith::materials::DruckerPragerPlaneStrain::p_beta + 1;
+
+// Indices of property database values (order must match dbProperties).
+const int pylith::materials::DruckerPragerPlaneStrain::db_density = 0;
+
+const int pylith::materials::DruckerPragerPlaneStrain::db_vs =
+ pylith::materials::DruckerPragerPlaneStrain::db_density + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::db_vp =
+ pylith::materials::DruckerPragerPlaneStrain::db_vs + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::db_frictionAngle =
+ pylith::materials::DruckerPragerPlaneStrain::db_vp + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::db_cohesion =
+ pylith::materials::DruckerPragerPlaneStrain::db_frictionAngle + 1;
+
+const int pylith::materials::DruckerPragerPlaneStrain::db_dilatationAngle =
+ pylith::materials::DruckerPragerPlaneStrain::db_cohesion + 1;
+
+// Indices of state variables.
+const int pylith::materials::DruckerPragerPlaneStrain::s_stressZZInitial = 0;
+
+const int pylith::materials::DruckerPragerPlaneStrain::s_plasticStrain =
+ pylith::materials::DruckerPragerPlaneStrain::s_stressZZInitial + 1;
+
+// Indices of state variable database values (order must match dbStateVars).
+const int pylith::materials::DruckerPragerPlaneStrain::db_stressZZInitial = 0;
+
+const int pylith::materials::DruckerPragerPlaneStrain::db_plasticStrain =
+ pylith::materials::DruckerPragerPlaneStrain::db_stressZZInitial + 1;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::materials::DruckerPragerPlaneStrain::DruckerPragerPlaneStrain(void) :
+ ElasticMaterial(_DruckerPragerPlaneStrain::dimension,
+ _DruckerPragerPlaneStrain::tensorSize,
+ _DruckerPragerPlaneStrain::numElasticConsts,
+ Metadata(_DruckerPragerPlaneStrain::properties,
+ _DruckerPragerPlaneStrain::numProperties,
+ _DruckerPragerPlaneStrain::dbProperties,
+ _DruckerPragerPlaneStrain::numDBProperties,
+ _DruckerPragerPlaneStrain::stateVars,
+ _DruckerPragerPlaneStrain::numStateVars,
+ _DruckerPragerPlaneStrain::dbStateVars,
+ _DruckerPragerPlaneStrain::numDBStateVars)),
+ _fitMohrCoulomb(MOHR_COULOMB_INSCRIBED),
+ _calcElasticConstsFn(0),
+ _calcStressFn(0),
+ _updateStateVarsFn(0)
+{ // constructor
+ useElasticBehavior(true);
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::materials::DruckerPragerPlaneStrain::~DruckerPragerPlaneStrain(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set fit to Mohr-Coulomb surface.
+void
+pylith::materials::DruckerPragerPlaneStrain::fitMohrCoulomb(
+ FitMohrCoulombEnum value)
+{ // fitMohrCoulomb
+ _fitMohrCoulomb = value;
+} // fitMohrCoulomb
+
+// ----------------------------------------------------------------------
+// Set whether elastic or inelastic constitutive relations are used.
+void
+pylith::materials::DruckerPragerPlaneStrain::useElasticBehavior(const bool flag)
+{ // useElasticBehavior
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::DruckerPragerPlaneStrain::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastic;
+ _updateStateVarsFn =
+ &pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastic;
+
+ } else {
+ _calcStressFn =
+ &pylith::materials::DruckerPragerPlaneStrain::_calcStressElastoplastic;
+ _calcElasticConstsFn =
+ &pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastoplastic;
+ _updateStateVarsFn =
+ &pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastoplastic;
+ } // if/else
+} // useElasticBehavior
+
+// ----------------------------------------------------------------------
+// Compute properties from values in spatial database.
+void
+pylith::materials::DruckerPragerPlaneStrain::_dbToProperties(
+ PylithScalar* const propValues,
+ const scalar_array& dbValues)
+{ // _dbToProperties
+ assert(0 != propValues);
+ const int numDBValues = dbValues.size();
+ assert(_DruckerPragerPlaneStrain::numDBProperties == numDBValues);
+
+ const PylithScalar density = dbValues[db_density];
+ const PylithScalar vs = dbValues[db_vs];
+ const PylithScalar vp = dbValues[db_vp];
+ const PylithScalar frictionAngle = dbValues[db_frictionAngle];
+ const PylithScalar cohesion = dbValues[db_cohesion];
+ const PylithScalar dilatationAngle = dbValues[db_dilatationAngle];
+
+ const PylithScalar pi = M_PI;
+
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || cohesion < 0.0 ||
+ frictionAngle < 0.0 || frictionAngle > pi/2 ||
+ dilatationAngle < 0.0 || dilatationAngle > pi/2 ||
+ frictionAngle < dilatationAngle) {
+ std::ostringstream msg;
+ msg << "Spatial database returned illegal value for physical "
+ << "properties.\n"
+ << "density: " << density << "\n"
+ << "vp: " << vp << "\n"
+ << "vs: " << vs << "\n"
+ << "frictionAngle: " << frictionAngle << "\n"
+ << "cohesion: " << cohesion << "\n"
+ << "dilatationAngle: " << dilatationAngle << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ if (fabs(frictionAngle - dilatationAngle) > 1.0e-6)
+ _isJacobianSymmetric = false;
+
+ const PylithScalar mu = density * vs*vs;
+ const PylithScalar lambda = density * vp*vp - 2.0*mu;
+
+ PylithScalar alphaYield = 0.0;
+ PylithScalar beta = 0.0;
+ PylithScalar alphaFlow = 0.0;
+ switch (_fitMohrCoulomb) { // switch
+ case MOHR_COULOMB_INSCRIBED: {
+ const PylithScalar denomFriction = sqrt(3.0) * (3.0 - sin(frictionAngle));
+ const PylithScalar denomDilatation = sqrt(3.0) *
+ (3.0 - sin(dilatationAngle));
+ alphaYield = 2.0 * sin(frictionAngle)/denomFriction;
+ beta = 6.0 * cohesion * cos(frictionAngle)/denomFriction;
+ alphaFlow = 2.0 * sin(dilatationAngle)/denomDilatation;
+ break;
+ } // MOHR_COULOMB_INSCRIBED
+ case MOHR_COULOMB_MIDDLE: {
+ alphaYield = sin(frictionAngle)/3.0;
+ beta = cohesion * cos(frictionAngle);
+ alphaFlow = sin(dilatationAngle)/3.0;
+ break;
+ } // MOHR_COULOMB_MIDDLE
+ case MOHR_COULOMB_CIRCUMSCRIBED: {
+ const PylithScalar denomFriction = sqrt(3.0) * (3.0 + sin(frictionAngle));
+ const PylithScalar denomDilatation = sqrt(3.0) *
+ (3.0 + sin(dilatationAngle));
+ alphaYield = 2.0 * sin(frictionAngle)/denomFriction;
+ beta = 6.0 * cohesion * cos(frictionAngle)/denomFriction;
+ alphaFlow = 2.0 * sin(dilatationAngle)/denomDilatation;
+ break;
+ } // MOHR_COULOMB_CIRCUMSCRIBED
+ default :
+ assert(0);
+ throw std::logic_error("Unknown Mohr-Coulomb fit.");
+ break;
+ } // switch
+
+ 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_mu] = mu;
+ propValues[p_lambda] = lambda;
+ propValues[p_alphaYield] = alphaYield;
+ propValues[p_beta] = beta;
+ propValues[p_alphaFlow] = alphaFlow;
+
+ PetscLogFlops(28);
+} // _dbToProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::DruckerPragerPlaneStrain::_nondimProperties(PylithScalar* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(_normalizer);
+ assert(values);
+ assert(nvalues == _numPropsQuadPt);
+
+ const PylithScalar densityScale = _normalizer->densityScale();
+ const PylithScalar pressureScale = _normalizer->pressureScale();
+
+ values[p_density] =
+ _normalizer->nondimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->nondimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->nondimensionalize(values[p_lambda], pressureScale);
+ values[p_beta] =
+ _normalizer->nondimensionalize(values[p_beta], pressureScale);
+
+ PetscLogFlops(4);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::DruckerPragerPlaneStrain::_dimProperties(PylithScalar* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(_normalizer);
+ assert(values);
+ assert(nvalues == _numPropsQuadPt);
+
+ const PylithScalar densityScale = _normalizer->densityScale();
+ const PylithScalar pressureScale = _normalizer->pressureScale();
+
+ values[p_density] =
+ _normalizer->dimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->dimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->dimensionalize(values[p_lambda], pressureScale);
+ values[p_beta] =
+ _normalizer->dimensionalize(values[p_beta], pressureScale);
+
+ PetscLogFlops(4);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Compute initial state variables from values in spatial database.
+void
+pylith::materials::DruckerPragerPlaneStrain::_dbToStateVars(
+ PylithScalar* const stateValues,
+ const scalar_array& dbValues)
+{ // _dbToStateVars
+ assert(stateValues);
+ const int numDBValues = dbValues.size();
+ assert(_DruckerPragerPlaneStrain::numDBStateVars == numDBValues);
+
+ const int totalSize = 5;
+ assert(totalSize == _numVarsQuadPt);
+ assert(totalSize == numDBValues);
+ memcpy(stateValues, &dbValues[0], totalSize*sizeof(PylithScalar));
+
+} // _dbToStateVars
+
+// ----------------------------------------------------------------------
+// Nondimensionalize state variables.
+void
+pylith::materials::DruckerPragerPlaneStrain::_nondimStateVars(
+ PylithScalar* const values,
+ const int nvalues) const
+{ // _nondimStateVars
+ assert(_normalizer);
+ assert(values);
+ assert(nvalues == _numVarsQuadPt);
+
+ const PylithScalar pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(&values[s_stressZZInitial], 1, pressureScale);
+
+ PetscLogFlops(1);
+} // _nondimStateVars
+
+// ----------------------------------------------------------------------
+// Dimensionalize state variables.
+void
+pylith::materials::DruckerPragerPlaneStrain::_dimStateVars(
+ PylithScalar* const values,
+ const int nvalues) const
+{ // _dimStateVars
+ assert(_normalizer);
+ assert(values);
+ assert(nvalues == _numVarsQuadPt);
+
+ const PylithScalar pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(&values[s_stressZZInitial], 1, pressureScale);
+
+ PetscLogFlops(1);
+} // _dimStateVars
+
+// ----------------------------------------------------------------------
+// Compute density at location from properties.
+void
+pylith::materials::DruckerPragerPlaneStrain::_calcDensity(
+ PylithScalar* const density,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars)
+{ // _calcDensity
+ assert(0 != density);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+
+ density[0] = properties[p_density];
+} // _calcDensity
+
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+PylithScalar
+pylith::materials::DruckerPragerPlaneStrain::_stableTimeStepImplicit(
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars) const
+{ // _stableTimeStepImplicit
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ // It's unclear what to do for an elasto-plastic material, which has no
+ // inherent time scale. For now, just set dtStable to a large value.
+ const PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
+
+ return dtStable;
+} // _stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties as an elastic
+// material.
+void
+pylith::materials::DruckerPragerPlaneStrain::_calcStressElastic(
+ PylithScalar* const stress,
+ const int stressSize,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars)
+{ // _calcStressElastic
+ assert(stress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == stressSize);
+ assert(properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(totalStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == strainSize);
+ assert(initialStress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStressSize);
+ assert(initialStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStrainSize);
+
+ const PylithScalar mu = properties[p_mu];
+ const PylithScalar lambda = properties[p_lambda];
+ const PylithScalar mu2 = 2.0 * mu;
+
+ const PylithScalar e11 = totalStrain[0] - initialStrain[0];
+ const PylithScalar e22 = totalStrain[1] - initialStrain[1];
+ const PylithScalar e12 = totalStrain[2] - initialStrain[2];
+
+ const PylithScalar 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 an elastoplastic
+// material.
+void
+pylith::materials::DruckerPragerPlaneStrain::_calcStressElastoplastic(
+ PylithScalar* const stress,
+ const int stressSize,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars)
+{ // _calcStressElastoplastic
+ assert(stress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == stressSize);
+ assert(properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(totalStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == strainSize);
+ assert(initialStress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStressSize);
+ assert(initialStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStrainSize);
+
+ const int tensorSize = 3;
+ const int tensorSizePS = 4;
+ assert(_tensorSize == tensorSize);
+ const PylithScalar mu = properties[p_mu];
+ const PylithScalar lambda = properties[p_lambda];
+
+ // We need to compute the plastic strain increment if state variables are
+ // from previous time step.
+ if (computeStateVars) {
+
+ const PylithScalar alphaYield = properties[p_alphaYield];
+ const PylithScalar beta = properties[p_beta];
+ const PylithScalar alphaFlow = properties[p_alphaFlow];
+ const PylithScalar mu2 = 2.0 * mu;
+ const PylithScalar bulkModulus = lambda + mu2/3.0;
+ const PylithScalar ae = 1.0/mu2;
+ const PylithScalar am = 1.0/(3.0 * bulkModulus);
+
+ const PylithScalar stressZZInitial = stateVars[s_stressZZInitial];
+
+ const PylithScalar plasticStrainT[tensorSizePS] = {
+ stateVars[s_plasticStrain],
+ stateVars[s_plasticStrain + 1],
+ stateVars[s_plasticStrain + 2],
+ stateVars[s_plasticStrain + 3]};
+ const PylithScalar meanPlasticStrainT = (plasticStrainT[0] +
+ plasticStrainT[1] +
+ plasticStrainT[2])/3.0;
+ PylithScalar devPlasticStrainT[tensorSizePS];
+ calcDeviatoric2DPS(devPlasticStrainT, plasticStrainT, meanPlasticStrainT);
+
+ const PylithScalar diag[tensorSizePS] = { 1.0, 1.0, 1.0, 0.0 };
+
+ // Initial stress values
+ const PylithScalar meanStressInitial = (initialStress[0] +
+ initialStress[1] +
+ stressZZInitial)/3.0;
+ const PylithScalar devStressInitial[tensorSizePS] = {
+ initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ stressZZInitial - meanStressInitial,
+ initialStress[2]};
+
+ // Initial strain values
+ const PylithScalar meanStrainInitial = (initialStrain[0] +
+ initialStrain[1])/3.0;
+ const PylithScalar devStrainInitial[tensorSizePS] = {
+ initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ -meanStrainInitial,
+ initialStrain[2]};
+
+ // Values for current time step
+ const PylithScalar meanStrainTpdt = (totalStrain[0] + totalStrain[1])/3.0;
+ const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT -
+ meanStrainInitial;
+
+ const PylithScalar strainPPTpdt[tensorSizePS] =
+ { totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[2] - devPlasticStrainT[3] - devStrainInitial[3] };
+
+ // Compute trial elastic stresses and yield function to see if yield should
+ // occur.
+ const PylithScalar trialDevStress[tensorSizePS] = {
+ strainPPTpdt[0]/ae + devStressInitial[0],
+ strainPPTpdt[1]/ae + devStressInitial[1],
+ strainPPTpdt[2]/ae + devStressInitial[2],
+ strainPPTpdt[3]/ae + devStressInitial[3]
+ };
+ const PylithScalar trialMeanStress = meanStrainPPTpdt/am +
+ meanStressInitial;
+ const PylithScalar stressInvar2 =
+ sqrt(0.5 * scalarProduct2DPS(trialDevStress, trialDevStress));
+ const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
+ stressInvar2 - beta;
+#if 0 // DEBUGGING
+ std::cout << "Function _calcStressElastoPlastic: elastic" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " trialMeanStress: " << trialMeanStress << std::endl;
+ std::cout << " stressInvar2: " << stressInvar2 << std::endl;
+ std::cout << " yieldFunction: " << yieldFunction << std::endl;
+#endif
+ PetscLogFlops(62);
+
+ // If yield function is greater than zero, compute elastoplastic stress.
+ if (yieldFunction >= 0.0) {
+ const PylithScalar devStressInitialProd =
+ scalarProduct2DPS(devStressInitial, devStressInitial);
+ const PylithScalar strainPPTpdtProd =
+ scalarProduct2DPS(strainPPTpdt, strainPPTpdt);
+ const PylithScalar d =
+ sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
+ scalarProduct2DPS(devStressInitial, strainPPTpdt) +
+ strainPPTpdtProd);
+ const PylithScalar plasticMult = 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ const PylithScalar meanStressTpdt =
+ (meanStrainPPTpdt - plasticMult * alphaFlow)/am + meanStressInitial;
+ PylithScalar deltaDevPlasticStrain = 0.0;
+ PylithScalar devStressTpdt = 0.0;
+ PylithScalar totalStress[tensorSizePS];
+ for (int iComp=0; iComp < tensorSizePS; ++iComp) {
+ deltaDevPlasticStrain = plasticMult *(strainPPTpdt[iComp] +
+ ae * devStressInitial[iComp])/
+ (sqrt(2.0) * d);
+ devStressTpdt = (strainPPTpdt[iComp] - deltaDevPlasticStrain)/ae +
+ devStressInitial[iComp];
+ totalStress[iComp] = devStressTpdt + diag[iComp] * meanStressTpdt;
+ } // for
+ stress[0] = totalStress[0];
+ stress[1] = totalStress[1];
+ stress[2] = totalStress[3];
+
+ PetscLogFlops(51 + 11 * tensorSizePS);
+
+ } else {
+ // No plastic strain.
+ const PylithScalar meanStressTpdt = meanStrainPPTpdt/am +
+ meanStressInitial;
+ stress[0] = strainPPTpdt[0]/ae + devStressInitial[0] + meanStressTpdt;
+ stress[1] = strainPPTpdt[1]/ae + devStressInitial[1] + meanStressTpdt;
+ stress[2] = strainPPTpdt[3]/ae + devStressInitial[3];
+
+ PetscLogFlops(10);
+ } // if
+
+ // If state variables have already been updated, the plastic strain for the
+ // time step has already been computed.
+ } else {
+ const PylithScalar mu2 = 2.0 * mu;
+ const PylithScalar plasticStrainTpdt[tensorSize] = {
+ stateVars[s_plasticStrain],
+ stateVars[s_plasticStrain + 1],
+ stateVars[s_plasticStrain + 3]
+ };
+
+ const PylithScalar e11 = totalStrain[0] - plasticStrainTpdt[0] -
+ initialStrain[0];
+ const PylithScalar e22 = totalStrain[1] - plasticStrainTpdt[1] -
+ initialStrain[1];
+ const PylithScalar e12 = totalStrain[2] - plasticStrainTpdt[2] -
+ initialStrain[2];
+
+ const PylithScalar traceStrainTpdt = e11 + e22;
+ const PylithScalar s12 = lambda * traceStrainTpdt;
+
+ stress[0] = s12 + mu2 * e11 + initialStress[0];
+ stress[1] = s12 + mu2 * e22 + initialStress[1];
+ stress[2] = mu2 * e12 + initialStress[2];
+
+ PetscLogFlops(17);
+
+ } // else
+#if 0 // DEBUGGING
+ //** This debugging code won't work for plane strain since we're missing
+ // stressZZ. Not sure if this can be fixed.
+ const PylithScalar alphaYield = properties[p_alphaYield];
+ const PylithScalar beta = properties[p_beta];
+ const PylithScalar alphaFlow = properties[p_alphaFlow];
+ const PylithScalar meanStressTest = (stress[0] + stress[1])/3.0;
+ const PylithScalar devStressTest[] = { stress[0] - meanStressTest,
+ stress[1] - meanStressTest,
+ stress[2] - meanStressTest,
+ stress[3],
+ stress[4],
+ stress[5]};
+ const PylithScalar stressInvar2Test =
+ sqrt(0.5 *
+ pylith::materials::ElasticMaterial::scalarProduct3D(devStressTest,
+ devStressTest));
+
+ const PylithScalar yieldFunctionTest = 3.0 * alphaYield * meanStressTest +
+ stressInvar2Test - beta;
+ std::cout << "Function _calcStressElastoPlastic: end" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " meanStressTest: " << meanStressTest << std::endl;
+ std::cout << " stressInvar2Test: " << stressInvar2Test << std::endl;
+ std::cout << " yieldFunctionTest: " << yieldFunctionTest << std::endl;
+#endif
+
+} // _calcStressElastoplastic
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from properties.
+void
+pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastic(
+ PylithScalar* const elasticConsts,
+ const int numElasticConsts,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize)
+{ // _calcElasticConstsElastic
+ assert(elasticConsts);
+ assert(_DruckerPragerPlaneStrain::numElasticConsts == numElasticConsts);
+ assert(properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(totalStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == strainSize);
+ assert(initialStress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStressSize);
+ assert(initialStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStrainSize);
+
+ const PylithScalar mu = properties[p_mu];
+ const PylithScalar lambda = properties[p_lambda];
+
+ const PylithScalar mu2 = 2.0 * mu;
+ const PylithScalar 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 an elastoplastic material.
+void
+pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastoplastic(
+ PylithScalar* const elasticConsts,
+ const int numElasticConsts,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize)
+{ // _calcElasticConstsElastoplastic
+ assert(elasticConsts);
+ assert(_DruckerPragerPlaneStrain::numElasticConsts == numElasticConsts);
+ assert(properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(totalStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == strainSize);
+ assert(initialStress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStressSize);
+ assert(initialStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStrainSize);
+
+ // Duplicate functionality of _calcStressElastoplastic
+ // Get properties
+ const int tensorSize = 3;
+ const int tensorSizePS = 4;
+ const PylithScalar mu = properties[p_mu];
+ const PylithScalar lambda = properties[p_lambda];
+ const PylithScalar alphaYield = properties[p_alphaYield];
+ const PylithScalar beta = properties[p_beta];
+ const PylithScalar alphaFlow = properties[p_alphaFlow];
+
+ const PylithScalar mu2 = 2.0 * mu;
+ const PylithScalar bulkModulus = lambda + mu2/3.0;
+ const PylithScalar ae = 1.0/mu2;
+ const PylithScalar am = 1.0/(3.0 * bulkModulus);
+
+ const PylithScalar stressZZInitial = stateVars[s_stressZZInitial];
+
+ // Get state variables from previous time step
+ const PylithScalar plasticStrainT[tensorSizePS] = {
+ stateVars[s_plasticStrain],
+ stateVars[s_plasticStrain + 1],
+ stateVars[s_plasticStrain + 2],
+ stateVars[s_plasticStrain + 3]};
+ const PylithScalar meanPlasticStrainT = (plasticStrainT[0] +
+ plasticStrainT[1] +
+ plasticStrainT[2])/3.0;
+ PylithScalar devPlasticStrainT[tensorSizePS];
+ calcDeviatoric2DPS(devPlasticStrainT, plasticStrainT, meanPlasticStrainT);
+
+ const PylithScalar diag[tensorSize] = { 1.0, 1.0, 0.0 };
+
+ // Initial stress values
+ const PylithScalar meanStressInitial = (initialStress[0] +
+ initialStress[1] +
+ stressZZInitial)/3.0;
+ const PylithScalar devStressInitial[tensorSizePS] = {
+ initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ stressZZInitial - meanStressInitial,
+ initialStress[2]};
+
+ // Initial strain values
+ const PylithScalar meanStrainInitial = (initialStrain[0] +
+ initialStrain[1])/3.0;
+ const PylithScalar devStrainInitial[tensorSizePS] = {
+ initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ -meanStrainInitial,
+ initialStrain[2]};
+
+ // Values for current time step
+ const PylithScalar meanStrainTpdt = (totalStrain[0] + totalStrain[1])/3.0;
+ const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT -
+ meanStrainInitial;
+
+ const PylithScalar strainPPTpdt[tensorSizePS] =
+ { totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ -meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[2] - devPlasticStrainT[3] - devStrainInitial[3]};
+
+ // Compute trial elastic stresses and yield function to see if yield should
+ // occur.
+ const PylithScalar trialDevStress[tensorSizePS] = {
+ strainPPTpdt[0]/ae + devStressInitial[0],
+ strainPPTpdt[1]/ae + devStressInitial[1],
+ strainPPTpdt[2]/ae + devStressInitial[2],
+ strainPPTpdt[3]/ae + devStressInitial[3]};
+ const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
+ const PylithScalar stressInvar2 =
+ sqrt(0.5 * scalarProduct2DPS(trialDevStress, trialDevStress));
+ const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
+ stressInvar2 - beta;
+#if 0 // DEBUGGING
+ std::cout << "Function _calcElasticConstsElastoPlastic:" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " trialMeanStress: " << trialMeanStress << std::endl;
+ std::cout << " stressInvar2: " << stressInvar2 << std::endl;
+ std::cout << " yieldFunction: " << yieldFunction << std::endl;
+#endif
+ PetscLogFlops(62);
+
+ // If yield function is greater than zero, compute elastoplastic stress and
+ // corresponding tangent matrix.
+ if (yieldFunction >= 0.0) {
+ const PylithScalar devStressInitialProd =
+ scalarProduct2DPS(devStressInitial, devStressInitial);
+ const PylithScalar strainPPTpdtProd =
+ scalarProduct2DPS(strainPPTpdt, strainPPTpdt);
+ const PylithScalar d =
+ sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
+ scalarProduct2DPS(devStressInitial, strainPPTpdt) +
+ strainPPTpdtProd);
+ const PylithScalar plasticFac = 2.0 * ae * am/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
+ const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
+ const PylithScalar plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+
+ // Define some constants, vectors, and matrices.
+ const PylithScalar third = 1.0/3.0;
+ const PylithScalar dEdEpsilon[3][3] = {
+ { 2.0 * third, -third, 0.0},
+ { -third, 2.0 * third, 0.0},
+ { 0.0, 0.0, 1.0}};
+ const PylithScalar vec1[3] = {strainPPTpdt[0] + ae * devStressInitial[0],
+ strainPPTpdt[1] + ae * devStressInitial[1],
+ strainPPTpdt[3] + ae * devStressInitial[3]};
+ const PylithScalar dDdEpsilon[3] = {vec1[0]/d,
+ vec1[1]/d,
+ 2.0 * vec1[2]/d};
+ const PylithScalar dLambdadEpsilon[3] = {
+ plasticFac * (alphaYield/am + dFac * dDdEpsilon[0]),
+ plasticFac * (alphaYield/am + dFac * dDdEpsilon[1]),
+ plasticFac * ( dFac * dDdEpsilon[2])};
+
+ const PylithScalar dFac2 = 1.0/(sqrt(2.0) * d);
+ PylithScalar dDeltaEdEpsilon = 0.0;
+
+ // Compute elasticity matrix.
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ for (int jComp=0; jComp < tensorSize; ++jComp) {
+ int iCount = jComp + tensorSize * iComp;
+ dDeltaEdEpsilon = dFac2 * (vec1[iComp] *
+ (dLambdadEpsilon[jComp] -
+ plasticMult * dDdEpsilon[jComp]/d) +
+ plasticMult * dEdEpsilon[iComp][jComp]);
+ elasticConsts[iCount] = (dEdEpsilon[iComp][jComp] -
+ dDeltaEdEpsilon)/ae +
+ diag[iComp] * (third * diag[jComp] -
+ alphaFlow * dLambdadEpsilon[jComp])/am;
+ } // for
+ } // for
+
+ PetscLogFlops(76 + tensorSize * tensorSize * 15);
+
+ } else {
+ // No plastic strain.
+ const PylithScalar 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(1);
+ } // else
+
+} // _calcElasticConstsElastoplastic
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastic(
+ PylithScalar* const stateVars,
+ const int numStateVars,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVarsElastic
+ assert(stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(totalStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == strainSize);
+ assert(initialStress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStressSize);
+ assert(initialStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStrainSize);
+
+ for (int iComp=0; iComp < 4; ++iComp) {
+ stateVars[s_plasticStrain+iComp] = 0.0;
+ } // for
+
+ _needNewJacobian = true;
+} // _updateStateVarsElastic
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastoplastic(
+ PylithScalar* const stateVars,
+ const int numStateVars,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVarsElastoplastic
+ assert(stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(totalStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == strainSize);
+ assert(initialStress);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStressSize);
+ assert(initialStrain);
+ assert(_DruckerPragerPlaneStrain::tensorSize == initialStrainSize);
+
+ const int stressSize = _tensorSize;
+
+ // For now, we are duplicating the functionality of _calcStressElastoplastic,
+ // since otherwise we would have to redo a lot of calculations.
+
+ const int tensorSize = _tensorSize;
+ const int tensorSizePS = 4;
+ const PylithScalar mu = properties[p_mu];
+ const PylithScalar lambda = properties[p_lambda];
+ const PylithScalar alphaYield = properties[p_alphaYield];
+ const PylithScalar beta = properties[p_beta];
+ const PylithScalar alphaFlow = properties[p_alphaFlow];
+
+ const PylithScalar mu2 = 2.0 * mu;
+ const PylithScalar bulkModulus = lambda + mu2/3.0;
+ const PylithScalar ae = 1.0/mu2;
+ const PylithScalar am = 1.0/(3.0 * bulkModulus);
+
+ const PylithScalar stressZZInitial = stateVars[s_stressZZInitial];
+
+ const PylithScalar plasticStrainT[tensorSizePS] = {
+ stateVars[s_plasticStrain],
+ stateVars[s_plasticStrain + 1],
+ stateVars[s_plasticStrain + 2],
+ stateVars[s_plasticStrain + 3]};
+ const PylithScalar meanPlasticStrainT = (plasticStrainT[0] +
+ plasticStrainT[1] +
+ plasticStrainT[2])/3.0;
+ PylithScalar devPlasticStrainT[tensorSizePS];
+ calcDeviatoric2DPS(devPlasticStrainT, plasticStrainT, meanPlasticStrainT);
+
+ const PylithScalar diag[] = { 1.0, 1.0, 1.0, 0.0 };
+
+ // Initial stress values
+ const PylithScalar meanStressInitial = (initialStress[0] +
+ initialStress[1] +
+ stressZZInitial)/3.0;
+ const PylithScalar devStressInitial[tensorSizePS] = {
+ initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ stressZZInitial - meanStressInitial,
+ initialStress[2]};
+
+ // Initial strain values
+ const PylithScalar meanStrainInitial = (initialStrain[0] +
+ initialStrain[1])/3.0;
+ const PylithScalar devStrainInitial[tensorSizePS] = {
+ initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ -meanStrainInitial,
+ initialStrain[2]};
+
+ // Values for current time step
+ const PylithScalar meanStrainTpdt = (totalStrain[0] + totalStrain[1])/3.0;
+ const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT -
+ meanStrainInitial;
+
+ const PylithScalar strainPPTpdt[tensorSizePS] =
+ { totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ -meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[2] - devPlasticStrainT[3] - devStrainInitial[3]};
+
+ // Compute trial elastic stresses and yield function to see if yield should
+ // occur.
+ const PylithScalar trialDevStress[tensorSizePS] = {
+ strainPPTpdt[0]/ae + devStressInitial[0],
+ strainPPTpdt[1]/ae + devStressInitial[1],
+ strainPPTpdt[2]/ae + devStressInitial[2],
+ strainPPTpdt[3]/ae + devStressInitial[3]};
+ const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
+ const PylithScalar stressInvar2 =
+ sqrt(0.5 * scalarProduct2DPS(trialDevStress, trialDevStress));
+ const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
+ stressInvar2 - beta;
+#if 0 // DEBUGGING
+ std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " trialMeanStress: " << trialMeanStress << std::endl;
+ std::cout << " stressInvar2: " << stressInvar2 << std::endl;
+ std::cout << " yieldFunction: " << yieldFunction << std::endl;
+#endif
+ PetscLogFlops(62);
+
+ // If yield function is greater than zero, compute plastic strains.
+ // Otherwise, plastic strains remain the same.
+ if (yieldFunction >= 0.0) {
+ const PylithScalar devStressInitialProd =
+ scalarProduct2DPS(devStressInitial, devStressInitial);
+ const PylithScalar strainPPTpdtProd =
+ scalarProduct2DPS(strainPPTpdt, strainPPTpdt);
+ const PylithScalar d =
+ sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
+ scalarProduct2DPS(devStressInitial, strainPPTpdt) +
+ strainPPTpdtProd);
+ const PylithScalar plasticMult = 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
+ PylithScalar deltaDevPlasticStrain = 0.0;
+ for (int iComp=0; iComp < tensorSizePS; ++iComp) {
+ deltaDevPlasticStrain = plasticMult *(strainPPTpdt[iComp] +
+ ae * devStressInitial[iComp])/
+ (sqrt(2.0) * d);
+ stateVars[s_plasticStrain+iComp] += deltaDevPlasticStrain +
+ diag[iComp] * deltaMeanPlasticStrain;
+ } // for
+
+ PetscLogFlops(48 + 9 * tensorSizePS);
+
+ } // if
+
+ _needNewJacobian = true;
+
+} // _updateStateVarsElastoplastic
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh (from rev 19598, short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,533 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/DruckerPragerPlaneStrain.hh
+ *
+ * @brief 2-D, plane strain, Drucker-Prager elastic/perfectly plastic material.
+ */
+
+#if !defined(pylith_materials_druckerpragerplanestrain_hh)
+#define pylith_materials_druckerpragerplanestrain_hh
+
+// Include directives ---------------------------------------------------
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
+
+// DruckerPragerPlaneStrain ---------------------------------------------
+/** @brief 2-D, plane strain, Drucker-Prager elastic/perfectly plastic material.
+ *
+ * The physical properties are specified using density, shear-wave
+ * speed, friction angle, cohesion, dilatation angle, and
+ * compressional-wave speed. The physical properties are stored
+ * internally using density, lambda, mu, which are directly related to
+ * the elasticity constants used in the finite-element
+ * integration. The plasticity information is retained as alpha_yield,
+ * beta, and alpha_flow.
+ */
+
+class pylith::materials::DruckerPragerPlaneStrain : public ElasticMaterial
+{ // class DruckerPragerPlaneStrain
+ friend class TestDruckerPragerPlaneStrain; // unit testing
+
+ // PUBLIC ENUMS ///////////////////////////////////////////////////////
+public :
+
+ enum FitMohrCoulombEnum {
+ MOHR_COULOMB_CIRCUMSCRIBED=0,
+ MOHR_COULOMB_MIDDLE=1,
+ MOHR_COULOMB_INSCRIBED=2,
+ }; // FitMohrCoulombType
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor
+ DruckerPragerPlaneStrain(void);
+
+ /// Destructor
+ ~DruckerPragerPlaneStrain(void);
+
+ /** Set fit to Mohr-Coulomb surface.
+ *
+ * @param value Mohr-Coulomb surface match type.
+ */
+ void fitMohrCoulomb(FitMohrCoulombEnum value);
+
+ /** Set current time step.
+ *
+ * @param dt Current time step.
+ */
+ void timeStep(const PylithScalar 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(PylithScalar* const propValues,
+ const scalar_array& dbValues);
+
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _nondimProperties(PylithScalar* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _dimProperties(PylithScalar* 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(PylithScalar* const stateValues,
+ const scalar_array& dbValues);
+
+ /** Nondimensionalize state variables..
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _nondimStateVars(PylithScalar* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize state variables.
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _dimStateVars(PylithScalar* 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 State variables at location.
+ * @param numStateVars Number of state variables.
+ */
+ void _calcDensity(PylithScalar* const density,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* 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(PylithScalar* const stress,
+ const int stressSize,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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(PylithScalar* const elasticConsts,
+ const int numElasticConsts,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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
+ */
+ PylithScalar _stableTimeStepImplicit(const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars) const;
+
+ /** 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(PylithScalar* const stateVars,
+ const int numStateVars,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize);
+
+ // PRIVATE TYPEDEFS ///////////////////////////////////////////////////
+private :
+
+ /// Member prototype for _calcStress()
+ typedef void (pylith::materials::DruckerPragerPlaneStrain::*calcStress_fn_type)
+ (PylithScalar* const,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const bool);
+
+ /// Member prototype for _calcElasticConsts()
+ typedef void (pylith::materials::DruckerPragerPlaneStrain::*calcElasticConsts_fn_type)
+ (PylithScalar* const,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int);
+
+ /// Member prototype for _updateStateVars()
+ typedef void (pylith::materials::DruckerPragerPlaneStrain::*updateStateVars_fn_type)
+ (PylithScalar* const,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ const int,
+ const PylithScalar*,
+ 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 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 vars.
+ */
+ void _calcStressElastic(PylithScalar* const stress,
+ const int stressSize,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars);
+
+ /** Compute stress tensor from properties as an elastoplastic 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 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 vars.
+ */
+ void _calcStressElastoplastic(PylithScalar* const stress,
+ const int stressSize,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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 locations.
+ * @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 _calcElasticConstsElastic(PylithScalar* const elasticConsts,
+ const int numElasticConsts,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize);
+
+ /** Compute derivatives of elasticity matrix from properties as an
+ * elastoplastic 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 locations.
+ * @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 _calcElasticConstsElastoplastic(PylithScalar* const elasticConsts,
+ const int numElasticConsts,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize);
+
+ /** Update state variables after solve as an elastic material.
+ *
+ * @param stateVars State variables at locations.
+ * @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 _updateStateVarsElastic(PylithScalar* const stateVars,
+ const int numStateVars,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize);
+
+ /** Update state variables after solve as an elastoplastic material.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialState Initial state values.
+ * @param initialStateSize Size of initial state array.
+ */
+ void _updateStateVarsElastoplastic(PylithScalar* const stateVars,
+ const int numStateVars,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize);
+
+ /** Compute scalar product, assuming vector form of a tensor.
+ *
+ * @param tensor1 First tensor.
+ * @param tensor2 Second tensor.
+ */
+ /*
+ PylithScalar _scalarProduct(const PylithScalar* tensor1,
+ const PylithScalar* tensor2) const;
+ */
+
+ /** Compute tensor mean, assuming vector form of a tensor.
+ *
+ * @param vec Tensor represented as a vector.
+ */
+ PylithScalar _calcMean(const PylithScalar* vec) const;
+
+ /** Compute deviatoric components, assuming vector form of a tensor.
+ *
+ * @param vec Tensor represented as a vector.
+ * @param vecMean Tensor mean.
+ */
+ PylithScalar _calcDeviatoric(const PylithScalar* vec,
+ const PylithScalar vecMean) const;
+
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ /// 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;
+
+ /// Fit to Mohr Coulomb surface
+ FitMohrCoulombEnum _fitMohrCoulomb;
+
+ static const int p_density;
+ static const int p_mu;
+ static const int p_lambda;
+ static const int p_alphaYield;
+ static const int p_beta;
+ static const int p_alphaFlow;
+ static const int db_density;
+ static const int db_vs;
+ static const int db_vp;
+ static const int db_frictionAngle;
+ static const int db_cohesion;
+ static const int db_dilatationAngle;
+
+ static const int s_stressZZInitial;
+ static const int s_plasticStrain;
+ static const int db_stressZZInitial;
+ static const int db_plasticStrain;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ DruckerPragerPlaneStrain(const DruckerPragerPlaneStrain&); ///< Not implemented
+ const DruckerPragerPlaneStrain& operator=(const DruckerPragerPlaneStrain&); ///< Not implemented
+
+}; // class DruckerPragerPlaneStrain
+
+#include "DruckerPragerPlaneStrain.icc" // inline methods
+
+#endif // pylith_materials_druckerpragerplanestrain_hh
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc (from rev 19598, short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,112 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_materials_druckerpragerplanestrain_hh)
+#error "DruckerPragerPlaneStrain.icc can only be included from DruckerPragerPlaneStrain.hh"
+#endif
+
+#include <cassert> // USES assert()
+#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+
+// Set current time step.
+inline
+void
+pylith::materials::DruckerPragerPlaneStrain::timeStep(const PylithScalar dt) {
+ // Not sure what to do here. If we are using full Newton the Jacobian will
+ // always need reforming, but SNES may opt not to reform it sometimes.
+ _needNewJacobian = true;
+ _dt = dt;
+} // timeStep
+
+// Compute stress tensor from parameters.
+inline
+void
+pylith::materials::DruckerPragerPlaneStrain::_calcStress(PylithScalar* const stress,
+ const int stressSize,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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::DruckerPragerPlaneStrain::_calcElasticConsts(
+ PylithScalar* const elasticConsts,
+ const int numElasticConsts,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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::DruckerPragerPlaneStrain::_updateStateVars(
+ PylithScalar* const stateVars,
+ const int numStateVars,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.hh 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.hh 2012-02-07 21:21:03 UTC (rev 19599)
@@ -329,6 +329,17 @@
void calcDeviatoric2D(PylithScalar* const deviatoric,
const PylithScalar* vec,
const PylithScalar vecMean);
+
+ /** Compute 2D deviatoric stress/strain from vector (length 4) and mean value.
+ *
+ * @param deviatoric Array for deviatoric tensor.
+ * @param vec Input tensor (as vector).
+ * @param vecMean Tensor trace divided by spatial_dimension.
+ */
+ static
+ void calcDeviatoric2DPS(PylithScalar* const deviatoric,
+ const PylithScalar* vec,
+ const PylithScalar vecMean);
/** Compute 3D deviatoric stress/strain from vector and mean value.
*
@@ -350,6 +361,16 @@
PylithScalar scalarProduct2D(const PylithScalar* tensor1,
const PylithScalar* tensor2);
+ /** Compute 2D scalar product of two tensors represented as vectors of
+ * length 4.
+ *
+ * @param tensor1 First tensor.
+ * @param tensor2 Second tensor.
+ */
+ static
+ PylithScalar scalarProduct2DPS(const PylithScalar* tensor1,
+ const PylithScalar* tensor2);
+
/** Compute 3D scalar product of two tensors represented as vectors.
*
* @param tensor1 First tensor.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.icc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.icc 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/ElasticMaterial.icc 2012-02-07 21:21:03 UTC (rev 19599)
@@ -71,6 +71,22 @@
deviatoric[2] = vec[2];
} // calcDeviatoric2D
+// Compute 2D deviatoric stress/strain from vector and mean value.
+// In this case all 3 normal components are present, but only one shear comp.
+// 3 FLOPs per call
+inline
+void
+pylith::materials::ElasticMaterial::calcDeviatoric2DPS(
+ PylithScalar* const deviatoric,
+ const PylithScalar* vec,
+ const PylithScalar vecMean)
+{
+ deviatoric[0] = vec[0] - vecMean;
+ deviatoric[1] = vec[1] - vecMean;
+ deviatoric[2] = vec[2] - vecMean;
+ deviatoric[3] = vec[3];
+} // calcDeviatoric2D
+
// Compute 3D deviatoric stress/strain from vector and mean value.
// 3 FLOPs per call
inline
@@ -105,6 +121,25 @@
} // scalarProduct2D
+// Compute 2DPS scalar product of two tensors represented as vectors.
+// In this case all 3 normal components are present, but only one shear comp.
+// 8 FLOPs per call.
+inline
+PylithScalar
+pylith::materials::ElasticMaterial::scalarProduct2DPS(
+ const PylithScalar* tensor1,
+ const PylithScalar* tensor2)
+{ // scalarProduct2DPS
+ const PylithScalar scalarProduct =
+ tensor1[0] * tensor2[0] +
+ tensor1[1] * tensor2[1] +
+ tensor1[2] * tensor2[2] +
+ 2.0 * tensor1[3] * tensor2[3];
+
+ return scalarProduct;
+} // scalarProduct2DPS
+
+
// Compute 3D scalar product of two tensors represented as vectors.
// 12 FLOPs per call.
inline
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/Makefile.am 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/Makefile.am 2012-02-07 21:21:03 UTC (rev 19599)
@@ -42,6 +42,8 @@
PowerLaw3D.icc \
DruckerPrager3D.hh \
DruckerPrager3D.icc \
+ DruckerPragerPlaneStrain.hh \
+ DruckerPragerPlaneStrain.icc \
Metadata.hh \
Metadata.icc \
Material.hh \
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/materialsfwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/materialsfwd.hh 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/materialsfwd.hh 2012-02-07 21:21:03 UTC (rev 19599)
@@ -48,6 +48,7 @@
class GenMaxwellQpQsIsotropic3D;
class PowerLaw3D;
class DruckerPrager3D;
+ class DruckerPragerPlaneStrain;
class EffectiveStress;
class ViscoelasticMaxwell;
Copied: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i (from rev 19598, short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,239 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file modulesrc/materials/DruckerPragerPlaneStrain.i
+ *
+ * @brief Python interface to C++ DruckerPragerPlaneStrain object.
+ */
+
+namespace pylith {
+ namespace materials {
+
+ class pylith::materials::DruckerPragerPlaneStrain : public ElasticMaterial
+ { // class DruckerPragerPlaneStrain
+
+ // PUBLIC ENUMS ///////////////////////////////////////////////////
+ public :
+
+ enum FitMohrCoulombEnum {
+ MOHR_COULOMB_CIRCUMSCRIBED=0,
+ MOHR_COULOMB_MIDDLE=1,
+ MOHR_COULOMB_INSCRIBED=2,
+ }; // FitMohrCoulombType
+
+ // PUBLIC METHODS /////////////////////////////////////////////////
+ public :
+
+ /// Default constructor
+ DruckerPragerPlaneStrain(void);
+
+ /// Destructor
+ ~DruckerPragerPlaneStrain(void);
+
+ /** Set fit to Mohr-Coulomb surface.
+ *
+ * @param value Mohr-Coulomb surface match type.
+ */
+ void fitMohrCoulomb(FitMohrCoulombEnum value);
+
+ /** Set current time step.
+ *
+ * @param dt Current time step.
+ */
+ void timeStep(const PylithScalar 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(PylithScalar* const propValues,
+ const pylith::scalar_array& dbValues);
+
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _nondimProperties(PylithScalar* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _dimProperties(PylithScalar* 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(PylithScalar* const stateValues,
+ const pylith::scalar_array& dbValues);
+
+ /** Nondimensionalize state variables..
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _nondimStateVars(PylithScalar* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize state variables.
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _dimStateVars(PylithScalar* 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 State variables at location.
+ * @param numStateVars Number of state variables.
+ */
+ void _calcDensity(PylithScalar* const density,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* 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(PylithScalar* const stress,
+ const int stressSize,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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(PylithScalar* const elasticConsts,
+ const int numElasticConsts,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* 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
+ */
+ PylithScalar _stableTimeStepImplicit(const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars) const;
+
+ /** 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(PylithScalar* const stateVars,
+ const int numStateVars,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* totalStrain,
+ const int strainSize,
+ const PylithScalar* initialStress,
+ const int initialStressSize,
+ const PylithScalar* initialStrain,
+ const int initialStrainSize);
+
+ }; // class DruckerPragerPlaneStrain
+
+ } // materials
+} // pylith
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/Makefile.am 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/Makefile.am 2012-02-07 21:21:03 UTC (rev 19599)
@@ -38,7 +38,8 @@
GenMaxwellPlaneStrain.i \
GenMaxwellQpQsIsotropic3D.i \
PowerLaw3D.i \
- DruckerPrager3D.i
+ DruckerPrager3D.i \
+ DruckerPragerPlaneStrain.i
swigincludedir = $(pkgdatadir)/swig/$(subpackage)
swiginclude_HEADERS = \
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/materials.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/materials.i 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/materials.i 2012-02-07 21:21:03 UTC (rev 19599)
@@ -39,6 +39,7 @@
#include "pylith/materials/GenMaxwellQpQsIsotropic3D.hh"
#include "pylith/materials/PowerLaw3D.hh"
#include "pylith/materials/DruckerPrager3D.hh"
+#include "pylith/materials/DruckerPragerPlaneStrain.hh"
#include "pylith/utils/arrayfwd.hh"
#include "pylith/utils/sievetypes.hh"
@@ -82,6 +83,7 @@
%include "GenMaxwellQpQsIsotropic3D.i"
%include "PowerLaw3D.i"
%include "DruckerPrager3D.i"
+%include "DruckerPragerPlaneStrain.i"
// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am 2012-02-07 21:21:03 UTC (rev 19599)
@@ -86,6 +86,7 @@
materials/MaxwellPlaneStrain.py \
materials/PowerLaw3D.py \
materials/DruckerPrager3D.py \
+ materials/DruckerPragerPlaneStrain.py \
meshio/__init__.py \
meshio/CellFilter.py \
meshio/CellFilterAvgMesh.py \
Copied: short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPragerPlaneStrain.py (from rev 19598, short/3D/PyLith/trunk/pylith/materials/DruckerPragerPlaneStrain.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPragerPlaneStrain.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPragerPlaneStrain.py 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,130 @@
+#!/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-2011 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/materials/DruckerPragerPlaneStrain.py
+##
+## @brief Python object implementing 2-D plane strain Drucker-Prager
+## elastoplastic material.
+##
+## Factory: material.
+
+from ElasticMaterial import ElasticMaterial
+from materials import DruckerPragerPlaneStrain as ModuleDruckerPragerPlaneStrain
+
+# Validator to fit to Mohr-Coulomb
+def validateFitMohrCoulomb(value):
+ """
+ Validate fit to Mohr-Coulomb yield surface.
+ """
+ if not value in ["inscribed", "middle", "circumscribed"]:
+ raise ValueError("Unknown fit to Mohr-Coulomb yield surface.")
+ return value
+
+
+# DruckerPragerPlaneStrain class
+class DruckerPragerPlaneStrain(ElasticMaterial, ModuleDruckerPragerPlaneStrain):
+ """
+ Python object implementing 2-D plane strain Drucker-Prager elastoplastic
+ material.
+
+ Factory: material.
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(ElasticMaterial.Inventory):
+ """
+ Python object for managing DruckerPragerPlaneStrain facilities and
+ properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing DruckerPragerPlaneStrain facilities and
+ ## properties.
+ ##
+ ## \b Properties
+ ## @li \b fit_mohr_coulomb Fit to Mohr-Coulomb yield surface.
+ ##
+ ## \b Facilities
+ ## @li None
+
+ import pyre.inventory
+
+ from pylith.meshio.OutputMatElastic import OutputMatElastic
+ fitMohrCoulomb = pyre.inventory.str("fit_mohr_coulomb", default="inscribed",
+ validator=validateFitMohrCoulomb)
+ fitMohrCoulomb.meta['tip'] = "Fit to Mohr-Coulomb yield surface."
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="druckerpragerplanestrain"):
+ """
+ Constructor.
+ """
+ ElasticMaterial.__init__(self, name)
+ self.availableFields = \
+ {'vertex': \
+ {'info': [],
+ 'data': []},
+ 'cell': \
+ {'info': ["mu", "lambda", "density",
+ "alpha_yield", "beta", "alpha_flow"],
+ 'data': ["total_strain", "stress", "stress_zz_initial",
+ "plastic_strain"]}}
+ self._loggingPrefix = "MaDP2D "
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Setup members using inventory.
+ """
+ ElasticMaterial._configure(self)
+ if self.inventory.fitMohrCoulomb == "inscribed":
+ fitEnum = ModuleDruckerPragerPlaneStrain.MOHR_COULOMB_INSCRIBED
+ elif self.inventory.fitMohrCoulomb == "middle":
+ fitEnum = ModuleDruckerPragerPlaneStrain.MOHR_COULOMB_MIDDLE
+ elif self.inventory.fitMohrCoulomb == "circumscribed":
+ fitEnum = ModuleDruckerPragerPlaneStrain.MOHR_COULOMB_CIRCUMSCRIBED
+ else:
+ raise ValueError("Unknown fit to Mohr-Coulomb yield surface.")
+ ModuleDruckerPragerPlaneStrain.fitMohrCoulomb(self, fitEnum)
+ return
+
+
+ def _createModuleObj(self):
+ """
+ Call constructor for module object for access to C++ object.
+ """
+ ModuleDruckerPragerPlaneStrain.__init__(self)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def material():
+ """
+ Factory associated with DruckerPragerPlaneStrain.
+ """
+ return DruckerPragerPlaneStrain()
+
+
+# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/Makefile.am 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/Makefile.am 2012-02-07 21:21:03 UTC (rev 19599)
@@ -42,6 +42,7 @@
TestGenMaxwellQpQsIsotropic3D.cc \
TestPowerLaw3D.cc \
TestDruckerPrager3D.cc \
+ TestDruckerPragerPlaneStrain.cc \
TestEffectiveStress.cc \
test_materials.cc
@@ -62,6 +63,7 @@
TestMaxwellPlaneStrain.hh \
TestPowerLaw3D.hh \
TestDruckerPrager3D.hh \
+ TestDruckerPragerPlaneStrain.hh \
TestEffectiveStress.hh
# Source files associated with testing data
@@ -86,7 +88,9 @@
data/PowerLaw3DElasticData.cc \
data/PowerLaw3DTimeDepData.cc \
data/DruckerPrager3DElasticData.cc \
- data/DruckerPrager3DTimeDepData.cc
+ data/DruckerPrager3DTimeDepData.cc \
+ data/DruckerPragerPlaneStrainElasticData.cc \
+ data/DruckerPragerPlaneStrainTimeDepData.cc
noinst_HEADERS += \
@@ -111,6 +115,8 @@
data/PowerLaw3DTimeDepData.hh \
data/DruckerPrager3DElasticData.hh \
data/DruckerPrager3DTimeDepData.hh \
+ data/DruckerPragerPlaneStrainElasticData.hh \
+ data/DruckerPragerPlaneStrainTimeDepData.hh \
data/header.hh
AM_CPPFLAGS = \
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,240 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestDruckerPragerPlaneStrain.hh" // Implementation of class methods
+
+#include "data/DruckerPragerPlaneStrainElasticData.hh" // USES DruckerPragerPlaneStrainElasticData
+#include "data/DruckerPragerPlaneStrainTimeDepData.hh" // USES DruckerPragerPlaneStrainTimeDepData
+
+#include "pylith/materials/DruckerPragerPlaneStrain.hh" // USES DruckerPragerPlaneStrain
+
+#include <cstring> // USES memcpy()
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestDruckerPragerPlaneStrain );
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::materials::TestDruckerPragerPlaneStrain::setUp(void)
+{ // setUp
+ _material = new DruckerPragerPlaneStrain();
+ _matElastic = new DruckerPragerPlaneStrain();
+ _data = new DruckerPragerPlaneStrainElasticData();
+ _dataElastic = new DruckerPragerPlaneStrainElasticData();
+ setupNormalizer();
+} // setUp
+
+// ----------------------------------------------------------------------
+// Test timeStep()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::testTimeStep(void)
+{ // testTimeStep
+ DruckerPragerPlaneStrain material;
+
+ CPPUNIT_ASSERT_EQUAL(false, material._needNewJacobian);
+
+ const PylithScalar dt1 = 1.0;
+ material.timeStep(dt1);
+ CPPUNIT_ASSERT_EQUAL(dt1, material.Material::timeStep());
+ CPPUNIT_ASSERT_EQUAL(true, material.needNewJacobian());
+
+ const PylithScalar 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::TestDruckerPragerPlaneStrain::testUseElasticBehavior(void)
+{ // testUseElasticBehavior
+ DruckerPragerPlaneStrain 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::DruckerPragerPlaneStrain::_calcStressElastic ==
+ material._calcStressFn) {
+ CPPUNIT_ASSERT(
+ &pylith::materials::DruckerPragerPlaneStrain::_calcStressElastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(
+ &pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(
+ &pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastic ==
+ material._updateStateVarsFn);
+
+ material.useElasticBehavior(false);
+ CPPUNIT_ASSERT(
+ &pylith::materials::DruckerPragerPlaneStrain::_calcStressElastoplastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(
+ &pylith::materials::DruckerPragerPlaneStrain::_calcElasticConstsElastoplastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(
+ &pylith::materials::DruckerPragerPlaneStrain::_updateStateVarsElastoplastic ==
+ material._updateStateVarsFn);
+ } // if
+} // testUseElasticBehavior
+
+// ----------------------------------------------------------------------
+// Test usesHasStateVars()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::testHasStateVars(void)
+{ // testHasStateVars
+ DruckerPragerPlaneStrain material;
+ CPPUNIT_ASSERT_EQUAL(true, material.hasStateVars());
+} // testHasStateVars
+
+// ----------------------------------------------------------------------
+// Test _calcStressElastic()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::test_calcStressElastic(void)
+{ // test_calcStressElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_calcStress();
+} // test_calcStressElastic
+
+// ----------------------------------------------------------------------
+// Test calcElasticConstsElastic()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::test_calcElasticConstsElastic(void)
+{ // test_calcElasticConstsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_calcElasticConsts();
+} // test_calcElasticConstsElastic
+
+// ----------------------------------------------------------------------
+// Test _updateStateVarsElastic()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::test_updateStateVarsElastic(void)
+{ // test_updateStateVarsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(true);
+
+ test_updateStateVars();
+} // test_updateStateVarsElastic
+
+// ----------------------------------------------------------------------
+// Test _calcStressTimeDep()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::test_calcStressTimeDep(void)
+{ // test_calcStressTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new DruckerPragerPlaneStrainTimeDepData();
+
+ PylithScalar dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcStress();
+} // test_calcStressTimeDep
+
+// ----------------------------------------------------------------------
+// Test _calcElasticConstsTimeDep()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::test_calcElasticConstsTimeDep(void)
+{ // test_calcElasticConstsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new DruckerPragerPlaneStrainTimeDepData();
+
+ PylithScalar dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcElasticConsts();
+} // test_calcElasticConstsTimeDep
+
+// ----------------------------------------------------------------------
+// Test _updateStateVarsTimeDep()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::test_updateStateVarsTimeDep(void)
+{ // test_updateStateVarsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useElasticBehavior(false);
+
+ delete _dataElastic; _dataElastic = new DruckerPragerPlaneStrainTimeDepData();
+
+ PylithScalar dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_updateStateVars();
+
+} // test_updateStateVarsTimeDep
+
+// ----------------------------------------------------------------------
+// Test _stableTimeStepImplicit()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::test_stableTimeStepImplicit(void)
+{ // test_stableTimeStepImplicit
+ CPPUNIT_ASSERT(0 != _matElastic);
+
+ delete _dataElastic; _dataElastic = new DruckerPragerPlaneStrainTimeDepData();
+
+ TestElasticMaterial::test_stableTimeStepImplicit();
+} // test_stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Test hasProperty().
+void
+pylith::materials::TestDruckerPragerPlaneStrain::testHasProperty(void)
+{ // testHasProperty
+ DruckerPragerPlaneStrain material;
+
+ CPPUNIT_ASSERT(material.hasProperty("mu"));
+ CPPUNIT_ASSERT(material.hasProperty("lambda"));
+ CPPUNIT_ASSERT(material.hasProperty("density"));
+ CPPUNIT_ASSERT(material.hasProperty("alpha_yield"));
+ CPPUNIT_ASSERT(material.hasProperty("beta"));
+ CPPUNIT_ASSERT(material.hasProperty("alpha_flow"));
+ CPPUNIT_ASSERT(!material.hasProperty("aaa"));
+} // testHasProperty
+
+// ----------------------------------------------------------------------
+// Test hasStateVar().
+void
+pylith::materials::TestDruckerPragerPlaneStrain::testHasStateVar(void)
+{ // testHasStateVar
+ DruckerPragerPlaneStrain material;
+
+ CPPUNIT_ASSERT(material.hasStateVar("stress_zz_initial"));
+ CPPUNIT_ASSERT(material.hasStateVar("plastic_strain"));
+} // testHasStateVar
+
+// ----------------------------------------------------------------------
+// Test _dbToProperties().
+void
+pylith::materials::TestDruckerPragerPlaneStrain::testDBToProperties(void)
+{ // testDBToProperties
+ CPPUNIT_ASSERT(0 != _material);
+
+ TestMaterial::testDBToProperties();
+
+ CPPUNIT_ASSERT_EQUAL(false, _material->isJacobianSymmetric());
+} // testDBToProperties
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,126 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh
+ *
+ * @brief C++ TestDruckerPragerPlaneStrain object
+ *
+ * C++ unit testing for DruckerPragerPlaneStrain.
+ */
+
+#if !defined(pylith_materials_testdruckerpragerplanestrain_hh)
+#define pylith_materials_testdruckerpragerplanestrain_hh
+
+#include "TestElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class TestDruckerPragerPlaneStrain;
+ class DruckerPragerPlaneStrainElasticData;
+ class DruckerPragerPlaneStrainTimeDepData;
+ } // materials
+} // pylith
+
+/// C++ unit testing for DruckerPragerPlaneStrain
+class pylith::materials::TestDruckerPragerPlaneStrain : public TestElasticMaterial
+{ // class TestDruckerPragerPlaneStrain
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestDruckerPragerPlaneStrain );
+
+ 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 Drucker-Prager elastoplastic 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( testHasProperty );
+ CPPUNIT_TEST( testHasStateVar );
+
+ 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);
+
+ /// Test hasProperty()
+ void testHasProperty(void);
+
+ /// Test hasStateVar()
+ void testHasStateVar(void);
+
+ /// Test _dbToProperties() and check flag for symmetry of Jacobian.
+ void testDBToProperties(void);
+
+}; // class TestDruckerPragerPlaneStrain
+
+#endif // pylith_materials_testdruckerpragerplanestrain_hh
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElastic.py (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainElastic.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElastic.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElastic.py 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,228 @@
+#!/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-2011 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/DruckerPragerPlaneStrainElastic.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ DruckerPragerPlaneStrain object with elastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+import math
+
+# ----------------------------------------------------------------------
+dimension = 2
+numElasticConsts = 9
+tensorSize = 3
+
+# DruckerPragerPlaneStrainElastic class
+class DruckerPragerPlaneStrainElastic(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ DruckerPragerPlaneStrain object with elastic behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="druckerpragerplanestrainelastic"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ # import pdb
+ # pdb.set_trace()
+ numLocs = 2
+
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "friction-angle", "cohesion",
+ "dilatation-angle"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+
+ self.dbStateVarValues = ["stress-zz-initial",
+ "plastic-strain-xx",
+ "plastic-strain-yy",
+ "plastic-strain-zz",
+ "plastic-strain-xy"
+ ]
+ self.numStateVarValues = numpy.array([1, 4], dtype=numpy.int32)
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ # First case has different values for friction angle and dilatation angle.
+ frictionAngleA = math.radians(30.0)
+ dilatationAngleA = math.radians(20.0)
+ cohesionA = 3.0e5
+ strainA = [1.1e-4, 1.2e-4, 1.4e-4]
+ initialStressA = [2.1e4, 2.2e4, 2.4e4]
+ initialStrainA = [3.1e-4, 3.2e-4, 3.4e-4]
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+ stressZZInitialA = numpy.array([1.075e4], dtype=numpy.float64)
+
+ denomFrictionA = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleA))
+ denomDilatationA = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleA))
+ alphaYieldA = 2.0 * math.sin(frictionAngleA)/denomFrictionA
+ betaA = 6.0 * cohesionA * math.cos(frictionAngleA)/denomFrictionA
+ alphaFlowA = 2.0 * math.sin(dilatationAngleA)/denomDilatationA
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ # Second case has same values for friction angle and dilatation angle.
+ frictionAngleB = math.radians(25.0)
+ dilatationAngleB = math.radians(25.0)
+ cohesionB = 1.0e5
+ strainB = [4.1e-4, 4.2e-4, 4.4e-4]
+ initialStressB = [5.1e4, 5.2e4, 5.4e4]
+ initialStrainB = [6.1e-4, 6.2e-4, 6.4e-4]
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ denomFrictionB = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleB))
+ denomDilatationB = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleB))
+ alphaYieldB = 2.0 * math.sin(frictionAngleB)/denomFrictionB
+ betaB = 6.0 * cohesionB * math.cos(frictionAngleB)/denomFrictionB
+ alphaFlowB = 2.0 * math.sin(dilatationAngleB)/denomDilatationB
+ stressZZInitialB = numpy.array([2.575e4], dtype=numpy.float64)
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.densityScale = 1.0e+3
+ self.timeScale = 1.0
+
+ self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
+ frictionAngleA, cohesionA, \
+ dilatationAngleA],
+ [densityB, vsB, vpB, \
+ frictionAngleB, cohesionB, \
+ dilatationAngleB] ],
+ dtype=numpy.float64)
+ self.properties = numpy.array([ [densityA, muA, lambdaA, \
+ alphaYieldA, betaA, \
+ alphaFlowA],
+ [densityB, muB, lambdaB, \
+ alphaYieldB, betaB, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ self.dbStateVars = numpy.zeros( (numLocs, 1 + 4), dtype=numpy.float64)
+ self.dbStateVars[0, 0] = stressZZInitialA
+ self.dbStateVars[1, 0] = stressZZInitialB
+
+ self.stateVars = numpy.zeros( (numLocs, 1 + 4), dtype=numpy.float64)
+ self.stateVars[0, 0] = stressZZInitialA
+ self.stateVars[1, 0] = stressZZInitialB
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, \
+ alphaYieldA, betaA/mu0, \
+ alphaFlowA],
+ [densityB/density0, muB/mu0, lambdaB/mu0, \
+ alphaYieldB, betaB/mu0, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ stressZZInitialANondim = stressZZInitialA/mu0
+ stressZZInitialBNondim = stressZZInitialB/mu0
+
+ self.stateVarsNondim = numpy.zeros( (numLocs, 1 + 4), dtype=numpy.float64)
+ self.stateVarsNondim[0, 0] = stressZZInitialANondim
+ self.stateVarsNondim[1, 0] = stressZZInitialBNondim
+
+ 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.elasticConsts[0,:], self.stress[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA, \
+ initialStressA, initialStrainA)
+ (self.elasticConsts[1,:], self.stress[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB, \
+ initialStressB, initialStrainB)
+
+ self.dtStableImplicit = 1.0e10
+
+ stateVarsUpdatedA = numpy.array([stressZZInitialA, 0.0, 0.0, 0.0, 0.0],
+ dtype=numpy.float64)
+ stateVarsUpdatedB = numpy.array([stressZZInitialB, 0.0, 0.0, 0.0, 0.0],
+ dtype=numpy.float64)
+
+ self.stateVarsUpdated = numpy.array( [stateVarsUpdatedA, stateVarsUpdatedB],
+ dtype=numpy.float64)
+
+ return
+
+
+ def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV):
+ """
+ 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)
+
+ strain = numpy.reshape(strainV, (tensorSize,1))
+ initialStress = numpy.reshape(initialStressV, (tensorSize,1))
+ initialStrain = numpy.reshape(initialStrainV, (tensorSize,1))
+ elastic = numpy.array([ [C1111, C1122, C1112],
+ [C2211, C2222, C2212],
+ [C1211, C1222, C1212] ], dtype=numpy.float64)
+ stress = numpy.dot(elastic, strain-initialStrain) + initialStress
+ return (elasticConsts, numpy.ravel(stress))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = DruckerPragerPlaneStrainElastic()
+ app.run()
+
+
+# End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.cc (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.cc)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.cc 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,278 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerpragerplanestrainelastic.
+
+#include "DruckerPragerPlaneStrainElasticData.hh"
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_dimension = 2;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numLocs = 2;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numProperties = 6;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numStateVars = 2;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numDBProperties = 6;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numDBStateVars = 5;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numPropsQuadPt = 6;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numVarsQuadPt = 5;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_lengthScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_timeScale = 1.00000000e+00;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_pressureScale = 2.25000000e+10;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_densityScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_dtStableImplicit = 1.00000000e+10;
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::DruckerPragerPlaneStrainElasticData::_numStateVarValues[] = {
+1,
+4,
+};
+
+const char* pylith::materials::DruckerPragerPlaneStrainElasticData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"friction-angle",
+"cohesion",
+"dilatation-angle",
+};
+
+const char* pylith::materials::DruckerPragerPlaneStrainElasticData::_dbStateVarValues[] = {
+"stress-zz-initial",
+"plastic-strain-xx",
+"plastic-strain-yy",
+"plastic-strain-zz",
+"plastic-strain-xy",
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 5.23598776e-01,
+ 3.00000000e+05,
+ 3.49065850e-01,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 4.36332313e-01,
+ 1.00000000e+05,
+ 4.36332313e-01,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_dbStateVars[] = {
+ 1.07500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.57500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 2.30940108e-01,
+ 3.60000000e+05,
+ 1.48583084e-01,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 1.89338478e-01,
+ 1.21811303e+05,
+ 1.89338478e-01,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_stateVars[] = {
+ 1.07500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.57500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 2.30940108e-01,
+ 1.60000000e-05,
+ 1.48583084e-01,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 1.89338478e-01,
+ 5.41383567e-06,
+ 1.89338478e-01,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_stateVarsNondim[] = {
+ 4.77777778e-07,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 1.14444444e-06,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_strain[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_stress[] = {
+ -1.79790000e+07,
+ -1.79780000e+07,
+ -8.97600000e+06,
+ -2.25300000e+06,
+ -2.25200000e+06,
+ -1.09800000e+06,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_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 PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.40000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.40000000e+04,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_initialStrain[] = {
+ 3.10000000e-04,
+ 3.20000000e-04,
+ 3.40000000e-04,
+ 6.10000000e-04,
+ 6.20000000e-04,
+ 6.40000000e-04,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainElasticData::_stateVarsUpdated[] = {
+ 1.07500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.57500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+pylith::materials::DruckerPragerPlaneStrainElasticData::DruckerPragerPlaneStrainElasticData(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<PylithScalar*>(_dbProperties);
+ dbStateVars = const_cast<PylithScalar*>(_dbStateVars);
+ properties = const_cast<PylithScalar*>(_properties);
+ stateVars = const_cast<PylithScalar*>(_stateVars);
+ propertiesNondim = const_cast<PylithScalar*>(_propertiesNondim);
+ stateVarsNondim = const_cast<PylithScalar*>(_stateVarsNondim);
+ density = const_cast<PylithScalar*>(_density);
+ strain = const_cast<PylithScalar*>(_strain);
+ stress = const_cast<PylithScalar*>(_stress);
+ elasticConsts = const_cast<PylithScalar*>(_elasticConsts);
+ initialStress = const_cast<PylithScalar*>(_initialStress);
+ initialStrain = const_cast<PylithScalar*>(_initialStrain);
+ stateVarsUpdated = const_cast<PylithScalar*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::DruckerPragerPlaneStrainElasticData::~DruckerPragerPlaneStrainElasticData(void)
+{}
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.hh (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.hh)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainElasticData.hh 2012-02-07 21:21:03 UTC (rev 19599)
@@ -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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerpragerplanestrainelastic.
+
+#if !defined(pylith_materials_druckerpragerplanestrainelasticdata_hh)
+#define pylith_materials_druckerpragerplanestrainelasticdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class DruckerPragerPlaneStrainElasticData;
+ } // pylith
+} // materials
+
+class pylith::materials::DruckerPragerPlaneStrainElasticData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ DruckerPragerPlaneStrainElasticData(void);
+
+ /// Destructor
+ ~DruckerPragerPlaneStrainElasticData(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 PylithScalar _lengthScale;
+
+ static const PylithScalar _timeScale;
+
+ static const PylithScalar _pressureScale;
+
+ static const PylithScalar _densityScale;
+
+ static const PylithScalar _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const PylithScalar _dbProperties[];
+
+ static const PylithScalar _dbStateVars[];
+
+ static const PylithScalar _properties[];
+
+ static const PylithScalar _stateVars[];
+
+ static const PylithScalar _propertiesNondim[];
+
+ static const PylithScalar _stateVarsNondim[];
+
+ static const PylithScalar _density[];
+
+ static const PylithScalar _strain[];
+
+ static const PylithScalar _stress[];
+
+ static const PylithScalar _elasticConsts[];
+
+ static const PylithScalar _initialStress[];
+
+ static const PylithScalar _initialStrain[];
+
+ static const PylithScalar _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_druckerpragerplanestrainelasticdata_hh
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,378 @@
+#!/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-2011 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDep.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ DruckerPragerPlaneStrain object with time dependent behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+import math
+
+# ----------------------------------------------------------------------
+dimension = 2
+numElasticConsts = 9
+tensorSize = 3
+
+# DruckerPragerPlaneStrainTimeDep class
+class DruckerPragerPlaneStrainTimeDep(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ DruckerPragerPlaneStrain object with time dependent behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="druckerpragerplanestraintimedep"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ # import pdb
+ # pdb.set_trace()
+ numLocs = 2
+
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "friction-angle", "cohesion",
+ "dilatation-angle"]
+ self.propertyValues = ["density", "mu", "lambda",
+ "alpha_yield", "beta", "alpha_flow"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+
+ self.dbStateVarValues = ["stress-zz-initial",
+ "plastic-strain-xx",
+ "plastic-strain-yy",
+ "plastic-strain-zz",
+ "plastic-strain-xy"
+ ]
+ self.numStateVarValues = numpy.array([1, 4], dtype=numpy.int32)
+
+ self.dt = 2.0e5
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ # First case has same values for friction angle and dilatation angle.
+ frictionAngleA = math.radians(30.0)
+ dilatationAngleA = math.radians(20.0)
+ cohesionA = 3.0e5
+ 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
+ stressZZInitialA = 1.075e4
+
+ denomFrictionA = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleA))
+ denomDilatationA = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleA))
+ alphaYieldA = 2.0 * math.sin(frictionAngleA)/denomFrictionA
+ betaA = 6.0 * cohesionA * math.cos(frictionAngleA)/denomFrictionA
+ alphaFlowA = 2.0 * math.sin(dilatationAngleA)/denomDilatationA
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ # Second case has different values for friction angle and dilatation angle.
+ frictionAngleB = math.radians(25.0)
+ dilatationAngleB = math.radians(25.0)
+ cohesionB = 1.0e4
+ 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
+ stressZZInitialB = 2.575e4
+
+ denomFrictionB = math.sqrt(3.0) * (3.0 - math.sin(frictionAngleB))
+ denomDilatationB = math.sqrt(3.0) * (3.0 - math.sin(dilatationAngleB))
+ alphaYieldB = 2.0 * math.sin(frictionAngleB)/denomFrictionB
+ betaB = 6.0 * cohesionB * math.cos(frictionAngleB)/denomFrictionB
+ alphaFlowB = 2.0 * math.sin(dilatationAngleB)/denomDilatationB
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.densityScale = 1.0e+3
+ self.timeScale = 1.0
+
+ self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
+ frictionAngleA, cohesionA, \
+ dilatationAngleA],
+ [densityB, vsB, vpB, \
+ frictionAngleB, cohesionB, \
+ dilatationAngleB] ],
+ dtype=numpy.float64)
+ self.properties = numpy.array([ [densityA, muA, lambdaA, \
+ alphaYieldA, betaA, \
+ alphaFlowA],
+ [densityB, muB, lambdaB, \
+ alphaYieldB, betaB, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ self.dbStateVars = numpy.zeros( (numLocs, 1 + 4), dtype=numpy.float64)
+ self.dbStateVars[0, 0] = stressZZInitialA
+ self.dbStateVars[1, 0] = stressZZInitialB
+
+ stateVarsA = [stressZZInitialA, 4.1e-5, 4.2e-5, 4.4e-5, 4.5e-5]
+ stateVarsB = [stressZZInitialB, 1.1e-5, 1.2e-5, 1.4e-5, 1.5e-5]
+ plasStrainA = stateVarsA[1:]
+ plasStrainB = stateVarsB[1:]
+ self.stateVars = numpy.array([[stateVarsA], [stateVarsB] ],
+ dtype=numpy.float64)
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, \
+ alphaYieldA, betaA/mu0, \
+ alphaFlowA],
+ [densityB/density0, muB/mu0, lambdaB/mu0, \
+ alphaYieldB, betaB/mu0, \
+ alphaFlowB] ],
+ dtype=numpy.float64)
+
+ stressZZInitialANondim = stressZZInitialA/mu0
+ stressZZInitialBNondim = stressZZInitialB/mu0
+ self.stateVarsNondim = numpy.array([[stateVarsA], [stateVarsB] ],
+ dtype=numpy.float64)
+ self.stateVarsNondim[0, 0] = stressZZInitialANondim
+ self.stateVarsNondim[1, 0] = stressZZInitialBNondim
+
+ 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.stateVarsUpdated = numpy.zeros( (numLocs, 5),
+ dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts), \
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA,
+ alphaYieldA, betaA, alphaFlowA,
+ plasStrainA,
+ initialStressA, initialStrainA,
+ stateVarsA)
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB,
+ alphaYieldB, betaB, alphaFlowB,
+ plasStrainB,
+ initialStressB, initialStrainB,
+ stateVarsB)
+
+ self.dtStableImplicit = 1.0e30
+
+ return
+
+
+ def _scalarProduct(self, tensor1, tensor2):
+ """
+ Compute the scalar product of two tensors stored in vector form (length 4).
+ """
+ scalarProduct = tensor1[0] * tensor2[0] + \
+ tensor1[1] * tensor2[1] + \
+ tensor1[2] * tensor2[2] + \
+ 2.0 * tensor1[3] * tensor2[3]
+ return scalarProduct
+
+
+ def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
+ muV, lambdaV, alphaYieldV, betaV, alphaFlowV,
+ plasStrainT, 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, visStrainTpdt = self._computeStress(strainTest, muV, lambdaV,
+ alphaYieldV, betaV,
+ alphaFlowV,
+ plasStrainT,
+ initialStress,
+ initialStrain,
+ stateVars)
+ return stressTpdt[stressComp]
+
+
+ def _computeStress(self, strainTpdt, muV, lambdaV, alphaYieldV, betaV,
+ alphaFlowV, plasStrainT, initialStress, initialStrain,
+ stateVars):
+ """
+ Function to compute stresses and plastic strains.
+ """
+
+ # Constants
+ mu2 = 2.0 * muV
+ lamPlusMu = lambdaV + muV
+ bulkModulus = lambdaV + mu2/3.0
+ ae = 1.0/mu2
+ am = 1.0/(3.0 * bulkModulus)
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0], dtype=numpy.float64)
+
+ # Values from previous time step
+ meanPlasStrainT = (plasStrainT[0] + plasStrainT[1] + plasStrainT[2])/3.0
+ devPlasStrainT = plasStrainT - meanPlasStrainT * diag
+
+ # Initial stress values
+ initialStress4 = numpy.array([initialStress[0], initialStress[1],
+ stateVars[0], initialStress[2]],
+ dtype=numpy.float64)
+ meanStressInitial = (initialStress[0] + initialStress[1] + stateVars[0])/3.0
+ devStressInitial = initialStress4 - meanStressInitial * diag
+
+ # Initial strain values
+ initialStrain4 = numpy.array([initialStrain[0], initialStrain[1],
+ 0.0, initialStrain[2]], dtype=numpy.float64)
+ meanStrainInitial = (initialStrain[0] + initialStrain[1])/3.0
+ devStrainInitial = initialStrain4 - meanStrainInitial * diag
+
+ # Values for current time step
+ strainTpdt4 = numpy.array([strainTpdt[0], strainTpdt[1],
+ 0.0, strainTpdt[2]], dtype=numpy.float64)
+ meanStrainTpdt = (strainTpdt[0] + strainTpdt[1])/3.0
+ meanStrainPPTpdt = meanStrainTpdt - meanPlasStrainT - meanStrainInitial
+ strainPPTpdt = strainTpdt4 - meanStrainTpdt * diag - \
+ devPlasStrainT - devStrainInitial
+
+ # Compute trial elastic stresses and yield function to see if yield should
+ # occur.
+ trialDevStress = strainPPTpdt/ae + devStressInitial
+ trialMeanStress = meanStrainPPTpdt/am + meanStressInitial
+ yieldFunction = 3.0 * alphaYieldV * trialMeanStress + \
+ math.sqrt(0.5 * self._scalarProduct(trialDevStress,
+ trialDevStress)) - \
+ betaV
+
+ # If yield function is greater than zero, compute elastoplastic stress.
+ if (yieldFunction >= 0.0):
+ devStressInitialProd = self._scalarProduct(devStressInitial,
+ devStressInitial)
+ strainPPTpdtProd = self._scalarProduct(strainPPTpdt, strainPPTpdt)
+ d = math.sqrt(ae * ae * devStressInitialProd + \
+ 2.0 * ae * self._scalarProduct(devStressInitial, \
+ strainPPTpdt) + \
+ strainPPTpdtProd)
+ plasticMult = 2.0 * ae * am * \
+ (3.0 * alphaYieldV * meanStrainPPTpdt/am + \
+ d/(math.sqrt(2.0) * ae) - betaV)/ \
+ (6.0 * alphaYieldV * alphaFlowV * ae + am)
+ deltaMeanPlasticStrain = plasticMult * alphaFlowV
+ meanStressTpdt = (meanStrainPPTpdt - deltaMeanPlasticStrain)/am + \
+ meanStressInitial
+
+ deltaDevPlasticStrain = 0.0
+ stressTpdt4 = numpy.zeros( (4), dtype=numpy.float64)
+ stressTpdt = numpy.zeros( (tensorSize), dtype=numpy.float64)
+ plasStrainTpdt = numpy.zeros( (4), dtype=numpy.float64)
+
+ for iComp in range(4):
+ deltaDevPlasticStrain = plasticMult * (strainPPTpdt[iComp] + \
+ ae * devStressInitial[iComp])/ \
+ (math.sqrt(2.0) * d)
+ devStressTpdt = (strainPPTpdt[iComp] - deltaDevPlasticStrain)/ae + \
+ devStressInitial[iComp]
+ stressTpdt4[iComp] = devStressTpdt + diag[iComp] * meanStressTpdt
+ plasStrainTpdt[iComp] = plasStrainT[iComp] + deltaDevPlasticStrain + \
+ diag[iComp] * deltaMeanPlasticStrain
+
+ stressTpdt[0] = stressTpdt4[0]
+ stressTpdt[1] = stressTpdt4[1]
+ stressTpdt[2] = stressTpdt4[3]
+
+ return stressTpdt, plasStrainTpdt
+
+
+ def _calcStress(self, strainV, muV, lambdaV, alphaYieldV, betaV, alphaFlowV,
+ plasStrainV, initialStressV, initialStrainV, stateVars):
+ """
+ Compute stress, updated state variables and derivative of elasticity matrix.
+ """
+ import scipy.misc
+
+ # Define some numpy arrays
+ strainTpdt = numpy.array(strainV, dtype=numpy.float64)
+ plasStrainT = numpy.array(plasStrainV, dtype=numpy.float64)
+ initialStress = numpy.array(initialStressV, dtype=numpy.float64)
+ initialStrain = numpy.array(initialStrainV, dtype=numpy.float64)
+ stressTpdt, plasStrainTpdt = self._computeStress(strainTpdt, muV, lambdaV,
+ alphaYieldV, betaV,
+ alphaFlowV, plasStrainT,
+ initialStress,
+ initialStrain,
+ stateVars)
+ stateVarsUpdated = numpy.array( (stateVars[0], plasStrainTpdt[0],
+ plasStrainTpdt[1], plasStrainTpdt[2],
+ plasStrainTpdt[3]))
+ # 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 = scipy.misc.derivative(self._calcStressComponent,
+ strainTpdt[strainComp],
+ dx=derivDx,
+ args=(strainComp,
+ stressComp,
+ strainTpdt, muV, lambdaV,
+ alphaYieldV, betaV,
+ alphaFlowV,
+ plasStrainT,
+ 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 = DruckerPragerPlaneStrainTimeDep()
+ app.run()
+
+
+# End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,278 @@
+// -*- 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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerpragerplanestraintimedep.
+
+#include "DruckerPragerPlaneStrainTimeDepData.hh"
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dimension = 2;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numLocs = 2;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numProperties = 6;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numStateVars = 2;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numDBProperties = 6;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numDBStateVars = 5;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numPropsQuadPt = 6;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numVarsQuadPt = 5;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_lengthScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_timeScale = 1.00000000e+00;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_pressureScale = 2.25000000e+10;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_densityScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dtStableImplicit = 1.00000000e+30;
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numStateVarValues[] = {
+1,
+4,
+};
+
+const char* pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"friction-angle",
+"cohesion",
+"dilatation-angle",
+};
+
+const char* pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dbStateVarValues[] = {
+"stress-zz-initial",
+"plastic-strain-xx",
+"plastic-strain-yy",
+"plastic-strain-zz",
+"plastic-strain-xy",
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 5.23598776e-01,
+ 3.00000000e+05,
+ 3.49065850e-01,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 4.36332313e-01,
+ 1.00000000e+04,
+ 4.36332313e-01,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dbStateVars[] = {
+ 1.07500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.57500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 2.30940108e-01,
+ 3.60000000e+05,
+ 1.48583084e-01,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 1.89338478e-01,
+ 1.21811303e+04,
+ 1.89338478e-01,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_stateVars[] = {
+ 1.07500000e+04,
+ 4.10000000e-05,
+ 4.20000000e-05,
+ 4.40000000e-05,
+ 4.50000000e-05,
+ 2.57500000e+04,
+ 1.10000000e-05,
+ 1.20000000e-05,
+ 1.40000000e-05,
+ 1.50000000e-05,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 2.30940108e-01,
+ 1.60000000e-05,
+ 1.48583084e-01,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 1.89338478e-01,
+ 5.41383567e-07,
+ 1.89338478e-01,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_stateVarsNondim[] = {
+ 4.77777778e-07,
+ 4.77777778e-07,
+ 4.77777778e-07,
+ 4.77777778e-07,
+ 4.77777778e-07,
+ 1.14444444e-06,
+ 1.14444444e-06,
+ 1.14444444e-06,
+ 1.14444444e-06,
+ 1.14444444e-06,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_strain[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_stress[] = {
+ -4.95807185e+05,
+ -3.82197636e+05,
+ 7.08863124e+05,
+ 5.51473672e+05,
+ 5.45142209e+05,
+ -2.85351554e+05,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_elasticConsts[] = {
+ 2.44487161e+10,
+ 1.17789183e+10,
+ -1.66473214e+10,
+ 1.09997245e+10,
+ 2.10335485e+10,
+ -1.62468306e+10,
+ -1.31854142e+10,
+ -1.29851687e+10,
+ 1.38346134e+10,
+ 1.73722726e+09,
+ 2.48498883e+09,
+ -2.42123101e+09,
+ 2.48498825e+09,
+ 1.68402528e+09,
+ -2.37432763e+09,
+ -1.21061568e+09,
+ -1.18716349e+09,
+ 1.33926043e+09,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.40000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.40000000e+04,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_initialStrain[] = {
+ 3.60000000e-05,
+ 3.50000000e-05,
+ 3.30000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.60000000e-05,
+};
+
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_stateVarsUpdated[] = {
+ 1.07500000e+04,
+ 7.52490579e-05,
+ 8.37466235e-05,
+ 2.04754109e-05,
+ 9.17808195e-05,
+ 2.57500000e+04,
+ 3.24452570e-04,
+ 3.34725393e-04,
+ -7.68586819e-05,
+ 4.32915200e-04,
+};
+
+pylith::materials::DruckerPragerPlaneStrainTimeDepData::DruckerPragerPlaneStrainTimeDepData(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<PylithScalar*>(_dbProperties);
+ dbStateVars = const_cast<PylithScalar*>(_dbStateVars);
+ properties = const_cast<PylithScalar*>(_properties);
+ stateVars = const_cast<PylithScalar*>(_stateVars);
+ propertiesNondim = const_cast<PylithScalar*>(_propertiesNondim);
+ stateVarsNondim = const_cast<PylithScalar*>(_stateVarsNondim);
+ density = const_cast<PylithScalar*>(_density);
+ strain = const_cast<PylithScalar*>(_strain);
+ stress = const_cast<PylithScalar*>(_stress);
+ elasticConsts = const_cast<PylithScalar*>(_elasticConsts);
+ initialStress = const_cast<PylithScalar*>(_initialStress);
+ initialStrain = const_cast<PylithScalar*>(_initialStrain);
+ stateVarsUpdated = const_cast<PylithScalar*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::DruckerPragerPlaneStrainTimeDepData::~DruckerPragerPlaneStrainTimeDepData(void)
+{}
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.hh (from rev 19598, short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.hh)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.hh 2012-02-07 21:21:03 UTC (rev 19599)
@@ -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-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application druckerpragerplanestraintimedep.
+
+#if !defined(pylith_materials_druckerpragerplanestraintimedepdata_hh)
+#define pylith_materials_druckerpragerplanestraintimedepdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class DruckerPragerPlaneStrainTimeDepData;
+ } // pylith
+} // materials
+
+class pylith::materials::DruckerPragerPlaneStrainTimeDepData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ DruckerPragerPlaneStrainTimeDepData(void);
+
+ /// Destructor
+ ~DruckerPragerPlaneStrainTimeDepData(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 PylithScalar _lengthScale;
+
+ static const PylithScalar _timeScale;
+
+ static const PylithScalar _pressureScale;
+
+ static const PylithScalar _densityScale;
+
+ static const PylithScalar _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const PylithScalar _dbProperties[];
+
+ static const PylithScalar _dbStateVars[];
+
+ static const PylithScalar _properties[];
+
+ static const PylithScalar _stateVars[];
+
+ static const PylithScalar _propertiesNondim[];
+
+ static const PylithScalar _stateVarsNondim[];
+
+ static const PylithScalar _density[];
+
+ static const PylithScalar _strain[];
+
+ static const PylithScalar _stress[];
+
+ static const PylithScalar _elasticConsts[];
+
+ static const PylithScalar _initialStress[];
+
+ static const PylithScalar _initialStrain[];
+
+ static const PylithScalar _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_druckerpragerplanestraintimedepdata_hh
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/generate.sh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/generate.sh 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/data/generate.sh 2012-02-07 21:21:03 UTC (rev 19599)
@@ -78,6 +78,11 @@
--data.object=GenMaxwellPlaneStrainElasticData \
--data.parent=ElasticMaterialData
+ python DruckerPragerPlaneStrainElastic.py \
+ --data.namespace=pylith,materials \
+ --data.object=DruckerPragerPlaneStrainElasticData \
+ --data.parent=ElasticMaterialData
+
# 1-D ----------------------------------------------------------------
python ElasticStrain1D.py \
@@ -130,6 +135,11 @@
--data.object=DruckerPrager3DTimeDepData \
--data.parent=ElasticMaterialData
+ python DruckerPragerPlaneStrainTimeDep.py \
+ --data.namespace=pylith,materials \
+ --data.object=DruckerPragerPlaneStrainTimeDepData \
+ --data.parent=ElasticMaterialData
+
fi
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/Makefile.am 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/Makefile.am 2012-02-07 21:21:03 UTC (rev 19599)
@@ -41,7 +41,8 @@
TestGenMaxwellPlaneStrain.py \
TestGenMaxwellQpQsIsotropic3D.py \
TestPowerLaw3D.py \
- TestDruckerPrager3D.py
+ TestDruckerPrager3D.py \
+ TestDruckerPragerPlaneStrain.py
# End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py (from rev 19598, short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py 2012-02-07 21:21:03 UTC (rev 19599)
@@ -0,0 +1,102 @@
+#!/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-2011 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/materials/TestDruckerPragerPlaneStrain.py
+
+## @brief Unit testing of DruckerPragerPlaneStrain object.
+
+import unittest
+
+from pylith.materials.DruckerPragerPlaneStrain import DruckerPragerPlaneStrain
+
+# ----------------------------------------------------------------------
+class TestDruckerPragerPlaneStrain(unittest.TestCase):
+ """
+ Unit testing of DruckerPragerPlaneStrain object.
+ """
+
+ def setUp(self):
+ """
+ Setup test subject.
+ """
+ self.material = DruckerPragerPlaneStrain()
+ return
+
+
+ def test_constructor(self):
+ """
+ Test constructor.
+ """
+ self.assertEqual(2, self.material.dimension())
+ return
+
+
+ def test_fitMohrCoulomb(self):
+ """
+ Test useElasticBehavior().
+ """
+ self.material.fitMohrCoulomb(self.material.MOHR_COULOMB_MIDDLE)
+ return
+
+
+ def test_useElasticBehavior(self):
+ """
+ Test useElasticBehavior().
+ """
+ self.material.useElasticBehavior(False)
+ return
+
+
+ def testHasStateVars(self):
+ self.failUnless(self.material.hasStateVars())
+ return
+
+
+ def testTensorSize(self):
+ self.assertEqual(3, self.material.tensorSize())
+ return
+
+
+ def testNeedNewJacobian(self):
+ """
+ Test needNewJacobian().
+ """
+ # Default should be False.
+ self.failIf(self.material.needNewJacobian())
+
+ # Should require a new Jacobian even if time step is the same.
+ self.material.timeStep(1.0)
+ self.failUnless(self.material.needNewJacobian())
+ self.material.timeStep(2.0)
+ self.failUnless(self.material.needNewJacobian())
+
+ self.material.timeStep(2.0)
+ self.failUnless(self.material.needNewJacobian())
+ return
+
+
+ def test_factory(self):
+ """
+ Test factory method.
+ """
+ from pylith.materials.DruckerPragerPlaneStrain import material
+ m = material()
+ return
+
+
+# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/testmaterials.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/testmaterials.py 2012-02-07 21:16:56 UTC (rev 19598)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/testmaterials.py 2012-02-07 21:21:03 UTC (rev 19599)
@@ -98,6 +98,9 @@
from TestDruckerPrager3D import TestDruckerPrager3D
suite.addTest(unittest.makeSuite(TestDruckerPrager3D))
+ from TestDruckerPragerPlaneStrain import TestDruckerPragerPlaneStrain
+ suite.addTest(unittest.makeSuite(TestDruckerPragerPlaneStrain))
+
from TestMaterial import TestMaterial
suite.addTest(unittest.makeSuite(TestMaterial))
More information about the CIG-COMMITS
mailing list