[cig-commits] r14953 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Sat May 9 18:33:18 PDT 2009
Author: willic3
Date: 2009-05-09 18:33:18 -0700 (Sat, 09 May 2009)
New Revision: 14953
Added:
short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
Modified:
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc
Log:
More work on PowerLaw3D, and started writing root-finding code for
effective stress function.
PowerLaw3D won't compile yet until I get functions passed correctly.
Effective stress function isn't complete yet.
Added: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc 2009-05-10 01:33:18 UTC (rev 14953)
@@ -0,0 +1,180 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "EffectiveStress.hh" // implementation of object methods
+
+#include "petsc.h" // USES PetscLogFlops
+
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Compute effective stress, given an initial guess, a vector of parameters,
+// and functions to compute the effective stress function and it's derivative.
+// The first entry in effStressParams should be a stress scaling factor
+// (e.g., mu) to provide a reasonable initial guess in the case where the
+// actual initial guess is zero.
+double
+pylith::materials::EffectiveStress::getEffStress(
+ const double effStressInitialGuess,
+ const double* effStressParams,
+ ***effStressFunc,
+ ***effStressFuncDFunc)
+{ // getEffStress
+ // Check parameters
+ if (effStressInitialGuess < 0.0)
+ throw std::runtime_error("Effective stress initial guess must be >= 0.");
+
+ // Bracket the root.
+ double x1 = 0.0;
+ double x2 = 0.0;
+ if (effStressInitialGuess > 0.0) {
+ x1 = effStressInitialGuess - 0.5 * effStressInitialGuess;
+ x2 = effStressInitialGuess + 0.5 * effStressInitialGuess;
+ } else {
+ x1 = effStressParams[0] - 0.5 * effStressParams[0];
+ x2 = effStressParams[0] + 0.5 * effStressParams[0];
+ } // else
+
+ PetscLogFlops(4);
+ _bracketEffStress(x1, x2, effStressParams, effStressFunc);
+
+ // Find effective stress using Newton's method with bisection.
+ const double effStress = _findEffStress(x1,
+ x2,
+ effStressParams,
+ effStressFunc,
+ effStressFuncDFunc);
+
+ return effStress;
+
+} // getEffStress
+
+// ----------------------------------------------------------------------
+// Bracket effective stress.
+void
+pylith::materials::EffectiveStress::_bracketEffStress(
+ double* const x1,
+ double* const x2,
+ const double* effStressParams,
+ ??? effStressFunc)
+{ // _bracketEffStress
+ // Arbitrary number of iterations to bracket the root
+ const int maxIterations = 50;
+
+ // Arbitrary factor by which to increase the brackets.
+ const double bracketFactor = 1.6;
+
+ double funcValue1 = effStressFunc(x1, effStressParams);
+ double funcValue2 = effStressFunc(x2, effStressParams);
+
+ int iteration = 0;
+ bool bracketed = false;
+
+ while ((iteration < maxIterations) && (bracketed == false)) {
+ if ((funcValue1 * funcValue2) < 0.0) {
+ bracketed = true;
+ } else {
+ if (std::abs(f1) < std::abs(f2)) {
+ x1 += bracketFactor * (x1 - x2);
+ x1 = std::max(x1, 0.0);
+ funcValue1 = effStressFunc(x1, effStressParams);
+ } else {
+ x2 += bracketFactor * (x1 - x2);
+ x2 = std::max(x2, 0.0);
+ funcValue2 = effStressFunc(x2, effStressParams);
+ } // else
+ } // else
+ ++iteration;
+ } // while
+
+ PetscLogFlops(5 * iteration);
+ if (bracketed == false)
+ throw std::runtime_error("Unable to bracket effective stress.");
+
+} // _bracketEffStress
+
+// ----------------------------------------------------------------------
+// Find root using Newton's method with bisection.
+void
+pylith::materials::EffectiveStress::_findEffStress(
+ double* const x1,
+ double* const x2,
+ const double* effStressParams,
+ ??? effStressFunc,
+ ??? effStressFuncDFunc)
+{ // _findEffStress
+ // Arbitrary number of iterations to find the root
+ const int maxIterations = 100;
+
+ // Desired accuracy for root. This is a bit arbitrary for now.
+ const double accuracy = 1.0e-10;
+
+ /// Determine if root has already been found, or if root is not bracketed.
+ // Otherwise, organize search so that effStressFunc(xLow) is less than zero.
+ double funcValueLow = effStressFunc(x1, effStressParams);
+ double funcValueHigh = effStressFunc(x2, effStressParams);
+ if (funcValueLow * funcValueHigh > 0.0)
+ throw std::runtime_error("Effective stress is not bracketed.");
+
+ double effStress = 0.0;
+ double xLow = 0.0;
+ double xHigh = 0.0;
+
+ if (std::abs(funcValueLow) < accuracy) {
+ effStress = x1;
+ return effStress;
+ } else if (std::abs(funcValueHigh) < accuracy) {
+ effStress = x2;
+ return effStress;
+ } else if (funcValueLow < 0.0) {
+ xLow = x1;
+ xHigh = x2;
+ } else {
+ xHigh = x1;
+ xLow = x2;
+ }
+
+ effStress = 0.5 * (x1 + x2);
+ double dxPrevious = std::abs(x2 - x1);
+ double dx = dxPrevious;
+ double funcValue = 0.0;
+ double funcDeriv = 0.0;
+ effStressFuncDFunc(effStress, effStressParams, funcValue, funcDeriv);
+ int iteration = 0;
+ bool converged = false;
+
+ // ******* finish fixing through here ***********
+ while ((iteration < maxIterations) && (converged == false)) {
+ if ((funcValue1 * funcValue2) < 0.0) {
+ bracketed = true;
+ } else {
+ if (std::abs(f1) < std::abs(f2)) {
+ x1 += bracketFactor * (x1 - x2);
+ x1 = std::max(x1, 0.0);
+ funcValue1 = effStressFunc(x1, effStressParams);
+ } else {
+ x2 += bracketFactor * (x1 - x2);
+ x2 = std::max(x2, 0.0);
+ funcValue2 = effStressFunc(x2, effStressParams);
+ } // else
+ } // else
+ } // while
+
+ if (bracketed == false)
+ throw std::runtime_error("Unable to bracket effective stress.");
+
+} // _bracketEffStress
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-10 00:16:38 UTC (rev 14952)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-10 01:33:18 UTC (rev 14953)
@@ -149,8 +149,6 @@
_updateStateVarsFn(0)
{ // constructor
useElasticBehavior(true);
- _viscousStrain.resize(_tensorSize);
- _stress.resize(_tensorSize);
} // constructor
// ----------------------------------------------------------------------
@@ -186,8 +184,8 @@
// Compute properties from values in spatial database.
void
pylith::materials::PowerLaw3D::_dbToProperties(
- double* const propValues,
- const double_array& dbValues) const
+ double* const propValues,
+ const double_array& dbValues) const
{ // _dbToProperties
assert(0 != propValues);
const int numDBValues = dbValues.size();
@@ -238,7 +236,7 @@
// Nondimensionalize properties.
void
pylith::materials::PowerLaw3D::_nondimProperties(double* const values,
- const int nvalues) const
+ const int nvalues) const
{ // _nondimProperties
assert(0 != _normalizer);
assert(0 != values);
@@ -272,7 +270,7 @@
{ // _dimProperties
assert(0 != _normalizer);
assert(0 != values);
- assert(nvalues == _totalPropsQuadPt);
+ assert(nvalues == _numPropsQuadPt);
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
@@ -317,6 +315,38 @@
} // _dbToStateVars
// ----------------------------------------------------------------------
+// Nondimensionalize state variables.
+void
+pylith::materials::PowerLaw3D::_nondimStateVars(double* const values,
+ const int nvalues) const
+{ // _nondimStateVars
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _numVarsQuadPt);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values[s_stress], _tensorSize, pressureScale);
+
+ PetscLogFlops(_tensorSize);
+} // _nondimStateVars
+
+// ----------------------------------------------------------------------
+// Dimensionalize state variables.
+void
+pylith::materials::PowerLaw3D::_dimStateVars(double* const values,
+ const int nvalues) const
+{ // _dimStateVars
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _numVarsQuadPt);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values[s_stress], _tensorSize, pressureScale);
+
+ PetscLogFlops(_tensorSize);
+} // _dimStateVars
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::PowerLaw3D::_calcDensity(double* const density,
@@ -347,7 +377,7 @@
assert(_numVarsQuadPt == numStateVars);
memcpy(&stress[0], &stateVars[s_stress],
- tensorSize * sizeof(double));
+ _tensorSize * sizeof(double));
const double meanStress = (stress[0] + stress[1] + stress[2])/3.0;
const double devStress[] = {stress[0] - meanStress,
stress[1] - meanStress,
@@ -361,6 +391,7 @@
dtStable = 0.1 * ((viscosityCoeff/effStressTpdt) ^ (powerLawExp - 1.0)) *
(viscosityCoeff/mu);
+ PetscLogFlops(20);
return dtStable;
} // _stableTimeStepImplicit
@@ -452,7 +483,7 @@
assert(0 != initialStrain);
assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const int tensorSize = _PowerLaw3D::tensorSize;
+ const int tensorSize = _tensorSize;
// We need to do root-finding method if state variables are from previous
// time step.
@@ -534,11 +565,13 @@
ae * _scalarProduct(devStressT, devStressInitial)) *
timeFac;
const double d = timeFac * effStressT;
+ const double stressScale = mu;
- PetscLogFlops(55);
+ PetscLogFlops(92);
// Put parameters into a vector and call root-finding algorithm.
// This could also be a struct.
- const double effStressParams[] = {ae,
+ const double effStressParams[] = {stressScale,
+ ae,
b,
c,
d,
@@ -551,10 +584,11 @@
// here. I would like the function to work something like:
const double effStressInitialGuess = effStressT;
double effStressTpdt =
- EffectiveStress::getEffStress(effStressInitialGuess,
- effStressParams,
- pylith::materials::PowerLaw3D::_effStressFunc,
- pylith::materials::PowerLaw3D::_effStressFuncDFunc);
+ EffectiveStress::getEffStress(
+ effStressInitialGuess,
+ effStressParams,
+ pylith::materials::PowerLaw3D::_effStressFunc,
+ pylith::materials::PowerLaw3D::_effStressFuncDFunc);
// Compute stresses from effective stress.
const double effStressTau = (1.0 - alpha) * effStressT +
@@ -586,9 +620,9 @@
// Effective stress function that computes effective stress function only
// (no derivative).
double
-pylith::materials::PowerLaw3D::_effStressFunc(double effStressTpdt,
- double *params)
-{ // _effStressFunc
+pylith::materials::PowerLaw3D::effStressFunc(double effStressTpdt,
+ double *params)
+{ // effStressFunc
double ae = params[0];
double b = params[1];
double c = params[2];
@@ -607,15 +641,15 @@
c * gammaTau - d * d * gammaTau * gammaTau;
PetscLogFlops(21);
return y;
-} // _effStressFunc
+} // effStressFunc
// ----------------------------------------------------------------------
// Effective stress function that computes effective stress function
// derivative only (no function value).
double
-pylith::materials::PowerLaw3D::_effStressDFunc(double effStressTpdt,
- double *params)
-{ // _effStressDFunc
+pylith::materials::PowerLaw3D::effStressDFunc(double effStressTpdt,
+ double *params)
+{ // effStressDFunc
double ae = params[0];
double c = params[2];
double d = params[3];
@@ -637,17 +671,17 @@
c - 2.0 * d * d * gammaTau);
PetscLogFlops(36);
return dy;
-} // _effStressDFunc
+} // effStressDFunc
// ----------------------------------------------------------------------
// Effective stress function that computes effective stress function
// and derivative.
void
-pylith::materials::PowerLaw3D::_effStressFuncDFunc(double effStressTpdt,
- double *params,
- double *y,
- double *dy)
-{ // _effStressFunc
+pylith::materials::PowerLaw3D::effStressFuncDFunc(double effStressTpdt,
+ double *params,
+ double *y,
+ double *dy)
+{ // effStressFunc
double ae = params[0];
double b = params[1];
double c = params[2];
@@ -671,7 +705,7 @@
(2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
c - 2.0 * d * d * gammaTau);
PetscLogFlops(46);
-} // _effStressFuncDFunc
+} // effStressFuncDFunc
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from properties.
@@ -766,7 +800,7 @@
assert(0 != initialStrain);
assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const int tensorSize = _PowerLaw3D::tensorSize;
+ const int tensorSize = _tensorSize;
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
const double viscosityCoeff = properties[p_viscosityCoeff];
@@ -811,7 +845,7 @@
elasticConsts[19] = 0; // C2313
elasticConsts[20] = elasticConsts[15]; // C1313
- PetscLogFlops(25);
+ PetscLogFlops(37);
} // _calcElasticConstsViscoelasticInitial
// ----------------------------------------------------------------------
@@ -846,7 +880,7 @@
assert(0 != initialStrain);
assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const int tensorSize = _PowerLaw3D::tensorSize;
+ const int tensorSize = _tensorSize;
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
@@ -871,20 +905,20 @@
/// Initial state.
// Initial stress values.
- const double meanStressInitial = (_stressInitial[0] + _stressInitial[1] +
- _stressInitial[2])/3.0;
- const double devStressInitial[] = { _stressInitial[0] - meanStressInitial,
- _stressInitial[1] - meanStressInitial,
- _stressInitial[2] - meanStressInitial,
- _stressInitial[3],
- _stressInitial[4],
- _stressInitial[5] };
+ const double meanStressInitial = (initialStress[0] + initialStress[1] +
+ initialStress[2])/3.0;
+ const double devStressInitial[] = { initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ initialStress[2] - meanStressInitial,
+ initialStress[3],
+ initialStress[4],
+ initialStress[5] };
const double stressInvar2Initial = 0.5 *
_scalarProduct(devStressInitial, devStressInitial);
// Initial strain values.
- const double meanStrainInitial = (_strainInitial[0] + _strainInitial[1] +
- _strainInitial[2])/3.0;
+ const double meanStrainInitial = (initialStrain[0] + initialStrain[1] +
+ initialStrain[2])/3.0;
/// Values for current time step
const double e11 = totalStrain[0];
@@ -907,13 +941,13 @@
_scalarProduct(strainPPTpdt, strainPPTpdt);
// 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 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);
@@ -925,11 +959,13 @@
ae * _scalarProduct(devStressT, devStressInitial)) *
timeFac;
const double d = timeFac * effStressT;
+ const double stressScale = mu;
- PetscLogFlops(45);
+ PetscLogFlops(92);
// Put parameters into a vector and call root-finding algorithm.
// This could also be a struct.
- const double effStressParams[] = {ae,
+ const double effStressParams[] = {stressScale,
+ ae,
b,
c,
d,
@@ -942,16 +978,16 @@
// here. I would like the function to work something like:
const double effStressInitialGuess = effStressT;
const double effStressTpdt =
- EffectiveStress::getEffStress(effStressInitialGuess,
- effStressParams,
- pylith::materials::PowerLaw3D::_effStressFunc,
- pylith::materials::PowerLaw3D::_effStressFuncDFunc);
+ EffectiveStress::getEffStress(
+ effStressInitialGuess,
+ effStressParams,
+ pylith::materials::PowerLaw3D::_effStressFunc,
+ pylith::materials::PowerLaw3D::_effStressFuncDFunc);
// Compute quantities at intermediate time tau used to compute values at
// end of time step.
const double effStressTau = (1.0 - alpha) * effStressT +
alpha * effStressTpdt;
-
const double gammaTau = 0.5 *
((effStressTau/viscosityCoeff)^(powerLawExp - 1.0))/viscosityCoeff;
const double factor1 = 1.0/(ae + alpha * _dt * gammaTau);
@@ -1048,128 +1084,242 @@
const double dStress13dStrain13 = factor1 * (1.0 - dGammaDStrain13 * dF13);
// Form elastic constants.
- elasticConsts[ 0] = bulkModulus + (2.0 * dStress11dStrain11 -
- dStress11dStrain22 -
- dStress11dStrain33)/3.0; // C1111
- elasticConsts[1] = bulkModulus - dStressDStrain/3.0; // C1122
- elasticConsts[2] = elasticConsts[2]; // C1133
+ elasticConsts[ 0] = bulkModulus + (2.0 * dStress11dStrain11
+ - dStress11dStrain22
+ - dStress11dStrain33)/3.0; // C1111
+ elasticConsts[ 1] = bulkModulus + (2.0 * dStress11dStrain22
+ - dStress11dStrain11
+ - dStress11dStrain33)/3.0; // C1122
+ elasticConsts[ 2] = bulkModulus + (2.0 * dStress11dStrain33
+ - dStress11dStrain22
+ - dStress11dStrain33)/3.0; // C1133
+ elasticConsts[ 3] = dStress11dStrain12; // C1112
+ elasticConsts[ 4] = dStress11dStrain23; // C1123
+ elasticConsts[ 5] = dStress11dStrain13; // C1113
+ elasticConsts[ 6] = bulkModulus + (2.0 * dStress22dStrain22
+ - dStress22dStrain11
+ - dStress22dStrain33)/3.0; // C2222
+ elasticConsts[ 7] = bulkModulus + (2.0 * dStress22dStrain33
+ - dStress22dStrain22
+ - dStress22dStrain33)/3.0; // C2233
+ elasticConsts[ 8] = dStress22dStrain12; // C2212
+ elasticConsts[ 9] = dStress22dStrain23; // C2223
+ elasticConsts[10] = dStress22dStrain13; // C2213
+ elasticConsts[11] = bulkModulus + (2.0 * dStress33dStrain33
+ - dStress33dStrain11
+ - dStress33dStrain22)/3.0; // C3333
+ elasticConsts[12] = dStress33dStrain12; // C3312
+ elasticConsts[13] = dStress33dStrain23; // C3323
+ elasticConsts[14] = dStress33dStrain13; // C3313
+ elasticConsts[15] = dStress12dStrain12; // C1212
+ elasticConsts[16] = dStress12dStrain23; // C1223
+ elasticConsts[17] = dStress12dStrain13; // C1213
+ elasticConsts[18] = dStress23dStrain23; // C2323
+ elasticConsts[19] = dStress23dStrain13; // C2313
+ elasticConsts[20] = dStress13dStrain13; // C1313
+ PetscLogFlops(265);
+} // _calcElasticConstsViscoelastic
- assert(0 != elasticConsts);
- assert(_PowerLaw3D::numElasticConsts == numElasticConsts);
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+pylith::materials::PowerLaw3D::_updateStateVarsElastic(
+ double* const stateVars,
+ const int numStateVars,
+ double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVarsElastic
+ assert(0 != stateVars);
+ assert(_totalPropsQuadPt == numStateVars);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_PowerLaw3D::tensorSize == strainSize);
- assert(_PowerLaw3D::tensorSize == initialStateSize);
-
- const double mu = properties[_PowerLaw3D::pidMu];
- const double lambda = properties[_PowerLaw3D::pidLambda];
- const double maxwelltime = properties[_PowerLaw3D::pidMaxwellTime];
+ assert(0 != initialStress);
+ assert(_PowerLaw3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const double mu2 = 2.0 * mu;
- const double bulkModulus = lambda + mu2/3.0;
+ 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);
- double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
+ for (int iComp=0; iComp < _tensorSize; ++iComp) {
+ stateVars[s_viscousStrain+iComp] = 0.0;
+ stateVars[s_stress+iComp] = stress[iComp];
+ } // for
- const double visFac = mu*dq/3.0;
- elasticConsts[ 0] = bulkModulus + 4.0*visFac; // C1111
- elasticConsts[ 1] = bulkModulus - 2.0*visFac; // C1122
- elasticConsts[ 2] = elasticConsts[1]; // C1133
- elasticConsts[ 3] = 0; // C1112
- elasticConsts[ 4] = 0; // C1123
- elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = elasticConsts[0]; // C2222
- elasticConsts[ 7] = elasticConsts[1]; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = elasticConsts[0]; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = 6.0 * visFac; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = elasticConsts[15]; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = elasticConsts[15]; // C1313
+ _needNewJacobian = true;
+} // _updateStateVarsElastic
- PetscLogFlops(10);
-} // _calcElasticConstsViscoelastic
-
// ----------------------------------------------------------------------
// Update state variables.
void
-pylith::materials::PowerLaw3D::_updatePropertiesElastic(
- double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _updatePropertiesElastic
+pylith::materials::PowerLaw3D::_updateStateVarsViscoelastic(
+ double* const stateVars,
+ const int numStateVars,
+ double* const properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVarsViscoelastic
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
assert(0 != totalStrain);
assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_PowerLaw3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const double maxwelltime = properties[_PowerLaw3D::pidMaxwellTime];
+ const int stressSize = _tensorSize;
+ // For now, we are duplicating the functionality of _calcStressViscoelastic,
+ // since otherwise we would have to redo a lot of calculations.
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+ const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double powerLawExp = properties[p_powerLawExponent];
+ memcpy(&visStrainT[0], &stateVars[s_viscousStrain],
+ tensorSize * sizeof(double));
+ memcpy(&stressT[0], &stateVars[s_stress],
+ tensorSize * sizeof(double));
+
+ 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 };
+
+ // 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);
+
+ // Initial stress values
+ const double meanStressInitial = (initialStress[0] + initialStress[1] +
+ initialStress[2])/3.0;
+ const double devStressInitial[] = { initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ initialStress[2] - meanStressInitial,
+ 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;
+
+ // 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 traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
+ // 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);
- const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+ // 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);
- for (int iComp=0; iComp < _PowerLaw3D::tensorSize; ++iComp) {
- properties[_PowerLaw3D::pidStrainT+iComp] = totalStrain[iComp];
- properties[_PowerLaw3D::pidVisStrainT+iComp] =
- totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
- } // for
- PetscLogFlops(3 + 2 * _PowerLaw3D::tensorSize);
+ // Finish defining parameters needed for root-finding algorithm.
+ const double b = strainInvar2Tpdt +
+ ae * _scalarProduct(strainPPTpdt, devStressInitial) +
+ ae * ae * stressInvar2Initial;
+ const double c = (_scalarProduct(strainPPTpdt, devStressT) +
+ ae * _scalarProduct(devStressT, devStressInitial)) *
+ timeFac;
+ const double d = timeFac * effStressT;
+ const double stressScale = mu;
- _needNewJacobian = true;
-} // _updatePropertiesElastic
+ PetscLogFlops(92);
+ // Put parameters into a vector and call root-finding algorithm.
+ // This could also be a struct.
+ const double effStressParams[] = {stressScale,
+ ae,
+ b,
+ c,
+ d,
+ alpha,
+ _dt,
+ effectiveStressT,
+ powerLawExp,
+ viscosityCoeff};
+ // I think the PETSc root-finding procedure is too involved for what we want
+ // here. I would like the function to work something like:
+ const double effStressInitialGuess = effStressT;
+ double effStressTpdt =
+ EffectiveStress::getEffStress(
+ effStressInitialGuess,
+ effStressParams,
+ pylith::materials::PowerLaw3D::_effStressFunc,
+ pylith::materials::PowerLaw3D::_effStressFuncDFunc);
-// ----------------------------------------------------------------------
-// Update state variables.
-void
-pylith::materials::PowerLaw3D::_updatePropertiesViscoelastic(
- double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize)
-{ // _updatePropertiesViscoelastic
- assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
- assert(0 != totalStrain);
- assert(_PowerLaw3D::tensorSize == strainSize);
- assert(_PowerLaw3D::tensorSize == initialStateSize);
+ // Compute stress and viscous strain and update appropriate state variables.
+ const double effStressTau = (1.0 - alpha) * effStressT +
+ alpha * effStressTpdt;
+ const double gammaTau = 0.5 *
+ ((effStressTau/viscosityCoeff)^(powerLawExp - 1.0))/viscosityCoeff;
+ const double factor1 = 1.0/(ae + alpha * _dt * gammaTau);
+ const double factor2 = timeFac * gammaTau;
+ double devStressTpdt = 0.0;
+ double devStressTau = 0.0;
- const int tensorSize = _PowerLaw3D::tensorSize;
+ 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
- pylith::materials::PowerLaw3D::_computeStateVars(properties,
- numProperties,
- totalStrain,
- strainSize,
- initialState,
- initialStateSize);
+ _needNewJacobian = true;
+ PetscLogFlops(14 + tensorSize * 14);
- memcpy(&properties[_PowerLaw3D::pidVisStrainT],
- &_visStrainT[0],
- tensorSize * sizeof(double));
- memcpy(&properties[_PowerLaw3D::pidStrainT],
- &totalStrain[0],
- tensorSize * sizeof(double));
-
- _needNewJacobian = false;
-
} // _updatePropertiesViscoelastic
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-05-10 00:16:38 UTC (rev 14952)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-05-10 01:33:18 UTC (rev 14953)
@@ -62,13 +62,6 @@
*/
void useElasticBehavior(const bool flag);
- /** Get flag indicating whether material implements an empty
- * _updateProperties() method.
- *
- * @returns False if _updateProperties() is empty, true otherwise.
- */
- bool usesUpdateProperties(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -99,20 +92,28 @@
void _dimProperties(double* const values,
const int nvalues) const;
- /** Nondimensionalize initial state.
+ /** Compute initial state variables from values in spatial database.
*
- * @param values Array of initial state values.
+ * @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 _nondimInitState(double* const values,
+ void _nondimStateVars(double* const values,
const int nvalues) const;
- /** Dimensionalize initial state.
+ /** Dimensionalize state variables.
*
- * @param values Array of initial state values.
+ * @param values Array of state variables.
* @param nvalues Number of values.
*/
- void _dimInitState(double* const values,
+ void _dimStateVars(double* const values,
const int nvalues) const;
/** Compute density from properties.
@@ -120,33 +121,46 @@
* @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 int numProperties,
+ const double* stateVars,
+ const int numStateVars);
- /** Compute stress tensor from properties. 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.
+ /** 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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties.
@@ -155,42 +169,65 @@
* @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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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* initialState,
- const int initialStateSize);
+ 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;
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const;
- /** Update properties (for next time step).
+ /** 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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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 _updateProperties(double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
+ 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 :
@@ -205,6 +242,10 @@
const int,
const double*,
const int,
+ const double*,
+ const int,
+ const double*,
+ const int,
const bool);
/// Member prototype for _calcElasticConsts()
@@ -216,56 +257,56 @@
const double*,
const int,
const double*,
+ const int,
+ const double*,
+ const int,
+ const double*,
const int);
- /// Member prototype for _updateProperties()
- typedef void (pylith::materials::PowerLaw3D::*updateProperties_fn_type)
+ /// Member prototype for _updateStateVars()
+ typedef void (pylith::materials::PowerLaw3D::*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 viscous strains (state variables) for the current time step.
- *
- * @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 _computeStateVars(const double* properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
-
/** 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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars);
/** Compute stress tensor from properties as an viscoelastic material.
@@ -274,20 +315,28 @@
* @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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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 _calcStressViscoelastic(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* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars);
/** Compute derivatives of elasticity matrix from properties as an
@@ -297,19 +346,27 @@
* @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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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* initialState,
- const int initialStateSize);
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
/** Compute derivatives of elasticity matrix from properties as a
* viscoelastic material.
@@ -318,35 +375,51 @@
* @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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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 _calcElasticConstsViscoelastic(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* initialState,
- const int initialStateSize);
+ 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 initialState Initial state values.
- * @param initialStateSize Size of initial state array.
+ * @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 _updatePropertiesElastic(double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
+ 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 a viscoelastic material.
*
@@ -357,37 +430,57 @@
* @param initialState Initial state values.
* @param initialStateSize Size of initial state array.
*/
- void _updatePropertiesViscoelastic(double* const properties,
- const int numProperties,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize);
+ void _updateStateVarsViscoelastic(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);
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
- /// Not implemented
- PowerLaw3D(const PowerLaw3D& m);
+ /// Method to use for _calcElasticConsts().
+ calcElasticConsts_fn_type _calcElasticConstsFn;
- /// Not implemented
- const PowerLaw3D& operator=(const PowerLaw3D& m);
+ /// Method to use for _calcStress().
+ calcStress_fn_type _calcStressFn;
+ /// Method to use for _updateStateVars().
+ updateProperties_fn_type _updateStateVarsFn;
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
- /// Viscous strain array.
- double_array _visStrain;
+ static const int p_density;
+ static const int p_mu;
+ static const int p_lambda;
+ static const int p_viscosityCoeff;
+ static const int p_powerLawExponent;
+ static const int db_density;
+ static const int db_vs;
+ static const int db_vp;
+ static const int db_viscosityCoeff;
+ static const int db_powerLawExponent;
- /// Method to use for _calcElasticConsts().
- calcElasticConsts_fn_type _calcElasticConstsFn;
+ static const int s_viscousStrain;
+ static const int s_stress;
+ static const int db_viscousStrain;
+ static const int db_stress;
- /// Method to use for _calcStress().
- calcStress_fn_type _calcStressFn;
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
- /// Method to use for _updateProperties().
- updateProperties_fn_type _updatePropertiesFn;
+ /// Not implemented
+ PowerLaw3D(const PowerLaw3D&);
+ /// Not implemented
+ const PowerLaw3D& operator=(const PowerLaw3D&);
+
}; // class PowerLaw3D
#include "PowerLaw3D.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc 2009-05-10 00:16:38 UTC (rev 14952)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc 2009-05-10 01:33:18 UTC (rev 14953)
@@ -14,7 +14,7 @@
#error "PowerLaw3D.icc can only be included from PowerLaw3D.hh"
#endif
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
// Set current time step.
@@ -27,57 +27,30 @@
_dt = dt;
} // timeStep
-// Set whether elastic or inelastic constitutive relations are used.
-inline
-void
-pylith::materials::PowerLaw3D::useElasticBehavior(const bool flag,
- const int iteration) {
- if (flag) {
- _calcStressFn =
- &pylith::materials::PowerLaw3D::_calcStressElastic;
- _calcElasticConstsFn =
- &pylith::materials::PowerLaw3D::_calcElasticConstsElastic;
- _updatePropertiesFn =
- &pylith::materials::PowerLaw3D::_updatePropertiesElastic;
- } else {
- _calcStressFn =
- &pylith::materials::PowerLaw3D::_calcStressViscoelastic;
- if (iteration == 0) {
- _calcElasticConstsFn =
- &pylith::materials::PowerLaw3D::_calcElasticConstsViscoelasticInitial;
- } else {
- _calcElasticConstsFn =
- &pylith::materials::PowerLaw3D::_calcElasticConstsViscoelastic;
- }
- _updatePropertiesFn =
- &pylith::materials::PowerLaw3D::_updatePropertiesViscoelastic;
- } // if/else
-} // useElasticBehavior
-
-// Get flag indicating whether material implements an empty
-inline
-bool
-pylith::materials::PowerLaw3D::usesUpdateProperties(void) const {
- return true;
-} // usesUpdateProperties
-
// Compute stress tensor from parameters.
inline
void
pylith::materials::PowerLaw3D::_calcStress(double* const stress,
- const int stressSize,
- const double* parameters,
- const int numParams,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize,
- const bool computeStateVars) {
+ 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,
- parameters, numParams,
+ properties, numProperties,
+ stateVars, numStateVars,
totalStrain, strainSize,
- initialState, initialStateSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize,
computeStateVars);
} // _calcStress
@@ -85,35 +58,49 @@
inline
void
pylith::materials::PowerLaw3D::_calcElasticConsts(
- double* const elasticConsts,
- const int numElasticConsts,
- const double* parameters,
- const int numParams,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize) {
+ 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,
- parameters, numParams,
+ properties, numProperties,
+ stateVars, numStateVars,
totalStrain, strainSize,
- initialState, initialStateSize);
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
} // _calcElasticConsts
// Update state variables after solve.
inline
void
-pylith::materials::PowerLaw3D::_updateProperties(double* const parameters,
- const int numParams,
- const double* totalStrain,
- const int strainSize,
- const double* initialState,
- const int initialStateSize) {
+pylith::materials::PowerLaw3D::_updateProperties(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 != _updatePropertiesFn);
- CALL_MEMBER_FN(*this, _updatePropertiesFn)(parameters, numParams,
+ CALL_MEMBER_FN(*this, _updatePropertiesFn)(stateVars, numStateVars,
+ properties, numProperties,
totalStrain, strainSize,
- initialState, initialStateSize);
-} // _updateProperties
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize);
+} // _updateStateVars
// Compute scalar product, assuming vector form of a tensor.
inline
@@ -127,7 +114,6 @@
2.0 * (tensor1[3] * tensor2[3] +
tensor1[4] * tensor2[4] +
tensor1[5] * tensor2[5]);
- PetscLogFlops(12);
return scalarProduct;
} // _scalarProduct
More information about the CIG-COMMITS
mailing list