[cig-commits] r16337 - short/3D/PyLith/trunk/libsrc/materials

[willic3 at geodynamics.org]

Author: willic3
Date: 2010-02-24 20:03:53 -0800 (Wed, 24 Feb 2010)
New Revision: 16337

Modified:
   short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
Log:
Fixed problem where root could not be bracketed when initial guess was
very small and the root was very small.



Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc	2010-02-25 03:35:08 UTC (rev 16336)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc	2010-02-25 04:03:53 UTC (rev 16337)
@@ -29,11 +29,13 @@
 { // getEffStress
   // Check parameters
   assert(effStressInitialGuess >= 0.0);
+  // If initial guess is too low, use stress scale instead.
+  const double xMin = 1.0e-10;
 
   // Bracket the root.
   double x1 = 0.0;
   double x2 = 0.0;
-  if (effStressInitialGuess > 0.0) {
+  if (effStressInitialGuess > xMin) {
     x1 = effStressInitialGuess - 0.5 * effStressInitialGuess;
     x2 = effStressInitialGuess + 0.5 * effStressInitialGuess;
   } else {
@@ -65,10 +67,9 @@
   // Arbitrary factor by which to increase the brackets.
   const double bracketFactor = 2;
   // Minimum allowed value for effective stress.
-  const double xMin = -1.0e-10;
+  const double xMin = 0.0;
   double x1 = *px1;
   double x2 = *px2;
-  double xLow = 1.0e-10;
 
   double funcValue1 = material->effStressFunc(x1);
   double funcValue2 = material->effStressFunc(x2);
@@ -83,15 +84,13 @@
 
     if (fabs(funcValue1) < fabs(funcValue2)) {
       x1 += bracketFactor * (x1 - x2);
-      x1 = std::max(x1, xLow);
+      x1 = std::max(x1, xMin);
       funcValue1 = material->effStressFunc(x1);
     } else {
       x2 += bracketFactor * (x1 - x2);
-      x2 = std::max(x2, xLow);
+      x2 = std::max(x2, xMin);
       funcValue2 = material->effStressFunc(x2);
     } // else
-    if (x1 == xLow || x2 == xLow)
-      xLow *= 0.5;
     ++iteration;
   } // while