Fix up RuinaFunction

dune/tectonic/nicefunction.hh

@@ -41,44 +41,25 @@ class RuinaFunction : public NiceFunction {
         h(h),
         rho(exp(-(mu + compound_state) / a)) {}
 
-  /*
-    If mu and sigma_n denote the coefficient of friction and the
-    normal stress, respectively, this function is given by
-
-    1/eta sigma_n r(max(eta id, rho)) + a rho + mu
-
-    with the constants a and b from Ruina's model and
-
-    r(beta) = beta (a (log beta - 1) + mu)
-
-    as well as
-
-    rho = exp(-mu/a)
-  */
+  // TODO: untested
   void virtual evaluate(double const &x, double &y) const {
+    assert(false);
     double const arg = eta * x / h;
-    if (arg <= rho) {
-      y = -std::log(rho); // TODO: We can write this out explicitly
-      return;
-    }
-    y = arg * (std::log(arg) - 1) - arg * std::log(rho) + rho - std::log(rho);
+    double const expstar = (arg == 0) ? 0 : arg * std::log(arg) - arg;
+    double const gamma =
+        normalStress * (mu + a * std::log(eta) + compound_state);
+    y = a * normalStress * expstar + gamma * arg;
     y *= coefficient * h * normalStress / eta;
   }
 
-  /*
-    (leaving some terms aside): with s > rho
-
-    1/eta d/dx [ a * (s log s - s) + mu s ] where s = eta x
-    = 1/eta [ a * (log (eta x) * eta) + eta mu ]
-    = a * log(eta x) + mu
-  */
   double virtual leftDifferential(double s) const {
     double const arg = eta * s / h;
-    if (arg <= rho)
+    if (arg == 0)
       return 0;
 
-    return coefficient * normalStress *
-           (a * std::log(arg) + mu + compound_state);
+    double const gamma =
+        normalStress * (mu + a * std::log(eta) + compound_state);
+    return coefficient * (a * normalStress * std::log(arg) + gamma);
   }
 
   /* see above */
@@ -86,21 +67,12 @@ class RuinaFunction : public NiceFunction {
     return leftDifferential(s);
   }
 
-  /*
-    d/dx a * log(eta x) + mu
-    = a * 1/(eta x) * eta
-    = a/x
-  */
   double virtual second_deriv(double s) const {
-    // includes the case eta * s = rho for which there is no second derivative
-    if (eta * s / h <= rho)
-      return 0;
-    else
-      return coefficient * normalStress * a / s;
+    return coefficient * a * normalStress / s;
   }
 
   double virtual regularity(double s) const {
-    if (eta * s == rho)
+    if (s == 0)
       return std::numeric_limits<double>::infinity();
 
     return std::abs(second_deriv(s));