[cig-commits] r16265 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon Feb 15 20:56:43 PST 2010
Author: willic3
Date: 2010-02-15 20:56:43 -0800 (Mon, 15 Feb 2010)
New Revision: 16265
Added:
short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.hh
short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.icc
Modified:
short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
Log:
Continued working on Drucker-Prager.
Everything is there except for computing elastoplastic elastic constants.
I still need to figure out how to deal with the non-symmetric constitutive
matrix (36 constants rather than 21). The present setup won't work with
the routines in feassemble.
Modified: short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.cc 2010-02-16 01:16:43 UTC (rev 16264)
+++ short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.cc 2010-02-16 04:56:43 UTC (rev 16265)
@@ -50,9 +50,9 @@
{ "density", 1, pylith::topology::FieldBase::SCALAR },
{ "mu", 1, pylith::topology::FieldBase::SCALAR },
{ "lambda", 1, pylith::topology::FieldBase::SCALAR },
- { "friction_angle", 1, pylith::topology::FieldBase::SCALAR },
- { "cohesion", 1, pylith::topology::FieldBase::SCALAR },
- { "dilatation_angle", 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
@@ -63,7 +63,7 @@
"dilatation-angle"};
/// Number of state variables.
- const int numStateVars = 2;
+ const int numStateVars = 1;
/// State variables.
const Metadata::ParamDescription stateVars[] = {
@@ -168,11 +168,11 @@
} else {
_calcStressFn =
- &pylith::materials::DruckerPragerEP3D::_calcStressViscoelastic;
+ &pylith::materials::DruckerPragerEP3D::_calcStressElastoplastic;
_calcElasticConstsFn =
- &pylith::materials::DruckerPragerEP3D::_calcElasticConstsViscoelastic;
+ &pylith::materials::DruckerPragerEP3D::_calcElasticConstsElastoplastic;
_updateStateVarsFn =
- &pylith::materials::DruckerPragerEP3D::_updateStateVarsViscoelastic;
+ &pylith::materials::DruckerPragerEP3D::_updateStateVarsElastoplastic;
} // if/else
} // useElasticBehavior
@@ -190,26 +190,34 @@
const double density = dbValues[db_density];
const double vs = dbValues[db_vs];
const double vp = dbValues[db_vp];
- const double referenceStrainRate = dbValues[db_referenceStrainRate];
- const double referenceStress = dbValues[db_referenceStress];
- const double powerLawExponent = dbValues[db_powerLawExponent];
+ const double frictionAngle = dbValues[db_frictionAngle];
+ const double cohesion = dbValues[db_cohesion];
+ const double dilatationAngle = dbValues[db_dilatationAngle];
- if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || referenceStrainRate <= 0.0
- || referenceStress <= 0.0 || powerLawExponent < 1.0) {
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || frictionAngle < 0.0
+ || cohesion <= 0.0 || dilatationAngle < 0.0
+ || frictionAngle < dilatationAngle) {
std::ostringstream msg;
msg << "Spatial database returned illegal value for physical "
<< "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n"
- << "referenceStrainRate: " << referenceStrainRate << "\n"
- << "referenceStress: " << referenceStress << "\n"
- << "powerLawExponent: " << powerLawExponent << "\n";
+ << "frictionAngle: " << frictionAngle << "\n"
+ << "cohesion: " << cohesion << "\n"
+ << "dilatationAngle: " << dilatationAngle << "\n";
throw std::runtime_error(msg.str());
} // if
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
+ const double alphaYield =
+ 2.0 * sin(frictionAngle)/(sqrt(3.0) * (3.0 - sin(frictionAngle)));
+ const double beta =
+ 6.0 * cohesion *
+ cos(frictionAngle)/(sqrt(3.0) * (3.0 - sin(frictionAngle)));
+ const double alphaFlow =
+ 2.0 * sin(dilatationAngle)/(sqrt(3.0) * (3.0 - sin(dilatationAngle)));
if (lambda <= 0.0) {
std::ostringstream msg;
@@ -224,11 +232,11 @@
propValues[p_density] = density;
propValues[p_mu] = mu;
propValues[p_lambda] = lambda;
- propValues[p_referenceStrainRate] = referenceStrainRate;
- propValues[p_referenceStress] = referenceStress;
- propValues[p_powerLawExponent] = powerLawExponent;
+ propValues[p_alphaYield] = alphaYield;
+ propValues[p_cohesion] = cohesion;
+ propValues[p_alphaFlow] = alphaFlow;
- PetscLogFlops(6);
+ PetscLogFlops(28);
} // _dbToProperties
// ----------------------------------------------------------------------
@@ -243,8 +251,6 @@
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
- const double timeScale = _normalizer->timeScale();
- const double strainRateScale = 1.0/timeScale;
values[p_density] =
_normalizer->nondimensionalize(values[p_density], densityScale);
@@ -252,13 +258,11 @@
_normalizer->nondimensionalize(values[p_mu], pressureScale);
values[p_lambda] =
_normalizer->nondimensionalize(values[p_lambda], pressureScale);
- values[p_referenceStrainRate] =
- _normalizer->nondimensionalize(values[p_referenceStrainRate],
- strainRateScale);
- values[p_referenceStress] =
- _normalizer->nondimensionalize(values[p_referenceStress], pressureScale);
+ values[p_beta] =
+ _normalizer->nondimensionalize(values[p_beta],
+ pressureScale);
- PetscLogFlops(6);
+ PetscLogFlops(4);
} // _nondimProperties
// ----------------------------------------------------------------------
@@ -273,8 +277,6 @@
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
- const double timeScale = _normalizer->timeScale();
- const double strainRateScale = 1.0/timeScale;
values[p_density] =
_normalizer->dimensionalize(values[p_density], densityScale);
@@ -282,12 +284,10 @@
_normalizer->dimensionalize(values[p_mu], pressureScale);
values[p_lambda] =
_normalizer->dimensionalize(values[p_lambda], pressureScale);
- values[p_referenceStrainRate] =
- _normalizer->dimensionalize(values[p_referenceStrainRate], strainRateScale);
- values[p_referenceStress] =
- _normalizer->dimensionalize(values[p_referenceStress], pressureScale);
+ values[p_beta] =
+ _normalizer->dimensionalize(values[p_beta], pressureScale);
- PetscLogFlops(6);
+ PetscLogFlops(4);
} // _dimProperties
// ----------------------------------------------------------------------
@@ -301,13 +301,11 @@
const int numDBValues = dbValues.size();
assert(_DruckerPragerEP3D::numDBStateVars == numDBValues);
- const int totalSize = 2 * _tensorSize;
+ const int totalSize = _tensorSize;
assert(totalSize == _numVarsQuadPt);
assert(totalSize == numDBValues);
- memcpy(&stateValues[s_viscousStrain], &dbValues[db_viscousStrain],
+ memcpy(&stateValues[s_plasticStrain], &dbValues[db_plasticStrain],
_tensorSize*sizeof(double));
- memcpy(&stateValues[s_stress], &dbValues[db_stress],
- _tensorSize*sizeof(double));
PetscLogFlops(0);
} // _dbToStateVars
@@ -322,10 +320,7 @@
assert(0 != values);
assert(nvalues == _numVarsQuadPt);
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->nondimensionalize(&values[s_stress], _tensorSize, pressureScale);
-
- PetscLogFlops(_tensorSize);
+ PetscLogFlops(0);
} // _nondimStateVars
// ----------------------------------------------------------------------
@@ -338,10 +333,7 @@
assert(0 != values);
assert(nvalues == _numVarsQuadPt);
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->dimensionalize(&values[s_stress], _tensorSize, pressureScale);
-
- PetscLogFlops(_tensorSize);
+ PetscLogFlops(0);
} // _dimStateVars
// ----------------------------------------------------------------------
@@ -373,38 +365,10 @@
assert(_numPropsQuadPt == numProperties);
assert(0 != stateVars);
assert(_numVarsQuadPt == numStateVars);
- const double mu = properties[p_mu];
- const double referenceStrainRate = properties[p_referenceStrainRate];
- const double referenceStress = properties[p_referenceStress];
- const double powerLawExp = properties[p_powerLawExponent];
-
- const double stress[] = {stateVars[s_stress],
- stateVars[s_stress + 1],
- stateVars[s_stress + 2],
- stateVars[s_stress + 3],
- stateVars[s_stress + 4],
- stateVars[s_stress + 5]};
- const double meanStress = (stress[0] + stress[1] + stress[2])/3.0;
- const double devStress[] = {stress[0] - meanStress,
- stress[1] - meanStress,
- stress[2] - meanStress,
- stress[3],
- stress[4],
- stress[5] };
- const double devStressProd = _scalarProduct(devStress, devStress);
- const double effStress = sqrt(0.5 * devStressProd);
- double dtTest = 1.0;
- if (effStress != 0.0)
- dtTest = 0.05 *
- pow((referenceStress/effStress), (powerLawExp - 1.0)) *
- (referenceStress/mu)/referenceStrainRate;
-
- const double dtStable = dtTest;
-#if 0 // DEBUGGING
- double maxwellTime = 10.0 * dtStable;
- std::cout << "Maxwell time: " << maxwellTime << std::endl;
-#endif
- PetscLogFlops(21);
+ // 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 double dtStable = 1.0e10;
+ PetscLogFlops(0);
return dtStable;
} // _stableTimeStepImplicit
@@ -465,10 +429,10 @@
} // _calcStressElastic
// ----------------------------------------------------------------------
-// Compute stress tensor at location from properties as a viscoelastic
+// Compute stress tensor at location from properties as an elastoplastic
// material.
void
-pylith::materials::DruckerPragerEP3D::_calcStressViscoelastic(
+pylith::materials::DruckerPragerEP3D::_calcStressElastoplastic(
double* const stress,
const int stressSize,
const double* properties,
@@ -482,7 +446,7 @@
const double* initialStrain,
const int initialStrainSize,
const bool computeStateVars)
-{ // _calcStressViscoelastic
+{ // _calcStressElastoplastic
assert(0 != stress);
assert(_DruckerPragerEP3D::tensorSize == stressSize);
assert(0 != properties);
@@ -497,42 +461,39 @@
assert(_DruckerPragerEP3D::tensorSize == initialStrainSize);
const int tensorSize = _tensorSize;
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
- // We need to do root-finding method if state variables are from previous
- // time step.
+ // We need to compute the plastic strain increment if state variables are
+ // from previous time step.
if (computeStateVars) {
- const double mu = properties[p_mu];
- const double lambda = properties[p_lambda];
- const double referenceStrainRate = properties[p_referenceStrainRate];
- const double referenceStress = properties[p_referenceStress];
- const double powerLawExp = properties[p_powerLawExponent];
- const double visStrainT[] = {stateVars[s_viscousStrain],
- stateVars[s_viscousStrain + 1],
- stateVars[s_viscousStrain + 2],
- stateVars[s_viscousStrain + 3],
- stateVars[s_viscousStrain + 4],
- stateVars[s_viscousStrain + 5]};
- const double stressT[] = {stateVars[s_stress],
- stateVars[s_stress + 1],
- stateVars[s_stress + 2],
- stateVars[s_stress + 3],
- stateVars[s_stress + 4],
- stateVars[s_stress + 5]};
-
+ const double alphaYield = properties[p_alphaYield];
+ const double beta = properties[p_beta];
+ const double alphaFlow = properties[p_alphaFlow];
const double mu2 = 2.0 * mu;
- const double lamPlusMu = lambda + mu;
const double bulkModulus = lambda + mu2/3.0;
const double ae = 1.0/mu2;
- const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+ const double am = 1.0/(3.0 * bulkModulus);
- // Need to figure out how time integration parameter alpha is going to be
- // specified. It should probably be specified in the problem definition and
- // then used only by the material types that use it. For now we are setting
- // it to 0.5, which should probably be the default value.
- const double alpha = 0.5;
- const double timeFac = _dt * (1.0 - alpha);
+ const double plasticStrainT[] = {stateVars[s_plasticStrain],
+ stateVars[s_plasticStrain + 1],
+ stateVars[s_plasticStrain + 2],
+ stateVars[s_plasticStrain + 3],
+ stateVars[s_plasticStrain + 4],
+ stateVars[s_plasticStrain + 5]};
+ const double meanPlasticStrainT = (plasticStrainT[0] +
+ plasticStrainT[1] +
+ plasticStrainT[2])/3.0;
+ const double devPlasticStrainT[] = { plasticStrainT[0] - meanPlasticStrainT,
+ plasticStrainT[1] - meanPlasticStrainT,
+ plasticStrainT[2] - meanPlasticStrainT,
+ plasticStrainT[3],
+ plasticStrainT[4],
+ plasticStrainT[5]};
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
// Initial stress values
const double meanStressInitial = (initialStress[0] +
initialStress[1] +
@@ -543,215 +504,124 @@
initialStress[3],
initialStress[4],
initialStress[5] };
- const double stressInvar2Initial = 0.5 *
- _scalarProduct(devStressInitial, devStressInitial);
// Initial strain values
const double meanStrainInitial = (initialStrain[0] +
initialStrain[1] +
initialStrain[2])/3.0;
+ const double devStrainInitial[] = { initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ initialStrain[2] - meanStrainInitial,
+ initialStrain[3],
+ initialStrain[4],
+ initialStrain[5] };
// Values for current time step
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
- const double meanStrainTpdt = (e11 + e22 + e33)/3.0 - meanStrainInitial;
- const double meanStressTpdt = 3.0 * bulkModulus * meanStrainTpdt;
+ const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
+ const double meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT -
+ meanStrainInitial;
- // Note that I use the initial strain rather than the deviatoric initial
- // strain since otherwise the initial mean strain would get used twice.
const double strainPPTpdt[] =
- { totalStrain[0] - meanStrainTpdt - visStrainT[0] - initialStrain[0],
- totalStrain[1] - meanStrainTpdt - visStrainT[1] - initialStrain[1],
- totalStrain[2] - meanStrainTpdt - visStrainT[2] - initialStrain[2],
- totalStrain[3] - visStrainT[3] - initialStrain[3],
- totalStrain[4] - visStrainT[4] - initialStrain[4],
- totalStrain[5] - visStrainT[5] - initialStrain[5] };
- const double strainPPInvar2Tpdt = 0.5 *
- _scalarProduct(strainPPTpdt, strainPPTpdt);
+ { 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] };
- // Values for previous time step
- const double meanStressT = (stressT[0] +
- stressT[1] +
- stressT[2])/3.0;
- const double devStressT[] = { stressT[0] - meanStressT,
- stressT[1] - meanStressT,
- stressT[2] - meanStressT,
- stressT[3],
- stressT[4],
- stressT[5] };
- const double stressInvar2T = 0.5 * _scalarProduct(devStressT, devStressT);
- const double effStressT = sqrt(stressInvar2T);
+ // Compute trial elastic stresses and yield function to see if yield should
+ // occur.
+ const double trialDevStress[] = { strainPPTpdt[0]/ae + devStressInitial[0],
+ strainPPTpdt[1]/ae + devStressInitial[1],
+ strainPPTpdt[2]/ae + devStressInitial[2],
+ strainPPTpdt[3]/ae + devStressInitial[3],
+ strainPPTpdt[4]/ae + devStressInitial[4],
+ strainPPTpdt[5]/ae + devStressInitial[5]};
+ const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
+ const double yieldFunction = 3.0* alphaYield * trialMeanStress +
+ _scalarProduct(trialDevStress, trialDevStress) - beta;
+ PetscLogFlops(74);
- // Finish defining parameters needed for root-finding algorithm.
- const double b = strainPPInvar2Tpdt +
- ae * _scalarProduct(strainPPTpdt, devStressInitial) +
- ae * ae * stressInvar2Initial;
- const double c = (_scalarProduct(strainPPTpdt, devStressT) +
- ae * _scalarProduct(devStressT, devStressInitial)) *
- timeFac;
- const double d = timeFac * effStressT;
+ // If yield function is greater than zero, compute elastoplastic stress.
+ if (yieldFunction >= 0.0) {
+ const double devStressInitialProd =
+ _scalarProduct(devStressInitial, devStressInitial);
+ const double strainPPTpdtProd =
+ _scalarProduct(strainPPTpdt, strainPPTpdt);
+ const double d = sqrt(ae * ae * devStressInitialProd +
+ 2.0 * ae *
+ _scalarProduct(devStressInitial, strainPPTpdt) +
+ strainPPTpdtProd);
+ plasticMult = 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am +
+ d/(sqrt(2.0) * ae) - beta)/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ const double meanStressTpdt =
+ (meanStrainPPTpdt - plasticMult * alphaFlow)/am + meanStressInitial;
+ double deltaDevPlasticStrain = 0.0;
+ double devStressTpdt = 0.0;
+ 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];
+ stress[iComp] = devStressTpdt + diag[iComp] * meanStressTpdt;
+ } // for
- PetscLogFlops(92);
+ PetscLogFlops(62 + 11 * tensorSize);
- // If b, c, and d are all zero, then the effective stress is zero and we
- // don't need a root-finding algorithm. Otherwise, use the algorithm to
- // find the effective stress.
- double effStressTpdt = 0.0;
- if (b != 0.0 || c != 0.0 || d != 0.0) {
- const double stressScale = mu;
-
- // Put parameters into a struct and call root-finding algorithm.
- _effStressParams.ae = ae;
- _effStressParams.b = b;
- _effStressParams.c = c;
- _effStressParams.d = d;
- _effStressParams.alpha = alpha;
- _effStressParams.dt = _dt;
- _effStressParams.effStressT = effStressT;
- _effStressParams.powerLawExp = powerLawExp;
- _effStressParams.referenceStrainRate = referenceStrainRate;
- _effStressParams.referenceStress = referenceStress;
-
- const double effStressInitialGuess = effStressT;
-
- double effStressTpdt =
- EffectiveStress::calculate<DruckerPragerEP3D>(effStressInitialGuess,
- stressScale, this);
+ } else {
+ // No plastic strain.
+ const double meanStressTpdt = meanStrainPPTpdt/am + meanStressInitial;
+ stress[0] = strainPPTpdt[0]/ae + devStressInitial[0] + meanStressTpdt;
+ stress[1] = strainPPTpdt[1]/ae + devStressInitial[1] + meanStressTpdt;
+ stress[2] = strainPPTpdt[2]/ae + devStressInitial[2] + meanStressTpdt;
+ stress[3] = strainPPTpdt[3]/ae + devStressInitial[3];
+ stress[4] = strainPPTpdt[4]/ae + devStressInitial[4];
+ stress[5] = strainPPTpdt[5]/ae + devStressInitial[5];
} // if
- // Compute stresses from effective stress.
- const double effStressTau = (1.0 - alpha) * effStressT +
- alpha * effStressTpdt;
- const double gammaTau = referenceStrainRate *
- pow((effStressTau/referenceStress),
- (powerLawExp - 1.0))/referenceStress;
- const double factor1 = 1.0/(ae + alpha * _dt * gammaTau);
- const double factor2 = timeFac * gammaTau;
- double devStressTpdt = 0.0;
-
- for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStressTpdt = factor1 *
- (strainPPTpdt[iComp] - factor2 * devStressT[iComp] +
- ae * devStressInitial[iComp]);
- stress[iComp] = devStressTpdt + diag[iComp] *
- (meanStressTpdt + meanStressInitial);
- } // for
- PetscLogFlops(14 + 8 * tensorSize);
-
- // If state variables have already been updated, current stress is already
- // contained in stress.
+ // If state variables have already been updated, the plastic strain for the
+ // time step has already been computed.
} else {
- memcpy(&stress[0], &stateVars[s_stress], tensorSize * sizeof(double));
- } // else
+ const double mu2 = 2.0 * mu;
+ const double plasticStrainTpdt[] = {stateVars[s_plasticStrain],
+ stateVars[s_plasticStrain + 1],
+ stateVars[s_plasticStrain + 2],
+ stateVars[s_plasticStrain + 3],
+ stateVars[s_plasticStrain + 4],
+ stateVars[s_plasticStrain + 5]};
-} // _calcStressViscoelastic
+ const double e11 = totalStrain[0] - plasticStrainTpdt[0] - initialStrain[0];
+ const double e22 = totalStrain[1] - plasticStrainTpdt[1] - initialStrain[1];
+ const double e33 = totalStrain[2] - plasticStrainTpdt[2] - initialStrain[2];
+ const double e12 = totalStrain[3] - plasticStrainTpdt[3] - initialStrain[3];
+ const double e23 = totalStrain[4] - plasticStrainTpdt[4] - initialStrain[4];
+ const double e13 = totalStrain[5] - plasticStrainTpdt[5] - initialStrain[5];
-// ----------------------------------------------------------------------
-// Effective stress function that computes effective stress function only
-// (no derivative).
-double
-pylith::materials::DruckerPragerEP3D::effStressFunc(const double effStressTpdt)
-{ // effStressFunc
- const double ae = _effStressParams.ae;
- const double b = _effStressParams.b;
- const double c = _effStressParams.c;
- const double d = _effStressParams.d;
- const double alpha = _effStressParams.alpha;
- const double dt = _effStressParams.dt;
- const double effStressT = _effStressParams.effStressT;
- const double powerLawExp = _effStressParams.powerLawExp;
- const double referenceStrainRate = _effStressParams.referenceStrainRate;
- const double referenceStress = _effStressParams.referenceStress;
- const double factor1 = 1.0-alpha;
- const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- const double gammaTau = referenceStrainRate *
- pow((effStressTau/referenceStress), (powerLawExp - 1.0))/referenceStress;
- const double a = ae + alpha * dt * gammaTau;
- const double y = a * a * effStressTpdt * effStressTpdt - b +
- c * gammaTau - d * d * gammaTau * gammaTau;
+ const double traceStrainTpdt = e11 + e22 + e33;
+ const double s123 = lambda * traceStrainTpdt;
- PetscLogFlops(21);
+ stress[0] = s123 + mu2 * e11 + initialStress[0];
+ stress[1] = s123 + mu2 * e22 + initialStress[1];
+ stress[2] = s123 + mu2 * e33 + initialStress[2];
+ stress[3] = mu2 * e12 + initialStress[3];
+ stress[4] = mu2 * e23 + initialStress[4];
+ stress[5] = mu2 * e13 + initialStress[5];
- return y;
-} // effStressFunc
+ PetscLogFlops(31);
-// ----------------------------------------------------------------------
-// Effective stress function that computes effective stress function
-// derivative only (no function value).
-double
-pylith::materials::DruckerPragerEP3D::effStressDerivFunc(const double effStressTpdt)
-{ // effStressDFunc
- const double ae = _effStressParams.ae;
- const double c = _effStressParams.c;
- const double d = _effStressParams.d;
- const double alpha = _effStressParams.alpha;
- const double dt = _effStressParams.dt;
- const double effStressT = _effStressParams.effStressT;
- const double powerLawExp = _effStressParams.powerLawExp;
- const double referenceStrainRate = _effStressParams.referenceStrainRate;
- const double referenceStress = _effStressParams.referenceStress;
- const double factor1 = 1.0-alpha;
- const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- const double gammaTau = referenceStrainRate *
- pow((effStressTau/referenceStress), (powerLawExp - 1.0))/referenceStress;
- const double a = ae + alpha * dt * gammaTau;
- const double dGammaTau = referenceStrainRate * alpha * (powerLawExp - 1.0) *
- pow((effStressTau/referenceStress), (powerLawExp - 2.0))/
- (referenceStress * referenceStress);
- const double dy = 2.0 * a * a * effStressTpdt + dGammaTau *
- (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
- c - 2.0 * d * d * gammaTau);
- PetscLogFlops(36);
+ } // else
- return dy;
-} // effStressDFunc
+} // _calcStressElastoplastic
// ----------------------------------------------------------------------
-// Effective stress function that computes effective stress function
-// and derivative.
-void
-pylith::materials::DruckerPragerEP3D::effStressFuncDerivFunc(double* func,
- double* dfunc,
- const double effStressTpdt)
-{ // effStressFuncDFunc
- double y = *func;
- double dy = *dfunc;
-
- const double ae = _effStressParams.ae;
- const double b = _effStressParams.b;
- const double c = _effStressParams.c;
- const double d = _effStressParams.d;
- const double alpha = _effStressParams.alpha;
- const double dt = _effStressParams.dt;
- const double effStressT = _effStressParams.effStressT;
- const double powerLawExp = _effStressParams.powerLawExp;
- const double referenceStrainRate = _effStressParams.referenceStrainRate;
- const double referenceStress = _effStressParams.referenceStress;
- const double factor1 = 1.0-alpha;
- const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- const double gammaTau = referenceStrainRate *
- pow((effStressTau/referenceStress), (powerLawExp - 1.0))/referenceStress;
- const double dGammaTau = referenceStrainRate * alpha * (powerLawExp - 1.0) *
- pow((effStressTau/referenceStress), (powerLawExp - 2.0))/
- (referenceStress * referenceStress);
- const double a = ae + alpha * dt * gammaTau;
- y = a * a * effStressTpdt * effStressTpdt -
- b +
- c * gammaTau -
- d * d * gammaTau * gammaTau;
- dy = 2.0 * a * a * effStressTpdt +
- dGammaTau *
- (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
- c - 2.0 * d * d * gammaTau);
-
- *func = y;
- *dfunc = dy;
-
- PetscLogFlops(46);
-} // effStressFuncDFunc
-
-// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from properties.
void
pylith::materials::DruckerPragerEP3D::_calcElasticConstsElastic(
@@ -793,30 +663,45 @@
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = lambda2mu; // C2222
- elasticConsts[ 7] = lambda; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = lambda2mu; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = mu2; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = mu2; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = mu2; // C1313
+ elasticConsts[ 6] = lambda; // C2211
+ elasticConsts[ 7] = lambda2mu; // C2222
+ elasticConsts[ 8] = lambda; // C2233
+ elasticConsts[ 9] = 0; // C2212
+ elasticConsts[10] = 0; // C2223
+ elasticConsts[11] = 0; // C2213
+ elasticConsts[12] = lambda; // C3311
+ elasticConsts[13] = lambda; // C3322
+ elasticConsts[14] = lambda2mu; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = mu2; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = mu2; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = mu2; // C1313
PetscLogFlops(2);
} // _calcElasticConstsElastic
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from properties
-// as a viscoelastic material.
+// as an elastoplastic material.
void
-pylith::materials::DruckerPragerEP3D::_calcElasticConstsViscoelastic(
+pylith::materials::DruckerPragerEP3D::_calcElasticConstsElastoplastic(
double* const elasticConsts,
const int numElasticConsts,
const double* properties,
@@ -829,7 +714,7 @@
const int initialStressSize,
const double* initialStrain,
const int initialStrainSize)
-{ // _calcElasticConstsViscoelastic
+{ // _calcElasticConstsElastoplastic
assert(0 != elasticConsts);
assert(_DruckerPragerEP3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
@@ -1052,7 +937,7 @@
PetscLogFlops(114);
} // else
-} // _calcElasticConstsViscoelastic
+} // _calcElasticConstsElastoplastic
// ----------------------------------------------------------------------
// Update state variables.
@@ -1080,20 +965,8 @@
assert(0 != initialStrain);
assert(_DruckerPragerEP3D::tensorSize == initialStrainSize);
- const bool computeStateVars = true;
- double stress[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
- const int stressSize = strainSize;
- _calcStressElastic(stress, stressSize,
- properties, numProperties,
- stateVars, numStateVars,
- totalStrain, strainSize,
- initialStress, initialStressSize,
- initialStrain, initialStrainSize,
- computeStateVars);
-
for (int iComp=0; iComp < _tensorSize; ++iComp) {
- stateVars[s_viscousStrain+iComp] = 0.0;
- stateVars[s_stress+iComp] = stress[iComp];
+ stateVars[s_plasticStrain+iComp] = 0.0;
} // for
_needNewJacobian = true;
@@ -1102,7 +975,7 @@
// ----------------------------------------------------------------------
// Update state variables.
void
-pylith::materials::DruckerPragerEP3D::_updateStateVarsViscoelastic(
+pylith::materials::DruckerPragerEP3D::_updateStateVarsElastoplastic(
double* const stateVars,
const int numStateVars,
const double* properties,
@@ -1113,7 +986,7 @@
const int initialStressSize,
const double* initialStrain,
const int initialStrainSize)
-{ // _updateStateVarsViscoelastic
+{ // _updateStateVarsElastoplastic
assert(0 != stateVars);
assert(_numVarsQuadPt == numStateVars);
assert(0 != properties);
@@ -1127,43 +1000,41 @@
const int stressSize = _tensorSize;
- // For now, we are duplicating the functionality of _calcStressViscoelastic,
+ // 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 double mu = properties[p_mu];
const double lambda = properties[p_lambda];
- const double referenceStrainRate = properties[p_referenceStrainRate];
- const double referenceStress = properties[p_referenceStress];
- const double powerLawExp = properties[p_powerLawExponent];
-
- const double visStrainT[] = {stateVars[s_viscousStrain],
- stateVars[s_viscousStrain + 1],
- stateVars[s_viscousStrain + 2],
- stateVars[s_viscousStrain + 3],
- stateVars[s_viscousStrain + 4],
- stateVars[s_viscousStrain + 5]};
-
- const double stressT[] = {stateVars[s_stress],
- stateVars[s_stress + 1],
- stateVars[s_stress + 2],
- stateVars[s_stress + 3],
- stateVars[s_stress + 4],
- stateVars[s_stress + 5]};
-
+ const double alphaYield = properties[p_alphaYield];
+ const double beta = properties[p_beta];
+ const double alphaFlow = properties[p_alphaFlow];
const double mu2 = 2.0 * mu;
- const double lamPlusMu = lambda + mu;
const double bulkModulus = lambda + mu2/3.0;
const double ae = 1.0/mu2;
- const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+ const double am = 1.0/(3.0 * bulkModulus);
- // Need to figure out how time integration parameter alpha is going to be
- // specified. It should probably be specified in the problem definition and
- // then used only by the material types that use it. For now we are setting
- // it to 0.5, which should probably be the default value.
- const double alpha = 0.5;
- const double timeFac = _dt * (1.0 - alpha);
+ const double plasticStrainT[] = {stateVars[s_plasticStrain],
+ stateVars[s_plasticStrain + 1],
+ stateVars[s_plasticStrain + 2],
+ stateVars[s_plasticStrain + 3],
+ stateVars[s_plasticStrain + 4],
+ stateVars[s_plasticStrain + 5]};
+ const double meanPlasticStrainT = (plasticStrainT[0] +
+ plasticStrainT[1] +
+ plasticStrainT[2])/3.0;
+ const double devPlasticStrainT[] = { plasticStrainT[0] - meanPlasticStrainT,
+ plasticStrainT[1] - meanPlasticStrainT,
+ plasticStrainT[2] - meanPlasticStrainT,
+ plasticStrainT[3],
+ plasticStrainT[4],
+ plasticStrainT[5]};
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
// Initial stress values
- const double meanStressInitial = (initialStress[0] + initialStress[1] +
+ const double meanStressInitial = (initialStress[0] +
+ initialStress[1] +
initialStress[2])/3.0;
const double devStressInitial[] = { initialStress[0] - meanStressInitial,
initialStress[1] - meanStressInitial,
@@ -1171,106 +1042,81 @@
initialStress[3],
initialStress[4],
initialStress[5] };
- const double stressInvar2Initial = 0.5 *
- _scalarProduct(devStressInitial, devStressInitial);
// Initial strain values
- const double meanStrainInitial = (initialStrain[0] + initialStrain[1] +
+ const double meanStrainInitial = (initialStrain[0] +
+ initialStrain[1] +
initialStrain[2])/3.0;
-
+ const double devStrainInitial[] = { initialStrain[0] - meanStrainInitial,
+ initialStrain[1] - meanStrainInitial,
+ initialStrain[2] - meanStrainInitial,
+ initialStrain[3],
+ initialStrain[4],
+ initialStrain[5] };
+
// Values for current time step
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
- const double meanStrainTpdt = (e11 + e22 + e33)/3.0 - meanStrainInitial;
- const double meanStressTpdt = 3.0 * bulkModulus * meanStrainTpdt;
+ const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
+ const double meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT -
+ meanStrainInitial;
- // Note that I use the initial strain rather than the deviatoric initial
- // strain since otherwise the initial mean strain would get used twice.
const double strainPPTpdt[] =
- { totalStrain[0] - meanStrainTpdt - visStrainT[0] - initialStrain[0],
- totalStrain[1] - meanStrainTpdt - visStrainT[1] - initialStrain[1],
- totalStrain[2] - meanStrainTpdt - visStrainT[2] - initialStrain[2],
- totalStrain[3] - visStrainT[3] - initialStrain[3],
- totalStrain[4] - visStrainT[4] - initialStrain[4],
- totalStrain[5] - visStrainT[5] - initialStrain[5] };
- const double strainPPInvar2Tpdt = 0.5 *
- _scalarProduct(strainPPTpdt, strainPPTpdt);
+ { 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] };
- // Values for previous time step
- const double meanStressT = (stressT[0] + stressT[1] + stressT[2])/3.0;
- const double devStressT[] = { stressT[0] - meanStressT,
- stressT[1] - meanStressT,
- stressT[2] - meanStressT,
- stressT[3],
- stressT[4],
- stressT[5] };
- const double stressInvar2T = 0.5 * _scalarProduct(devStressT, devStressT);
- const double effStressT = sqrt(stressInvar2T);
+ // Compute trial elastic stresses and yield function to see if yield should
+ // occur.
+ const double trialDevStress[] = { strainPPTpdt[0]/ae + devStressInitial[0],
+ strainPPTpdt[1]/ae + devStressInitial[1],
+ strainPPTpdt[2]/ae + devStressInitial[2],
+ strainPPTpdt[3]/ae + devStressInitial[3],
+ strainPPTpdt[4]/ae + devStressInitial[4],
+ strainPPTpdt[5]/ae + devStressInitial[5]};
+ const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
+ const double yieldFunction = 3.0* alphaYield * trialMeanStress +
+ _scalarProduct(trialDevStress, trialDevStress) - beta;
+ PetscLogFlops(74);
- // Finish defining parameters needed for root-finding algorithm.
- const double b = strainPPInvar2Tpdt +
- ae * _scalarProduct(strainPPTpdt, devStressInitial) +
- ae * ae * stressInvar2Initial;
- const double c = (_scalarProduct(strainPPTpdt, devStressT) +
- ae * _scalarProduct(devStressT, devStressInitial)) *
- timeFac;
- const double d = timeFac * effStressT;
- PetscLogFlops(92);
+ // If yield function is greater than zero, compute plastic strains.
+ // Otherwise, plastic strains remain the same.
+ if (yieldFunction >= 0.0) {
+ const double devStressInitialProd =
+ _scalarProduct(devStressInitial, devStressInitial);
+ const double strainPPTpdtProd =
+ _scalarProduct(strainPPTpdt, strainPPTpdt);
+ const double d = sqrt(ae * ae * devStressInitialProd +
+ 2.0 * ae *
+ _scalarProduct(devStressInitial, strainPPTpdt) +
+ strainPPTpdtProd);
+ plasticMult = 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am +
+ d/(sqrt(2.0) * ae) - beta)/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ const double deltaMeanPlasticStrain = plasticMult * alphaFlow;
+ double deltaDevPlasticStrain = 0.0;
+ for (int iComp=0; iComp < tensorSize; ++iComp) {
+ deltaDevPlasticStrain = plasticMult *(strainPPTpdt[iComp] +
+ ae * devStressInitial[iComp])/
+ (sqrt(2.0) * d);
+ stateVars[s_plasticStrain+iComp] += deltaDevPlasticStrain +
+ diag[iComp] * deltaMeanPlasticStrain;
+ } // for
- // If b, c, and d are all zero, then the effective stress is zero and we
- // don't need a root-finding algorithm. Otherwise, use the algorithm to
- // find the effective stress.
- double effStressTpdt = 0.0;
- if (b != 0.0 || c != 0.0 || d != 0.0) {
- const double stressScale = mu;
+ PetscLogFlops(60 + 9 * tensorSize);
- // Put parameters into a struct and call root-finding algorithm.
- _effStressParams.ae = ae;
- _effStressParams.b = b;
- _effStressParams.c = c;
- _effStressParams.d = d;
- _effStressParams.alpha = alpha;
- _effStressParams.dt = _dt;
- _effStressParams.effStressT = effStressT;
- _effStressParams.powerLawExp = powerLawExp;
- _effStressParams.referenceStrainRate = referenceStrainRate;
- _effStressParams.referenceStress = referenceStress;
-
- const double effStressInitialGuess = effStressT;
-
- double effStressTpdt =
- EffectiveStress::calculate<DruckerPragerEP3D>(effStressInitialGuess,
- stressScale, this);
-
} // if
- // Compute stress and viscous strain and update appropriate state variables.
- const double effStressTau = (1.0 - alpha) * effStressT +
- alpha * effStressTpdt;
- const double gammaTau = referenceStrainRate *
- pow((effStressTau/referenceStress),
- (powerLawExp - 1.0))/referenceStress;
- const double factor1 = 1.0/(ae + alpha * _dt * gammaTau);
- const double factor2 = timeFac * gammaTau;
- double devStressTpdt = 0.0;
- double devStressTau = 0.0;
- double deltaVisStrain = 0.0;
-
- for (int iComp=0; iComp < _tensorSize; ++iComp) {
- devStressTpdt = factor1 *
- (strainPPTpdt[iComp] - factor2 * devStressT[iComp] +
- ae * devStressInitial[iComp]);
- stateVars[s_stress+iComp] = devStressTpdt + diag[iComp] *
- (meanStressTpdt + meanStressInitial);
- devStressTau = (1.0 - alpha) * devStressT[iComp] + alpha * devStressTpdt;
- stateVars[s_viscousStrain+iComp] += _dt * gammaTau * devStressTau;
- } // for
-
_needNewJacobian = true;
- PetscLogFlops(14 + _tensorSize * 15);
-} // _updateStateVarsViscoelastic
+} // _updateStateVarsElastoplastic
// ----------------------------------------------------------------------
// Compute scalar product of two tensors.
Added: short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.hh 2010-02-16 04:56:43 UTC (rev 16265)
@@ -0,0 +1,491 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/DruckerPragerEP3D.hh
+ *
+ * @brief 3-D, isotropic, Drucker-Prager elastic/perfectly plastic material.
+ */
+
+#if !defined(pylith_materials_druckerpragerep3d_hh)
+#define pylith_materials_druckerpragerep3d_hh
+
+// Include directives ---------------------------------------------------
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
+
+// Powerlaw3D -----------------------------------------------------------
+/** @brief 3-D, isotropic, 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::DruckerPragerEP3D : public ElasticMaterial
+{ // class DruckerPragerEP3D
+ friend class TestDruckerPragerEP3D; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor
+ DruckerPragerEP3D(void);
+
+ /// Destructor
+ ~DruckerPragerEP3D(void);
+
+ /** Set current time step.
+ *
+ * @param dt Current time step.
+ */
+ void timeStep(const double 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(double* const propValues,
+ const double_array& dbValues) const;
+
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _nondimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _dimProperties(double* 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(double* const stateValues,
+ const double_array& dbValues) const;
+
+ /** Nondimensionalize state variables..
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _nondimStateVars(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize state variables.
+ *
+ * @param values Array of state variables.
+ * @param nvalues Number of values.
+ */
+ void _dimStateVars(double* 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(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* 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(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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
+ */
+ double _stableTimeStepImplicit(const double* properties,
+ const int numProperties,
+ const double* 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(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ // PRIVATE TYPEDEFS ///////////////////////////////////////////////////
+private :
+
+ /// Member prototype for _calcStress()
+ typedef void (pylith::materials::DruckerPragerEP3D::*calcStress_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const bool);
+
+ /// Member prototype for _calcElasticConsts()
+ typedef void (pylith::materials::DruckerPragerEP3D::*calcElasticConsts_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int);
+
+ /// Member prototype for _updateStateVars()
+ typedef void (pylith::materials::DruckerPragerEP3D::*updateStateVars_fn_type)
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ 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(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** Compute scalar product, assuming vector form of a tensor.
+ *
+ * @param tensor1 First tensor.
+ * @param tensor2 Second tensor.
+ */
+ double _scalarProduct(const double* tensor1,
+ const double* tensor2) 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;
+
+ 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_plasticStrain;
+ static const int db_plasticStrain;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ DruckerPragerEP3D(const DruckerPragerEP3D&); ///< Not implemented
+ const DruckerPragerEP3D& operator=(const DruckerPragerEP3D&); ///< Not implemented
+
+}; // class DruckerPragerEP3D
+
+#include "DruckerPragerEP3D.icc" // inline methods
+
+#endif // pylith_materials_druckerpragerep3d_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/DruckerPragerEP3D.icc 2010-02-16 04:56:43 UTC (rev 16265)
@@ -0,0 +1,105 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_materials_druckerpragerep3d_hh)
+#error "DruckerPragerEP3D.icc can only be included from DruckerPragerEP3D.hh"
+#endif
+
+#include <cassert> // USES assert()
+#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+
+// Set current time step.
+inline
+void
+pylith::materials::DruckerPragerEP3D::timeStep(const double 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::DruckerPragerEP3D::_calcStress(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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::DruckerPragerEP3D::_calcElasticConsts(
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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::DruckerPragerEP3D::_updateStateVars(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* 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/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2010-02-16 01:16:43 UTC (rev 16264)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2010-02-16 04:56:43 UTC (rev 16265)
@@ -537,7 +537,6 @@
stateVars[s_stress + 5]};
const double mu2 = 2.0 * mu;
- const double lamPlusMu = lambda + mu;
const double bulkModulus = lambda + mu2/3.0;
const double ae = 1.0/mu2;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
@@ -882,7 +881,6 @@
stateVars[s_stress + 5]};
const double mu2 = 2.0 * mu;
- const double lamPlusMu = lambda + mu;
const double bulkModulus = lambda + mu2/3.0;
const double ae = 1.0/mu2;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
@@ -1166,7 +1164,6 @@
stateVars[s_stress + 5]};
const double mu2 = 2.0 * mu;
- const double lamPlusMu = lambda + mu;
const double bulkModulus = lambda + mu2/3.0;
const double ae = 1.0/mu2;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
More information about the CIG-COMMITS
mailing list