[cig-commits] r19772 - in short/3D/PyLith/trunk: libsrc/pylith/materials modulesrc/materials pylith/materials unittests/libtests/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Mar 13 17:18:25 PDT 2012
Author: willic3
Date: 2012-03-13 17:18:25 -0700 (Tue, 13 Mar 2012)
New Revision: 19772
Modified:
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.hh
short/3D/PyLith/trunk/modulesrc/materials/DruckerPrager3D.i
short/3D/PyLith/trunk/pylith/materials/DruckerPrager3D.py
short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh
Log:
Added switch to control whether tensile yield is allowed or not.
Need to do this for plane strain and need to update Python tests.
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-13 16:41:07 UTC (rev 19771)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-14 00:18:25 UTC (rev 19772)
@@ -154,6 +154,7 @@
_DruckerPrager3D::dbStateVars,
_DruckerPrager3D::numDBStateVars)),
_fitMohrCoulomb(MOHR_COULOMB_INSCRIBED),
+ _allowTensileYield(false),
_calcElasticConstsFn(0),
_calcStressFn(0),
_updateStateVarsFn(0)
@@ -168,6 +169,14 @@
} // destructor
// ----------------------------------------------------------------------
+// Set boolean for whether tensile yield is allowed.
+void
+pylith::materials::DruckerPrager3D::allowTensileYield(const bool flag)
+{ // allowTensileYield
+ _allowTensileYield = flag;
+} // allowTensileYield
+
+// ----------------------------------------------------------------------
// Set fit to Mohr-Coulomb surface.
void
pylith::materials::DruckerPrager3D::fitMohrCoulomb(FitMohrCoulombEnum value)
@@ -246,7 +255,8 @@
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));
+ 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;
@@ -260,7 +270,8 @@
} // 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));
+ 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;
@@ -560,12 +571,16 @@
const PylithScalar e22 = totalStrain[1];
const PylithScalar e33 = totalStrain[2];
const PylithScalar meanStrainTpdt = (e11 + e22 + e33)/3.0;
- const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
+ const PylithScalar meanStrainPPTpdt =
+ meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSize] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5],
@@ -581,9 +596,12 @@
strainPPTpdt[4]/ae + devStressInitial[4],
strainPPTpdt[5]/ae + devStressInitial[5],
};
- const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
- const PylithScalar stressInvar2 = sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+ const PylithScalar trialMeanStress =
+ meanStrainPPTpdt/am + meanStressInitial;
+ const PylithScalar stressInvar2 =
+ sqrt(0.5 * scalarProduct3D(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;
@@ -603,18 +621,33 @@
const PylithScalar d =
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
- const PylithScalar plasticMult =
- std::min(sqrt(2.0)*d,
- 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
- (6.0 * alphaYield * alphaFlow * ae + am));
+
+ PylithScalar plasticMult = 0.0;
+ if(_allowTensileYield) {
+ plasticMult =
+ std::min(sqrt(2.0)*d,
+ 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae)
+ - beta)/ (6.0 * alphaYield * alphaFlow * ae + am));
+ } else {
+ plasticMult =
+ 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ } // if/else
+
const PylithScalar meanStressTpdt =
(meanStrainPPTpdt - plasticMult * alphaFlow)/am + meanStressInitial;
PylithScalar deltaDevPlasticStrain = 0.0;
PylithScalar devStressTpdt = 0.0;
- if (d > 0.0) {
+ if (d > 0.0 || !_allowTensileYield) {
for (int iComp=0; iComp < tensorSize; ++iComp) {
- deltaDevPlasticStrain = plasticMult * (strainPPTpdt[iComp] + ae * devStressInitial[iComp]) / (sqrt(2.0) * d);
- devStressTpdt = (strainPPTpdt[iComp] - deltaDevPlasticStrain) / ae + devStressInitial[iComp];
+ deltaDevPlasticStrain =
+ plasticMult * (strainPPTpdt[iComp] + ae * devStressInitial[iComp])/
+ (sqrt(2.0) * d);
+ devStressTpdt =
+ (strainPPTpdt[iComp] - deltaDevPlasticStrain)/ae +
+ devStressInitial[iComp];
stress[iComp] = devStressTpdt + diag[iComp] * meanStressTpdt;
} // for
} else {
@@ -650,12 +683,18 @@
stateVars[s_plasticStrain+5],
};
- const PylithScalar e11 = totalStrain[0] - plasticStrainTpdt[0] - initialStrain[0];
- const PylithScalar e22 = totalStrain[1] - plasticStrainTpdt[1] - initialStrain[1];
- const PylithScalar e33 = totalStrain[2] - plasticStrainTpdt[2] - initialStrain[2];
- const PylithScalar e12 = totalStrain[3] - plasticStrainTpdt[3] - initialStrain[3];
- const PylithScalar e23 = totalStrain[4] - plasticStrainTpdt[4] - initialStrain[4];
- const PylithScalar e13 = totalStrain[5] - plasticStrainTpdt[5] - initialStrain[5];
+ const PylithScalar e11 =
+ totalStrain[0] - plasticStrainTpdt[0] - initialStrain[0];
+ const PylithScalar e22 =
+ totalStrain[1] - plasticStrainTpdt[1] - initialStrain[1];
+ const PylithScalar e33 =
+ totalStrain[2] - plasticStrainTpdt[2] - initialStrain[2];
+ const PylithScalar e12 =
+ totalStrain[3] - plasticStrainTpdt[3] - initialStrain[3];
+ const PylithScalar e23 =
+ totalStrain[4] - plasticStrainTpdt[4] - initialStrain[4];
+ const PylithScalar e13 =
+ totalStrain[5] - plasticStrainTpdt[5] - initialStrain[5];
const PylithScalar traceStrainTpdt = e11 + e22 + e33;
const PylithScalar s123 = lambda * traceStrainTpdt;
@@ -857,9 +896,12 @@
meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSize] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5],
@@ -904,11 +946,18 @@
(6.0 * alphaYield * alphaFlow * ae + am);
const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
- const PylithScalar testMult = plasticFac *
- (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
- const bool tensileYield = (sqrt(2.0)*d < testMult) ? true : false;
- const PylithScalar plasticMult = (tensileYield) ? sqrt(2.0)*d : testMult;
+ PylithScalar plasticMult = 0.0;
+ bool tensileYield = false;
+ if (_allowTensileYield) {
+ const PylithScalar testMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ tensileYield = (sqrt(2.0)*d < testMult) ? true : false;
+ plasticMult = (tensileYield) ? sqrt(2.0)*d : testMult;
+ } else {
+ plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ } // if/else
// Define some constants, vectors, and matrices.
const PylithScalar third = 1.0/3.0;
@@ -1148,9 +1197,12 @@
meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSize] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5]
@@ -1167,8 +1219,29 @@
strainPPTpdt[5]/ae + devStressInitial[5]
};
const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
- const PylithScalar stressInvar2 = sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+ const PylithScalar stressInvar2 =
+ sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
+ const PylithScalar yieldFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+
+ if (!_allowTensileYield) {
+ const PylithScalar feasibleFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+ // If mean stress is too negative to project back to the yield surface,
+ // throw an exception.
+ if (feasibleFunction < 0.0) {
+ std::ostringstream msg;
+ msg << "Infeasible stress state - cannot project back to yield surface.\n"
+ << " alphaYield: " << alphaYield << "\n"
+ << " alphaFlow: " << alphaFlow << "\n"
+ << " beta: " << beta << "\n"
+ << " trialMeanStress: " << trialMeanStress << "\n"
+ << " stressInvar2: " << stressInvar2 << "\n"
+ << " yieldFunction: " << yieldFunction << "\n"
+ << " feasibleFunction: " << feasibleFunction << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
#if 0 // DEBUGGING
std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
std::cout << " alphaYield: " << alphaYield << std::endl;
@@ -1189,13 +1262,25 @@
const PylithScalar d =
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
- const PylithScalar plasticMult =
- std::min(sqrt(2.0)*d,
- 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
- (6.0 * alphaYield * alphaFlow * ae + am));
+ 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);
+
+ PylithScalar plasticMult = 0.0;
+ if (_allowTensileYield) {
+ plasticMult = std::min(sqrt(2.0) * d,
+ plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d -
+ beta));
+ } else {
+ plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ } // if/else
+
const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
PylithScalar deltaDevPlasticStrain = 0.0;
- if (d > 0.0) {
+ if (d > 0.0 || !_allowTensileYield) {
for (int iComp=0; iComp < tensorSize; ++iComp) {
deltaDevPlasticStrain = plasticMult *(strainPPTpdt[iComp] +
ae * devStressInitial[iComp])/
@@ -1205,7 +1290,8 @@
} // for
} else {
for (int iComp=0; iComp < tensorSize; ++iComp) {
- stateVars[s_plasticStrain+iComp] += diag[iComp] * deltaMeanPlasticStrain;
+ stateVars[s_plasticStrain+iComp] +=
+ diag[iComp] * deltaMeanPlasticStrain;
} // for
} // if/else
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.hh 2012-03-13 16:41:07 UTC (rev 19771)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.hh 2012-03-14 00:18:25 UTC (rev 19772)
@@ -27,7 +27,7 @@
// Include directives ---------------------------------------------------
#include "ElasticMaterial.hh" // ISA ElasticMaterial
-// Powerlaw3D -----------------------------------------------------------
+// DruckerPrager3D ------------------------------------------------------
/** @brief 3-D, isotropic, Drucker-Prager elastic/perfectly plastic material.
*
* The physical properties are specified using density, shear-wave
@@ -67,6 +67,12 @@
*/
void fitMohrCoulomb(FitMohrCoulombEnum value);
+ /** Set flag for whether to allow tensile yield.
+ *
+ * @param flag True if tensile yield is allowed.
+ */
+ void allowTensileYield(const bool flag);
+
/** Set current time step.
*
* @param dt Current time step.
@@ -499,6 +505,9 @@
/// Fit to Mohr Coulomb surface
FitMohrCoulombEnum _fitMohrCoulomb;
+ /// Whether to allow tensile yield
+ bool _allowTensileYield;
+
static const int p_density;
static const int p_mu;
static const int p_lambda;
Modified: short/3D/PyLith/trunk/modulesrc/materials/DruckerPrager3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/DruckerPrager3D.i 2012-03-13 16:41:07 UTC (rev 19771)
+++ short/3D/PyLith/trunk/modulesrc/materials/DruckerPrager3D.i 2012-03-14 00:18:25 UTC (rev 19772)
@@ -50,6 +50,12 @@
* @param value Mohr-Coulomb surface match type.
*/
void fitMohrCoulomb(FitMohrCoulombEnum value);
+
+ /** Set boolean for whether tensile yield is allowed.
+ *
+ * @param flag True if tensile yield is allowed, false otherwise.
+ */
+ void allowTensileYield(const bool flag);
/** Set current time step.
*
Modified: short/3D/PyLith/trunk/pylith/materials/DruckerPrager3D.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/DruckerPrager3D.py 2012-03-13 16:41:07 UTC (rev 19771)
+++ short/3D/PyLith/trunk/pylith/materials/DruckerPrager3D.py 2012-03-14 00:18:25 UTC (rev 19772)
@@ -57,6 +57,8 @@
##
## \b Properties
## @li \b fit_mohr_coulomb Fit to Mohr-Coulomb yield surface.
+ ## @li \b allow_tensile_yield If true, allow yield beyond tensile strength;
+ ## otherwise an exception occurs for excessive tensile sttress.
##
## \b Facilities
## @li None
@@ -67,6 +69,10 @@
fitMohrCoulomb = pyre.inventory.str("fit_mohr_coulomb", default="inscribed",
validator=validateFitMohrCoulomb)
fitMohrCoulomb.meta['tip'] = "Fit to Mohr-Coulomb yield surface."
+ allowTensileYield = pyre.inventory.bool("allow_tensile_yield",
+ default=False)
+ allowTensileYield.meta['tip'] = "If true, allow yield beyond tensile " \
+ "strength; otherwise, an exception is thrown."
# PUBLIC METHODS /////////////////////////////////////////////////////
@@ -104,6 +110,8 @@
else:
raise ValueError("Unknown fit to Mohr-Coulomb yield surface.")
ModuleDruckerPrager3D.fitMohrCoulomb(self, fitEnum)
+ ModuleDruckerPrager3D.allowTensileYield(self,
+ self.inventory.allowTensileYield)
return
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc 2012-03-13 16:41:07 UTC (rev 19771)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc 2012-03-14 00:18:25 UTC (rev 19772)
@@ -36,7 +36,12 @@
pylith::materials::TestDruckerPrager3D::setUp(void)
{ // setUp
_material = new DruckerPrager3D();
- _matElastic = new DruckerPrager3D();
+ // _matElastic = new DruckerPrager3D();
+ DruckerPrager3D* matTmp = new DruckerPrager3D();
+ CPPUNIT_ASSERT(matTmp);
+ matTmp->allowTensileYield(true);
+ _matElastic = matTmp;
+
_data = new DruckerPrager3DElasticData();
_dataElastic = new DruckerPrager3DElasticData();
setupNormalizer();
@@ -92,6 +97,18 @@
} // testUseElasticBehavior
// ----------------------------------------------------------------------
+// Test allowTensileYield()
+void
+pylith::materials::TestDruckerPrager3D::testAllowTensileYield(void)
+{ // testAllowTensileYield
+ DruckerPrager3D material;
+ CPPUNIT_ASSERT(!material._allowTensileYield);
+
+ material.allowTensileYield(true);
+ CPPUNIT_ASSERT(material._allowTensileYield);
+} // testAllowTensileYield
+
+// ----------------------------------------------------------------------
// Test usesHasStateVars()
void
pylith::materials::TestDruckerPrager3D::testHasStateVars(void)
@@ -140,6 +157,9 @@
{ // test_calcStressTimeDep
CPPUNIT_ASSERT(0 != _matElastic);
_matElastic->useElasticBehavior(false);
+ // pylith::materials::DruckerPrager3D->allowTensileYield(true);
+ // _matElastic->allowTensileYield(true);
+ // _matElastic._allowTensileYield = true;
delete _dataElastic; _dataElastic = new DruckerPrager3DTimeDepData();
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh 2012-03-13 16:41:07 UTC (rev 19771)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.hh 2012-03-14 00:18:25 UTC (rev 19772)
@@ -59,6 +59,7 @@
// Need to test Drucker-Prager elastoplastic specific behavior.
CPPUNIT_TEST( testTimeStep );
CPPUNIT_TEST( testUseElasticBehavior );
+ CPPUNIT_TEST( testAllowTensileYield );
CPPUNIT_TEST( testHasStateVars );
CPPUNIT_TEST( test_calcStressElastic );
@@ -85,6 +86,9 @@
/// Test useElasticBehavior()
void testUseElasticBehavior(void);
+ /// Test allowTensileYield()
+ void testAllowTensileYield(void);
+
/// Test hasStateVars()
void testHasStateVars(void);
More information about the CIG-COMMITS
mailing list