[cig-commits] r19724 - short/3D/PyLith/trunk/libsrc/pylith/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Sun Mar 4 18:50:36 PST 2012
Author: willic3
Date: 2012-03-04 18:50:35 -0800 (Sun, 04 Mar 2012)
New Revision: 19724
Modified:
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
Log:
Changed elastic constants computations when tensile failure occurs.
We still need to put in a switch as to whether we allow this or not.
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-05 02:48:50 UTC (rev 19723)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-05 02:50:35 UTC (rev 19724)
@@ -904,9 +904,13 @@
(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 =
- std::min(sqrt(2.0)*d,
- plasticFac * (meanStrainFac * meanStrainPPTpdt + dFac * d - beta));
+ const PylithScalar testMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ const PylithScalar plasticMult = std::min(sqrt(2.0)*d, testMult);
+
+ bool tensileYield = false;
+ if (plasticMult == sqrt(2.0) * d)
+ tensileYield = true;
// Define some constants, vectors, and matrices.
const PylithScalar third = 1.0/3.0;
@@ -925,6 +929,7 @@
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.
@@ -935,15 +940,25 @@
2.0 * vec1[3]/d,
2.0 * vec1[4]/d,
2.0 * vec1[5]/d};
- const PylithScalar dLambdadEpsilon[tensorSize] = {
- plasticFac * (alphaYield/am + dFac * dDdEpsilon[0]),
- plasticFac * (alphaYield/am + dFac * dDdEpsilon[1]),
- plasticFac * (alphaYield/am + dFac * dDdEpsilon[2]),
- plasticFac * ( dFac * dDdEpsilon[3]),
- plasticFac * ( dFac * dDdEpsilon[4]),
- plasticFac * ( dFac * dDdEpsilon[5])
- };
- const PylithScalar dFac2 = 1.0/(sqrt(2.0) * d);
+ PylithScalar dLambdadEpsilon[tensorSize];
+ if (tensileYield) {
+ dLambdadEpsilon[0] = sqrt(2.0) * dDdEpsilon[0];
+ dLambdadEpsilon[1] = sqrt(2.0) * dDdEpsilon[1];
+ dLambdadEpsilon[2] = sqrt(2.0) * dDdEpsilon[2];
+ dLambdadEpsilon[3] = sqrt(2.0) * dDdEpsilon[3];
+ dLambdadEpsilon[4] = sqrt(2.0) * dDdEpsilon[4];
+ dLambdadEpsilon[5] = sqrt(2.0) * dDdEpsilon[5];
+ } else {
+ dLambdadEpsilon[0] = plasticFac *
+ (alphaYield/am + dFac * dDdEpsilon[0]);
+ dLambdadEpsilon[1] = plasticFac *
+ (alphaYield/am + dFac * dDdEpsilon[1]);
+ dLambdadEpsilon[2] = plasticFac *
+ (alphaYield/am + dFac * dDdEpsilon[2]);
+ dLambdadEpsilon[3] = plasticFac * dFac * dDdEpsilon[3];
+ dLambdadEpsilon[4] = plasticFac * dFac * dDdEpsilon[4];
+ dLambdadEpsilon[5] = plasticFac * dFac * dDdEpsilon[5];
+ } // else
for (int iComp=0; iComp < tensorSize; ++iComp) {
for (int jComp=0; jComp < tensorSize; ++jComp) {
int iCount = jComp + tensorSize * iComp;
@@ -958,14 +973,12 @@
} // for
} // for
} else {
- const PylithScalar dLambdadEpsilon[tensorSize] = {
- plasticFac * (alphaYield/am),
- plasticFac * (alphaYield/am),
- plasticFac * (alphaYield/am),
- 0.0,
- 0.0,
- 0.0,
- };
+ const PylithScalar dLambdadEpsilon[tensorSize] = {0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.0};
for (int iComp=0; iComp < tensorSize; ++iComp) {
for (int jComp=0; jComp < tensorSize; ++jComp) {
int iCount = jComp + tensorSize * iComp;
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-03-05 02:48:50 UTC (rev 19723)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-03-05 02:50:35 UTC (rev 19724)
@@ -926,10 +926,14 @@
(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 =
- std::min(sqrt(2.0)*d,
- plasticFac * (meanStrainFac * meanStrainPPTpdt + dFac * d - beta));
+ const PylithScalar testMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ const PylithScalar plasticMult = std::min(sqrt(2.0)*d, testMult);
+ bool tensileYield = false;
+ if (plasticMult == sqrt(2.0) * d)
+ tensileYield = true;
+
// Define some constants, vectors, and matrices.
const PylithScalar third = 1.0/3.0;
const PylithScalar dEdEpsilon[3][3] = {
@@ -950,29 +954,42 @@
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])};
+
+ PylithScalar dLambdadEpsilon[tensorSize];
+ if (tensileYield) {
+ dLambdadEpsilon[0] = sqrt(2.0) * dDdEpsilon[0];
+ dLambdadEpsilon[1] = sqrt(2.0) * dDdEpsilon[1];
+ dLambdadEpsilon[2] = sqrt(2.0) * dDdEpsilon[2];
+ } else {
+ dLambdadEpsilon[0] = plasticFac *
+ (alphaYield/am + dFac * dDdEpsilon[0]);
+ dLambdadEpsilon[1] = plasticFac *
+ (alphaYield/am + dFac * dDdEpsilon[1]);
+ dLambdadEpsilon[2] = plasticFac * dFac * dDdEpsilon[2];
+ } // else
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) +
+ 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;
+ elasticConsts[iCount] = (dEdEpsilon[iComp][jComp] -
+ dDeltaEdEpsilon)/ae +
+ diag[iComp] * (third * diag[jComp] -
+ alphaFlow * dLambdadEpsilon[jComp])/am;
} // for
} // for
} else {
- const PylithScalar dLambdadEpsilon[3] = {
- plasticFac * (alphaYield/am),
- plasticFac * (alphaYield/am),
- 0.0};
+ const PylithScalar dLambdadEpsilon[3] = {0.0,
+ 0.0,
+ 0.0};
for (int iComp=0; iComp < tensorSize; ++iComp) {
for (int jComp=0; jComp < tensorSize; ++jComp) {
int iCount = jComp + tensorSize * iComp;
elasticConsts[iCount] = (dEdEpsilon[iComp][jComp])/ae +
- diag[iComp] * (third * diag[jComp] - alphaFlow * dLambdadEpsilon[jComp])/am;
+ diag[iComp] * (third * diag[jComp] -
+ alphaFlow * dLambdadEpsilon[jComp])/am;
} // for
} // for
} // if/else
More information about the CIG-COMMITS
mailing list