[Problem] Allow mu to be cut off (from below)

dune/tectonic/frictiondata.hh

@@ -9,6 +9,7 @@ struct FrictionData {
         a(parset.get<double>("a")),
         b(parset.get<double>("b")),
         mu0(parset.get<double>("mu0")),
+        mumin(parset.get<double>("mumin")),
         normalStress(normalStress) {}
 
   double const L;
@@ -16,6 +17,7 @@ struct FrictionData {
   double const a;
   double const b;
   double const mu0;
+  double const mumin;
   double const normalStress;
 };
 #endif

dune/tectonic/frictionpotential.hh

@@ -34,23 +34,15 @@ class FrictionPotential : public FrictionPotentialWrapper {
   FrictionPotential(double coefficient, FrictionData const &fd)
       : fd(fd), weightTimesNormalStress(coefficient * fd.normalStress) {}
 
-  double evaluate(double V) const {
-    assert(V >= 0);
-    if (V <= V_m)
-      return 0;
-
-    // V log(V/V_m) - V + V_m
-    return weightTimesNormalStress * fd.a * (V * std::log(V / V_m) - V + V_m);
-  }
-
   // log(V/V_m)  if V >= V_0
   // 0           otherwise
   double differential(double V) const {
-    assert(V >= 0);
-    if (V <= V_m)
-      return 0;
+    assert(V >= 0.0);
+    if (V <= V_cutoff) {
+      return fd.mumin;
+    }
 
-    return weightTimesNormalStress * fd.a * std::log(V / V_m);
+    return weightTimesNormalStress * fd.a * (std::log(V) - logV_m);
   }
 
   // 1/V        if V >  V_0
@@ -58,7 +50,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
   // 0          if V <  V_0
   double second_deriv(double V) const {
     assert(V >= 0);
-    if (V <= V_m)
+    if (V <= V_cutoff)
       return 0;
 
     return weightTimesNormalStress * fd.a / V;
@@ -67,7 +59,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
   double regularity(double V) const {
     assert(V >= 0);
     // TODO: Make this controllable
-    if (std::abs(V - V_m) < 1e-14)
+    if (std::abs(V - V_cutoff) < 1e-14)
       return std::numeric_limits<double>::infinity();
 
     return std::abs(second_deriv(V));
@@ -77,14 +69,16 @@ class FrictionPotential : public FrictionPotentialWrapper {
   // with K = -a log(V_m / V_0) = mu_0 + b [ alpha + log(V_0 / L) ]
   void updateState(double state) {
     // state is assumed to be logarithmic
-    V_m = fd.V0 *
-          std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a);
+    logV_m = std::log(fd.V0) +
+             (-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a);
+    V_cutoff = std::exp(logV_m + fd.mumin / fd.a);
   }
 
 private:
   FrictionData const &fd;
   double const weightTimesNormalStress;
-  double V_m;
+  double logV_m;
+  double V_cutoff; // velocity at which mu falls short of mumin
 };
 
 class TrivialFunction : public FrictionPotentialWrapper {

src/one-body-sample.parset

@@ -51,6 +51,7 @@ requiredResidual = 1e-12
 
 [boundary.friction]
 mu0 = 0.6
+mumin = 0.0
 a = 0.010
 b = 0.015
 V0 = 1e-6