Clean up RuinaFunction

Drop members a, V_0, K
Move a into coefficientProduct
Move V_0 into V_m

Add comments that show what we're computing

dune/tectonic/nicefunction.hh

@@ -38,56 +38,62 @@ class NiceFunction {
   }
 };
 
+// phi(V) = V log(V/V_m) - V + V_m  if V >= V_m
+//        = 0                       otherwise
+// with V_m = V_0 exp(-K/a),
+// i.e.   K = -a log(V_m / V_0) = mu_0 + b log(V_0 / L) + b alpha
 class RuinaFunction : public NiceFunction {
 public:
   RuinaFunction(double coefficient, double a, double mu_0, double V_0,
                 double normalStress, double b, double state, double L)
       : NiceFunction(),
-        a(a),
-        V_0(V_0),
-        coefficientProduct(coefficient * normalStress),
+        coefficientProduct(coefficient * a * normalStress),
         // state is assumed to be logarithmic
-        K(mu_0 + b * (state + std::log(V_0 / L))),
-        V_m(std::exp(-K / a)) {}
+        V_m(V_0 * std::exp(-(mu_0 + b * (state + std::log(V_0 / L))) / a))
+        // We could also compute V_m as
+        // V_0 * std::exp(-(mu_0 + b * state)/a)
+        //     * std::pow(V_0 / L, -b/a)
+        // which would avoid the std::exp(std::log())
+  {}
 
   double virtual evaluate(double V) const {
-    double const arg = V / V_0;
-    if (arg <= V_m)
+    assert(V >= 0);
+    if (V <= V_m)
       return 0;
 
-    return coefficientProduct *
-           (+a * arg * (std::log(arg) - 1) + K * arg + a * V_m);
+    // V log(V/V_m) - V + V_m
+    return coefficientProduct * (V * std::log(V / V_m) - V + V_m);
   }
 
+  // log(V/V_m)  if V >= V_0
+  // 0           otherwise
   double virtual leftDifferential(double V) const {
-    double const arg = V / V_0;
-    if (arg <= V_m)
+    assert(V >= 0);
+    if (V <= V_m)
       return 0;
 
-    double const ret = a * std::log(arg) + K;
-    return coefficientProduct * ret;
+    return coefficientProduct * std::log(V / V_m);
   }
 
   double virtual rightDifferential(double V) const {
     return leftDifferential(V);
   }
 
+  // 1/V        if V >  V_0
+  // undefined  if V == V_0
+  // 0          if V <  V_0
   double virtual second_deriv(double V) const {
     assert(V >= 0);
-    assert(V_0 > 0);
-    double const arg = V / V_0;
-    if (arg <= V_m)
+    if (V <= V_m)
       return 0;
 
-    return coefficientProduct * a / V;
+    return coefficientProduct / V;
   }
 
   double virtual regularity(double V) const {
     assert(V >= 0);
-    assert(V_0 > 0);
-    double const arg = V / V_0;
     // TODO: Make this controllable
-    if (std::abs(arg - V_m) < 1e-14)
+    if (std::abs(V - V_m) < 1e-14)
       return std::numeric_limits<double>::infinity();
 
     return std::abs(second_deriv(V));
@@ -96,10 +102,7 @@ class RuinaFunction : public NiceFunction {
   bool virtual smoothesNorm() const { return true; }
 
 private:
-  double const a;
-  double const V_0;
   double const coefficientProduct;
-  double const K;
   double const V_m;
 };