Make RuinaFunction convex(!)

dune/tectonic/nicefunction.hh

@@ -27,7 +27,9 @@ class NiceFunction {
   // Whether H(|.|) is smooth at zero
   bool virtual smoothesNorm() const { return false; }
 
-  void virtual evaluate(double const &x, double &y) const = 0;
+  void virtual evaluate(double const &x, double &y) const {
+    DUNE_THROW(NotImplemented, "evaluation not implemented");
+  }
 };
 
 class RuinaFunction : public NiceFunction {
@@ -35,30 +37,20 @@ class RuinaFunction : public NiceFunction {
   RuinaFunction(double coefficient, double a, double mu, double eta,
                 double normalStress, double b, double state, double L, double h)
       : a(a),
+        mu(mu),
+        eta(eta),
         h(h),
         coefficientProduct(coefficient * normalStress),
-        c(mu + (a - b) * std::log(eta) - b * std::log(L)),
-        K(b * state) // state is assumed to be logarithmic
-  {}
-
-  void virtual evaluate(double const &x, double &y) const {
-    double const arg = x / h;
-    if (arg == 0) { // TODO: Make this controllable
-      y = 0;
-      return;
-    }
-
-    double const expstar = arg * std::log(arg) - arg;
-    y = a * expstar + (c + K) * arg;
-    y *= coefficientProduct * h;
-  }
+        K(b *
+          (state - std::log(eta * L))), // state is assumed to be logarithmic
+        rho(std::exp(-(mu + K) / a)) {}
 
   double virtual leftDifferential(double x) const {
-    double const arg = x / h;
-    if (arg < 1e-14) // TODO: Make this controllable
+    double const arg = x / h * eta;
+    if (arg <= rho)
       return 0;
 
-    double const ret = (a * std::log(arg) + c + K);
+    double const ret = a * std::log(arg) + mu + K;
     return coefficientProduct * ret;
   }
 
@@ -66,26 +58,31 @@ class RuinaFunction : public NiceFunction {
     return leftDifferential(s);
   }
 
-  double virtual second_deriv(double s) const {
-    if (s < 1e-14) // TODO: Make this controllable
+  double virtual second_deriv(double x) const {
+    double const arg = x / h * eta;
+    if (arg <= rho)
       return 0;
 
-    return coefficientProduct * a / s;
+    return coefficientProduct * a / x;
   }
 
-  double virtual regularity(double s) const {
-    if (s == 0)
+  double virtual regularity(double x) const {
+    double const arg = x / h * eta;
+    // TODO: Make this controllable
+    if (arg < 1e-14 || std::abs(arg - rho) < 1e-14)
       return std::numeric_limits<double>::infinity();
 
-    return std::abs(second_deriv(s));
+    return std::abs(second_deriv(x));
   }
 
 private:
   double const a;
+  double const mu;
+  double const eta;
   double const h;
   double const coefficientProduct;
-  double const c;
   double const K;
+  double const rho;
 };
 
 class LinearFunction : public NiceFunction {