[cig-commits] r19775 - in short/3D/PyLith/trunk: libsrc/pylith/materials modulesrc/materials unittests/libtests/materials unittests/pytests/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Mar 13 18:47:03 PDT 2012
Author: willic3
Date: 2012-03-13 18:47:03 -0700 (Tue, 13 Mar 2012)
New Revision: 19775
Modified:
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh
short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i
short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh
short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPrager3D.py
short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py
Log:
Put in switch for whether to allow tensile yield for DPPlaneStrain.
A little cleanup of DP3D, and also put in Python tests.
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-14 01:47:03 UTC (rev 19775)
@@ -949,14 +949,17 @@
PylithScalar plasticMult = 0.0;
bool tensileYield = false;
+ PylithScalar dFac2 = 0.0;
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;
+ dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
} else {
plasticMult = plasticFac *
(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ dFac2 = 1.0/(sqrt(2.0) * d);
} // if/else
// Define some constants, vectors, and matrices.
@@ -976,7 +979,6 @@
strainPPTpdt[4] + ae * devStressInitial[4],
strainPPTpdt[5] + ae * devStressInitial[5]
};
- const PylithScalar dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
PylithScalar dDeltaEdEpsilon = 0.0;
// Compute elasticity matrix.
@@ -1242,6 +1244,7 @@
throw std::runtime_error(msg.str());
} // if
} // if
+
#if 0 // DEBUGGING
std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
std::cout << " alphaYield: " << alphaYield << std::endl;
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-03-14 01:47:03 UTC (rev 19775)
@@ -160,6 +160,7 @@
_DruckerPragerPlaneStrain::dbStateVars,
_DruckerPragerPlaneStrain::numDBStateVars)),
_fitMohrCoulomb(MOHR_COULOMB_INSCRIBED),
+ _allowTensileYield(false),
_calcElasticConstsFn(0),
_calcStressFn(0),
_updateStateVarsFn(0)
@@ -174,6 +175,14 @@
} // destructor
// ----------------------------------------------------------------------
+// Set boolean for whether tensile yield is allowed.
+void
+pylith::materials::DruckerPragerPlaneStrain::allowTensileYield(const bool flag)
+{ // allowTensileYield
+ _allowTensileYield = flag;
+} // allowTensileYield
+
+// ----------------------------------------------------------------------
// Set fit to Mohr-Coulomb surface.
void
pylith::materials::DruckerPragerPlaneStrain::fitMohrCoulomb(
@@ -304,8 +313,9 @@
// ----------------------------------------------------------------------
// Nondimensionalize properties.
void
-pylith::materials::DruckerPragerPlaneStrain::_nondimProperties(PylithScalar* const values,
- const int nvalues) const
+pylith::materials::DruckerPragerPlaneStrain::_nondimProperties(
+ PylithScalar* const values,
+ const int nvalues) const
{ // _nondimProperties
assert(_normalizer);
assert(values);
@@ -329,8 +339,9 @@
// ----------------------------------------------------------------------
// Dimensionalize properties.
void
-pylith::materials::DruckerPragerPlaneStrain::_dimProperties(PylithScalar* const values,
- const int nvalues) const
+pylith::materials::DruckerPragerPlaneStrain::_dimProperties(
+ PylithScalar* const values,
+ const int nvalues) const
{ // _dimProperties
assert(_normalizer);
assert(values);
@@ -578,13 +589,16 @@
// Values for current time step
const PylithScalar meanStrainTpdt = (totalStrain[0] + totalStrain[1])/3.0;
- const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
+ const PylithScalar meanStrainPPTpdt =
+ meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
// devStrainPPTpdt
const PylithScalar strainPPTpdt[tensorSizePS] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ 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]
};
@@ -596,11 +610,12 @@
strainPPTpdt[2]/ae + devStressInitial[2],
strainPPTpdt[3]/ae + devStressInitial[3]
};
- const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
+ const PylithScalar trialMeanStress =
+ meanStrainPPTpdt/am + meanStressInitial;
const PylithScalar stressInvar2 =
sqrt(0.5 * scalarProduct2DPS(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
- stressInvar2 - beta;
+ 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;
@@ -621,11 +636,23 @@
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct2DPS(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 meanStressTpdt = (meanStrainPPTpdt - plasticMult * alphaFlow) / am + meanStressInitial;
+
+ 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;
PylithScalar totalStress[tensorSizePS];
@@ -653,10 +680,14 @@
<< std::endl;
#endif
- if (d > 0.0) {
+ if (d > 0.0 || !_allowTensileYield) {
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];
+ 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
} else {
@@ -883,9 +914,12 @@
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[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[2] - devPlasticStrainT[3] - devStrainInitial[3],
};
@@ -926,11 +960,22 @@
(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;
+ PylithScalar dFac2 = 0.0;
+ 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;
+ dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
+ } else {
+ plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ dFac2 = 1.0/(sqrt(2.0) * d);
+ } // if/else
+
// Define some constants, vectors, and matrices.
const PylithScalar third = 1.0/3.0;
const PylithScalar dEdEpsilon[3][3] = {
@@ -943,7 +988,6 @@
strainPPTpdt[3] + ae * devStressInitial[3],
};
- const PylithScalar dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
PylithScalar dDeltaEdEpsilon = 0.0;
// Compute elasticity matrix.
@@ -1129,9 +1173,11 @@
meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSizePS] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- 0.0 - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ 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]
};
@@ -1145,8 +1191,28 @@
const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
const PylithScalar stressInvar2 =
sqrt(0.5 * scalarProduct2DPS(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
- stressInvar2 - beta;
+ 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;
@@ -1168,13 +1234,25 @@
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct2DPS(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 < tensorSizePS; ++iComp) {
deltaDevPlasticStrain = plasticMult *(strainPPTpdt[iComp] +
ae * devStressInitial[iComp])/
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh 2012-03-14 01:47:03 UTC (rev 19775)
@@ -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/DruckerPragerPlaneStrain.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/modulesrc/materials/DruckerPragerPlaneStrain.i 2012-03-14 01:47:03 UTC (rev 19775)
@@ -51,6 +51,12 @@
*/
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.
*
* @param dt Current time step.
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPrager3D.cc 2012-03-14 01:47:03 UTC (rev 19775)
@@ -36,7 +36,6 @@
pylith::materials::TestDruckerPrager3D::setUp(void)
{ // setUp
_material = new DruckerPrager3D();
- // _matElastic = new DruckerPrager3D();
DruckerPrager3D* matTmp = new DruckerPrager3D();
CPPUNIT_ASSERT(matTmp);
matTmp->allowTensileYield(true);
@@ -157,9 +156,6 @@
{ // 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/TestDruckerPragerPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc 2012-03-14 01:47:03 UTC (rev 19775)
@@ -36,7 +36,11 @@
pylith::materials::TestDruckerPragerPlaneStrain::setUp(void)
{ // setUp
_material = new DruckerPragerPlaneStrain();
- _matElastic = new DruckerPragerPlaneStrain();
+ DruckerPragerPlaneStrain* matTmp = new DruckerPragerPlaneStrain();
+ CPPUNIT_ASSERT(matTmp);
+ matTmp->allowTensileYield(true);
+ _matElastic = matTmp;
+
_data = new DruckerPragerPlaneStrainElasticData();
_dataElastic = new DruckerPragerPlaneStrainElasticData();
setupNormalizer();
@@ -98,6 +102,18 @@
} // testUseElasticBehavior
// ----------------------------------------------------------------------
+// Test allowTensileYield()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::testAllowTensileYield(void)
+{ // testAllowTensileYield
+ DruckerPragerPlaneStrain material;
+ CPPUNIT_ASSERT(!material._allowTensileYield);
+
+ material.allowTensileYield(true);
+ CPPUNIT_ASSERT(material._allowTensileYield);
+} // testAllowTensileYield
+
+// ----------------------------------------------------------------------
// Test usesHasStateVars()
void
pylith::materials::TestDruckerPragerPlaneStrain::testHasStateVars(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh 2012-03-14 01:47:03 UTC (rev 19775)
@@ -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);
Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPrager3D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPrager3D.py 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPrager3D.py 2012-03-14 01:47:03 UTC (rev 19775)
@@ -46,6 +46,14 @@
return
+ def test_allowTensileYield(self):
+ """
+ Test allowTensileYield().
+ """
+ self.material.allowTensileYield(True)
+ return
+
+
def test_fitMohrCoulomb(self):
"""
Test useElasticBehavior().
Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py 2012-03-14 01:07:01 UTC (rev 19774)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py 2012-03-14 01:47:03 UTC (rev 19775)
@@ -46,6 +46,14 @@
return
+ def test_allowTensileYield(self):
+ """
+ Test allowTensileYield().
+ """
+ self.material.allowTensileYield(True)
+ return
+
+
def test_fitMohrCoulomb(self):
"""
Test useElasticBehavior().
More information about the CIG-COMMITS
mailing list