[cig-commits] r15086 - in short/3D/PyLith/trunk/libsrc: . materials
brad at geodynamics.org
brad at geodynamics.org
Fri May 29 16:54:24 PDT 2009
Author: brad
Date: 2009-05-29 16:54:23 -0700 (Fri, 29 May 2009)
New Revision: 15086
Added:
short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
Removed:
short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
Log:
Use template functions in EffectiveStress to generalize interface.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2009-05-29 18:45:01 UTC (rev 15085)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2009-05-29 23:54:23 UTC (rev 15086)
@@ -77,7 +77,6 @@
materials/ViscoelasticMaxwell.cc \
materials/MaxwellIsotropic3D.cc \
materials/PowerLaw3D.cc \
- materials/EffectiveStress.cc \
meshio/BinaryIO.cc \
meshio/GMVFile.cc \
meshio/GMVFileAscii.cc \
Deleted: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc 2009-05-29 18:45:01 UTC (rev 15085)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc 2009-05-29 23:54:23 UTC (rev 15086)
@@ -1,203 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "EffectiveStress.hh" // implementation of object methods
-
-#include "petsc.h" // USES PetscLogFlops
-
-#include <cmath> // USES fabs()
-#include <cassert> // USES assert()
-#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 stressScale,
- const EffStressStruct& effStressParams,
- effStressFunc_fn_type effStressFunc,
- effStressFuncDFunc_fn_type effStressFuncDFunc)
-{ // getEffStress
- // Check parameters
- assert(effStressInitialGuess >= 0.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 = stressScale - 0.5 * stressScale;
- x2 = stressScale + 0.5 * stressScale;
- } // 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* px1,
- double* px2,
- const EffStressStruct& effStressParams,
- effStressFunc_fn_type 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 x1 = *px1;
- double x2 = *px2;
-
- double funcValue1 = effStressFunc(x1, effStressParams);
- double funcValue2 = effStressFunc(x2, effStressParams);
-
- int iteration = 0;
- bool bracketed = false;
- while (iteration < maxIterations) {
- if ((funcValue1 * funcValue2) < 0.0) {
- bracketed = true;
- break;
- } // if
-
- if (fabs(funcValue1) < fabs(funcValue2)) {
- 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
- ++iteration;
- } // while
-
- *px1 = x1;
- *px2 = x2;
-
- PetscLogFlops(5 * iteration);
- if (!bracketed)
- throw std::runtime_error("Unable to bracket effective stress.");
-} // _bracketEffStress
-
-// ----------------------------------------------------------------------
-// Find root using Newton's method with bisection.
-double
-pylith::materials::EffectiveStress::_findEffStress(
- const double x1,
- const double x2,
- const EffStressStruct& effStressParams,
- effStressFunc_fn_type effStressFunc,
- effStressFuncDFunc_fn_type 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-6;
-
- /// 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);
- assert(funcValueLow * funcValueHigh <= 0.0);
-
- double effStress = 0.0;
- double xLow = 0.0;
- double xHigh = 0.0;
- bool converged = false;
-
- if (fabs(funcValueLow) < accuracy) {
- effStress = x1;
- converged = true;
- return effStress;
- } else if (fabs(funcValueHigh) < accuracy) {
- effStress = x2;
- converged = true;
- return effStress;
- } else if (funcValueLow < 0.0) {
- xLow = x1;
- xHigh = x2;
- } else {
- xHigh = x1;
- xLow = x2;
- }
-
- effStress = 0.5 * (x1 + x2);
- double dxPrevious = fabs(x2 - x1);
- double dx = dxPrevious;
- double funcValue = 0.0;
- double funcDeriv = 0.0;
- double funcXHigh = 0.0;
- double funcXLow = 0.0;
- effStressFuncDFunc(effStress, effStressParams, &funcValue, &funcDeriv);
- int iteration = 0;
-
- while (iteration < maxIterations) {
- funcXHigh = (effStress - xHigh) * funcDeriv - funcValue;
- funcXLow = (effStress - xLow) * funcDeriv - funcValue;
- // Use bisection if solution goes out of bounds or is not converging
- // fast enough.
- if ( (funcXHigh * funcXLow >= 0.0) ||
- (fabs(2.0 * funcValue) > fabs(dxPrevious * funcDeriv))) {
- dxPrevious = dx;
- dx = 0.5 * (xHigh - xLow);
- effStress = xLow + dx;
- } else {
- dxPrevious = dx;
- dx = funcValue/funcDeriv;
- effStress = effStress - dx;
- } // else
- if (fabs(dx) < accuracy) {
- converged = true;
- break;
- } // if
- effStressFuncDFunc(effStress, effStressParams, &funcValue, &funcDeriv);
- if (funcValue < 0.0) {
- xLow = effStress;
- } else {
- xHigh = effStress;
- } // else
- ++iteration;
- } // while
-
- if (converged == false)
- throw std::runtime_error("Cannot find root of effective stress function.");
-
- PetscLogFlops(5 + 15 * iteration);
- return effStress;
-
- return effStress;
-} // _findEffStress
-
-
-// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh 2009-05-29 18:45:01 UTC (rev 15085)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh 2009-05-29 23:54:23 UTC (rev 15086)
@@ -14,8 +14,8 @@
*
* @brief C++ EffectiveStress object.
*
- * This class contains bracketing and root-finding functions for materials that
- * use an effective stress formulation.
+ * This class contains bracketing and root-finding functions for
+ * materials that use an effective stress formulation.
*/
#if !defined(pylith_materials_effectivestress_hh)
@@ -33,110 +33,61 @@
class pylith::materials::EffectiveStress
{ // class EffectiveStress
- // PUBLIC STRUCTS /////////////////////////////////////////////////////
-public :
-
- struct EffStressStruct {
- double ae;
- double b;
- double c;
- double d;
- double alpha;
- double dt;
- double effStressT;
- double powerLawExp;
- double viscosityCoeff;
- };
-
- // PUBLIC TYPEDEFS ///////////////////////////////////////////////////
-public :
-
- /// Member prototype for effStressFunc()
- typedef double (*effStressFunc_fn_type)
- (const double,
- const EffStressStruct&);
-
- /// Member prototype for effStressDFunc()
- typedef double (*effStressDFunc_fn_type)
- (const double,
- const EffStressStruct&);
-
- /// Member prototype for effStressFuncDFunc()
- typedef void (*effStressFuncDFunc_fn_type)
- (const double,
- const EffStressStruct&,
- double*,
- double*);
-
- // PUBLIC MEMBERS ///////////////////////////////////////////////////
-public :
-
- /// Metod to use for effStressFunc().
- effStressFunc_fn_type effStressFunc;
-
- /// Metod to use for effStressDFunc().
- effStressDFunc_fn_type effStressDFunc;
-
- /// Metod to use for effStressFuncDFunc().
- effStressFuncDFunc_fn_type effStressFuncDFunc;
-
-
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
/** Get effective stress from initial guess.
*
+ * The stressScale argument should provide a reasonable initial
+ * guess in the case where the
+ * actual initial guess is zero.
+ *
* @param effStressInitialGuess Initial guess for effective stress.
* @param effStressParams Parameters used in computing effective stress.
- * @param effStressFunc Function to compute effective stress only.
- * @param effStressFuncDFunc Function to compute effective stress and derivative.
+ * @param material Material with effective stress function.
*
* @returns Computed effective stress.
*/
+ template<typename material_type>
static
double getEffStress(const double effStressInitialGuess,
const double stressScale,
- const EffStressStruct& effStressParams,
- effStressFunc_fn_type effStressFunc,
- effStressFuncDFunc_fn_type effStressFuncDFunc);
+ material_type* const material);
// PRIVATE METHODS /////////////////////////////////////////////////////
private :
/** Bracket effective stress.
*
- * @param x1 Initial guess for first bracket.
- * @param x2 Initial guess for second bracket.
- * @param effStressParams Parameters used in computing effective stress.
- * @param effStressFunc Function to compute effective stress only.
+ * @param px1 Initial guess for first bracket.
+ * @param px2 Initial guess for second bracket.
+ * @param material Material with effective stress function.
*
*/
+ template<typename material_type>
static
void _bracketEffStress(double* px1,
double* px2,
- const EffStressStruct& effStressParams,
- effStressFunc_fn_type effStressFunc);
+ material_type* const material);
/** Solve for effective stress using Newton's method with bisection.
*
* @param x1 Initial guess for first bracket.
* @param x2 Initial guess for second bracket.
- * @param effStressParams Parameters used in computing effective stress.
- * @param effStressFunc Function to compute effective stress only.
- * @param effStressFuncDFunc Function to compute effective stress and derivative.
+ * @param material Material with effective stress function.
*
* @returns Computed effective stress.
*/
+ template<typename material_type>
static
double _findEffStress(double x1,
double x2,
- const EffStressStruct& effStressParams,
- effStressFunc_fn_type effStressFunc,
- effStressFuncDFunc_fn_type effStressFuncDFunc);
+ material_type* const material);
}; // class EffectiveStress
#endif // pylith_materials_effectivestress_hh
+#include "EffectiveStress.icc" // template methods
// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc 2009-05-29 23:54:23 UTC (rev 15086)
@@ -0,0 +1,191 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "petsc.h" // USES PetscLogFlops
+
+#include <cmath> // USES fabs()
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Get effective stress from initial guess.
+template<typename material_type>
+double
+pylith::materials::EffectiveStress::getEffStress(
+ const double effStressInitialGuess,
+ const double stressScale,
+ material_type* const material)
+{ // getEffStress
+ // Check parameters
+ assert(effStressInitialGuess >= 0.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 = stressScale - 0.5 * stressScale;
+ x2 = stressScale + 0.5 * stressScale;
+ } // else
+
+ _bracketEffStress(&x1, &x2, material);
+
+ // Find effective stress using Newton's method with bisection.
+ const double effStress = _findEffStress(x1, x2, material);
+
+ PetscLogFlops(4); // Log flops
+
+ return effStress;
+} // getEffStress
+
+// ----------------------------------------------------------------------
+// Bracket effective stress.
+template<typename material_type>
+void
+pylith::materials::EffectiveStress::_bracketEffStress(
+ double* px1,
+ double* px2,
+ material_type* const material)
+{ // _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 x1 = *px1;
+ double x2 = *px2;
+
+ double funcValue1 = material->effStressFunc(x1);
+ double funcValue2 = material->effStressFunc(x2);
+
+ int iteration = 0;
+ bool bracketed = false;
+ while (iteration < maxIterations) {
+ if ((funcValue1 * funcValue2) < 0.0) {
+ bracketed = true;
+ break;
+ } // if
+
+ if (fabs(funcValue1) < fabs(funcValue2)) {
+ x1 += bracketFactor * (x1 - x2);
+ x1 = std::max(x1, 0.0);
+ funcValue1 = material->effStressFunc(x1);
+ } else {
+ x2 += bracketFactor * (x1 - x2);
+ x2 = std::max(x2, 0.0);
+ funcValue2 = material->effStressFunc(x2);
+ } // else
+ ++iteration;
+ } // while
+
+ *px1 = x1;
+ *px2 = x2;
+
+ PetscLogFlops(5 * iteration);
+ if (!bracketed)
+ throw std::runtime_error("Unable to bracket effective stress.");
+} // _bracketEffStress
+
+// ----------------------------------------------------------------------
+// Find root using Newton's method with bisection.
+template<typename material_type>
+double
+pylith::materials::EffectiveStress::_findEffStress(
+ const double x1,
+ const double x2,
+ material_type* const material)
+{ // _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-6;
+
+ /// 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 = material->effStressFunc(x1);
+ double funcValueHigh = material->effStressFunc(x2);
+ assert(funcValueLow * funcValueHigh <= 0.0);
+
+ double effStress = 0.0;
+ double xLow = 0.0;
+ double xHigh = 0.0;
+ bool converged = false;
+
+ if (fabs(funcValueLow) < accuracy) {
+ effStress = x1;
+ converged = true;
+ return effStress;
+ } else if (fabs(funcValueHigh) < accuracy) {
+ effStress = x2;
+ converged = true;
+ return effStress;
+ } else if (funcValueLow < 0.0) {
+ xLow = x1;
+ xHigh = x2;
+ } else {
+ xHigh = x1;
+ xLow = x2;
+ }
+
+ effStress = 0.5 * (x1 + x2);
+ double dxPrevious = fabs(x2 - x1);
+ double dx = dxPrevious;
+ double funcValue = 0.0;
+ double funcDeriv = 0.0;
+ double funcXHigh = 0.0;
+ double funcXLow = 0.0;
+ material->effStressFuncDerivFunc(&funcValue, &funcDeriv, effStress);
+ int iteration = 0;
+
+ while (iteration < maxIterations) {
+ funcXHigh = (effStress - xHigh) * funcDeriv - funcValue;
+ funcXLow = (effStress - xLow) * funcDeriv - funcValue;
+ // Use bisection if solution goes out of bounds or is not converging
+ // fast enough.
+ if ( (funcXHigh * funcXLow >= 0.0) ||
+ (fabs(2.0 * funcValue) > fabs(dxPrevious * funcDeriv))) {
+ dxPrevious = dx;
+ dx = 0.5 * (xHigh - xLow);
+ effStress = xLow + dx;
+ } else {
+ dxPrevious = dx;
+ dx = funcValue / funcDeriv;
+ effStress = effStress - dx;
+ } // else
+ if (fabs(dx) < accuracy) {
+ converged = true;
+ break;
+ } // if
+ material->effStressFuncDerivFunc(&funcValue, &funcDeriv, effStress);
+ if (funcValue < 0.0) {
+ xLow = effStress;
+ } else {
+ xHigh = effStress;
+ } // else
+ ++iteration;
+ } // while
+
+ if (converged == false)
+ throw std::runtime_error("Cannot find root of effective stress function.");
+
+ PetscLogFlops(5 + 15 * iteration); // Log flops
+
+ return effStress;
+} // _findEffStress
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-29 18:45:01 UTC (rev 15085)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-29 23:54:23 UTC (rev 15086)
@@ -15,6 +15,7 @@
#include "PowerLaw3D.hh" // implementation of object methods
#include "Metadata.hh" // USES Metadata
+#include "EffectiveStress.hh" // USES EffectiveStress
#include "pylith/utils/array.hh" // USES double_array
@@ -602,12 +603,8 @@
const double effStressInitialGuess = effStressT;
double effStressTpdt =
- pylith::materials::EffectiveStress::getEffStress(
- effStressInitialGuess,
- stressScale,
- _effStressParams,
- pylith::materials::PowerLaw3D::effStressFunc,
- pylith::materials::PowerLaw3D::effStressFuncDFunc);
+ EffectiveStress::getEffStress<PowerLaw3D>(effStressInitialGuess,
+ stressScale, this);
// Compute stresses from effective stress.
const double effStressTau = (1.0 - alpha) * effStressT +
@@ -639,19 +636,17 @@
// Effective stress function that computes effective stress function only
// (no derivative).
double
-pylith::materials::PowerLaw3D::effStressFunc(
- const double effStressTpdt,
- const EffectiveStress::EffStressStruct& effStressParams)
+pylith::materials::PowerLaw3D::effStressFunc(const double effStressTpdt)
{ // effStressFunc
- double ae = effStressParams.ae;
- double b = effStressParams.b;
- double c = effStressParams.c;
- double d = effStressParams.d;
- double alpha = effStressParams.alpha;
- double dt = effStressParams.dt;
- double effStressT = effStressParams.effStressT;
- double powerLawExp = effStressParams.powerLawExp;
- double viscosityCoeff = effStressParams.viscosityCoeff;
+ double ae = _effStressParams.ae;
+ double b = _effStressParams.b;
+ double c = _effStressParams.c;
+ double d = _effStressParams.d;
+ double alpha = _effStressParams.alpha;
+ double dt = _effStressParams.dt;
+ double effStressT = _effStressParams.effStressT;
+ double powerLawExp = _effStressParams.powerLawExp;
+ double viscosityCoeff = _effStressParams.viscosityCoeff;
double factor1 = 1.0-alpha;
double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
@@ -667,17 +662,16 @@
// Effective stress function that computes effective stress function
// derivative only (no function value).
double
-pylith::materials::PowerLaw3D::effStressDFunc(const double effStressTpdt,
- const EffectiveStress::EffStressStruct& effStressParams)
+pylith::materials::PowerLaw3D::effStressDerivFunc(const double effStressTpdt)
{ // effStressDFunc
- double ae = effStressParams.ae;
- double c = effStressParams.c;
- double d = effStressParams.d;
- double alpha = effStressParams.alpha;
- double dt = effStressParams.dt;
- double effStressT = effStressParams.effStressT;
- double powerLawExp = effStressParams.powerLawExp;
- double viscosityCoeff = effStressParams.viscosityCoeff;
+ double ae = _effStressParams.ae;
+ double c = _effStressParams.c;
+ double d = _effStressParams.d;
+ double alpha = _effStressParams.alpha;
+ double dt = _effStressParams.dt;
+ double effStressT = _effStressParams.effStressT;
+ double powerLawExp = _effStressParams.powerLawExp;
+ double viscosityCoeff = _effStressParams.viscosityCoeff;
double factor1 = 1.0-alpha;
double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
@@ -690,6 +684,7 @@
(2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
c - 2.0 * d * d * gammaTau);
PetscLogFlops(36);
+
return dy;
} // effStressDFunc
@@ -697,24 +692,22 @@
// Effective stress function that computes effective stress function
// and derivative.
void
-pylith::materials::PowerLaw3D::effStressFuncDFunc(
- const double effStressTpdt,
- const EffectiveStress::EffStressStruct& effStressParams,
- double* py,
- double* pdy)
+pylith::materials::PowerLaw3D::effStressFuncDerivFunc(double* func,
+ double* dfunc,
+ const double effStressTpdt)
{ // effStressFuncDFunc
- double y = *py;
- double dy = *pdy;
+ double y = *func;
+ double dy = *dfunc;
- double ae = effStressParams.ae;
- double b = effStressParams.b;
- double c = effStressParams.c;
- double d = effStressParams.d;
- double alpha = effStressParams.alpha;
- double dt = effStressParams.dt;
- double effStressT = effStressParams.effStressT;
- double powerLawExp = effStressParams.powerLawExp;
- double viscosityCoeff = effStressParams.viscosityCoeff;
+ double ae = _effStressParams.ae;
+ double b = _effStressParams.b;
+ double c = _effStressParams.c;
+ double d = _effStressParams.d;
+ double alpha = _effStressParams.alpha;
+ double dt = _effStressParams.dt;
+ double effStressT = _effStressParams.effStressT;
+ double powerLawExp = _effStressParams.powerLawExp;
+ double viscosityCoeff = _effStressParams.viscosityCoeff;
double factor1 = 1.0-alpha;
double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
@@ -732,8 +725,8 @@
(2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
c - 2.0 * d * d * gammaTau);
- *py = y;
- *pdy = dy;
+ *func = y;
+ *dfunc = dy;
PetscLogFlops(46);
} // effStressFuncDFunc
@@ -1022,12 +1015,8 @@
const double effStressInitialGuess = effStressT;
const double effStressTpdt =
- EffectiveStress::getEffStress(
- effStressInitialGuess,
- stressScale,
- _effStressParams,
- pylith::materials::PowerLaw3D::effStressFunc,
- pylith::materials::PowerLaw3D::effStressFuncDFunc);
+ EffectiveStress::getEffStress<PowerLaw3D>(effStressInitialGuess,
+ stressScale, this);
// Compute quantities at intermediate time tau used to compute values at
// end of time step.
@@ -1345,12 +1334,8 @@
const double effStressInitialGuess = effStressT;
double effStressTpdt =
- EffectiveStress::getEffStress(
- effStressInitialGuess,
- stressScale,
- _effStressParams,
- pylith::materials::PowerLaw3D::effStressFunc,
- pylith::materials::PowerLaw3D::effStressFuncDFunc);
+ EffectiveStress::getEffStress<PowerLaw3D>(effStressInitialGuess,
+ stressScale, this);
// Compute stress and viscous strain and update appropriate state variables.
const double effStressTau = (1.0 - alpha) * effStressT +
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-05-29 18:45:01 UTC (rev 15085)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-05-29 23:54:23 UTC (rev 15086)
@@ -27,7 +27,6 @@
#define pylith_materials_powerlaw3d_hh
#include "ElasticMaterial.hh" // ISA ElasticMaterial
-#include "EffectiveStress.hh" // Uses EffectiveStress
/// 3-D, isotropic, linear Maxwell viscoelastic material.
class pylith::materials::PowerLaw3D : public ElasticMaterial
@@ -58,38 +57,29 @@
/** Compute effective stress function.
*
* @param effStressTpdt Effective stress value.
- * @param effStressParams Effective stress parameters.
*
* @returns Effective stress function value.
*/
- static double effStressFunc(
- const double effStressTpdt,
- const EffectiveStress::EffStressStruct& effStressParams);
+ double effStressFunc(const double effStressTpdt);
/** Compute effective stress function derivative.
*
* @param effStressTpdt Effective stress value.
- * @param effStressParams Effective stress parameters.
*
* @returns Effective stress function derivative value.
*/
- static double effStressDFunc(
- const double effStressTpdt,
- const EffectiveStress::EffStressStruct& effStressParams);
+ double effStressDerivFunc(const double effStressTpdt);
/** Compute effective stress function and derivative.
*
+ * @param func Returned effective stress function value.
+ * @param dfunc Returned effective stress function derivative value.
* @param effStressTpdt Effective stress value.
- * @param effStressParams Effective stress parameters.
- * @param py Returned effective stress function value.
- * @param pdy Returned effective stress function derivative value.
*
*/
- static void effStressFuncDFunc(
- const double effStressTpdt,
- const EffectiveStress::EffStressStruct& effStressParams,
- double* py,
- double* pdy);
+ void effStressFuncDerivFunc(double* func,
+ double* dfunc,
+ const double effStressTpdt);
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -507,11 +497,26 @@
double _scalarProduct(const double* tensor1,
const double* tensor2) const;
+ // PRIVATE STRUCTS ////////////////////////////////////////////////////
+private :
+
+ struct EffStressStruct {
+ double ae;
+ double b;
+ double c;
+ double d;
+ double alpha;
+ double dt;
+ double effStressT;
+ double powerLawExp;
+ double viscosityCoeff;
+ };
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
/// Structure to hold parameters for effective stress computation.
- pylith::materials::EffectiveStress::EffStressStruct _effStressParams;
+ EffStressStruct _effStressParams;
/// Method to use for _calcElasticConsts().
calcElasticConsts_fn_type _calcElasticConstsFn;
@@ -522,9 +527,6 @@
/// Method to use for _updateStateVars().
updateStateVars_fn_type _updateStateVarsFn;
- // PRIVATE MEMBERS ////////////////////////////////////////////////////
-private :
-
static const int p_density;
static const int p_mu;
static const int p_lambda;
@@ -544,12 +546,9 @@
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
- /// Not implemented
- PowerLaw3D(const PowerLaw3D&);
+ PowerLaw3D(const PowerLaw3D&); ///< Not implemented
+ const PowerLaw3D& operator=(const PowerLaw3D&); ///< Not implemented
- /// Not implemented
- const PowerLaw3D& operator=(const PowerLaw3D&);
-
}; // class PowerLaw3D
#include "PowerLaw3D.icc" // inline methods
More information about the CIG-COMMITS
mailing list