[cig-commits] r19585 - in short/3D/PyLith/trunk: libsrc/pylith libsrc/pylith/materials modulesrc/materials pylith pylith/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu Feb 2 16:53:04 PST 2012
Author: willic3
Date: 2012-02-02 16:53:03 -0800 (Thu, 02 Feb 2012)
New Revision: 19585
Added:
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc
short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i
short/3D/PyLith/trunk/pylith/materials/DruckerPragerPlaneStrain.py
Modified:
short/3D/PyLith/trunk/libsrc/pylith/Makefile.am
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.hh
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.icc
short/3D/PyLith/trunk/libsrc/pylith/materials/Makefile.am
short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh
short/3D/PyLith/trunk/modulesrc/materials/Makefile.am
short/3D/PyLith/trunk/pylith/Makefile.am
Log:
Initial version of Drucker-Prager plane strain.
Everything compiles, but I need to do unit tests and full-scale tests.
Modified: short/3D/PyLith/trunk/libsrc/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/Makefile.am 2012-02-02 22:14:58 UTC (rev 19584)
+++ short/3D/PyLith/trunk/libsrc/pylith/Makefile.am 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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 \
Added: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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 double 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[3] - initialStrain[3];
+
+ 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
Added: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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
Added: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.icc 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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/trunk/libsrc/pylith/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.hh 2012-02-02 22:14:58 UTC (rev 19584)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.hh 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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/trunk/libsrc/pylith/materials/ElasticMaterial.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.icc 2012-02-02 22:14:58 UTC (rev 19584)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.icc 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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/trunk/libsrc/pylith/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Makefile.am 2012-02-02 22:14:58 UTC (rev 19584)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Makefile.am 2012-02-03 00:53:03 UTC (rev 19585)
@@ -42,6 +42,8 @@
PowerLaw3D.icc \
DruckerPrager3D.hh \
DruckerPrager3D.icc \
+ DruckerPragerPlaneStrain.hh \
+ DruckerPragerPlaneStrain.icc \
Metadata.hh \
Metadata.icc \
Material.hh \
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh 2012-02-02 22:14:58 UTC (rev 19584)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/materialsfwd.hh 2012-02-03 00:53:03 UTC (rev 19585)
@@ -48,6 +48,7 @@
class GenMaxwellQpQsIsotropic3D;
class PowerLaw3D;
class DruckerPrager3D;
+ class DruckerPragerPlaneStrain;
class EffectiveStress;
class ViscoelasticMaxwell;
Added: short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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/trunk/modulesrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2012-02-02 22:14:58 UTC (rev 19584)
+++ short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2012-02-02 22:14:58 UTC (rev 19584)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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 \
Added: short/3D/PyLith/trunk/pylith/materials/DruckerPragerPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/DruckerPragerPlaneStrain.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/materials/DruckerPragerPlaneStrain.py 2012-02-03 00:53:03 UTC (rev 19585)
@@ -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
More information about the CIG-COMMITS
mailing list