fix Tresca friction

dune/tectonic/data-structures/friction/frictionpotential.hh

@@ -34,14 +34,14 @@ class TruncatedRateState : public FrictionPotential {
       : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
 
   double coefficientOfFriction(double V) const override {
-    if (V <= Vmin or Vmin == 0.0)
+    if (V <= Vmin or regularity(V)>10e8)
       return 0.0;
 
     //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
 
     auto res = fd.a * std::log(V / Vmin);
     if (std::isinf(res)) {
-        std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
+        //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
     }
 
     return res;
@@ -49,16 +49,19 @@ class TruncatedRateState : public FrictionPotential {
 
   double differential(double V) const override {
     //std::cout << "differential: " <<  weight * fd.C - weightedNormalStress * coefficientOfFriction(V) << std::endl;
+    if (V <= Vmin or regularity(V)>10e8)
+      return 0.0;
+
     return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
   }
 
   double second_deriv(double V) const override {
     //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : -weightedNormalStress * (fd.a / V)) << std::endl;
 
-    if (V <= Vmin or Vmin == 0.0)
+    if (V <= Vmin)
       return 0.0;
 
-    return -weightedNormalStress * (fd.a / V);
+    return - weightedNormalStress * (fd.a / V);
   }
 
   double regularity(double V) const override {
@@ -71,10 +74,7 @@ class TruncatedRateState : public FrictionPotential {
   }
 
   double evaluate(double V) const override {
-
-    //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1)) << std::endl;
-
-    if (V <= Vmin or Vmin == 0.0)
+    if (V <= Vmin or regularity(V)>10e8)
         return 0.0;
 
       return weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1);
@@ -83,9 +83,7 @@ class TruncatedRateState : public FrictionPotential {
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
-
-    if (Vmin == 0.0)
-        std::cout << " logrest: " << logrest << " alpha: " << alpha << std::endl;
+    //Vmin = 0.0;
   }
 
 private:
@@ -143,6 +141,8 @@ class Tresca : public FrictionPotential {
   }
 
   double differential(double V) const override {
+    //if (std::isinf(regularity(V)))
+    //    return 0.0;
     return - weightedNormalStress * coefficientOfFriction(V);
   }
 
@@ -151,11 +151,14 @@ class Tresca : public FrictionPotential {
   }
 
   double regularity(double V) const override {
+    if (std::abs(V) < 1e-14) // TODO
+        return std::numeric_limits<double>::infinity();
+
     return std::abs(second_deriv(V));
   }
 
   double evaluate(double V) const override {
-      return 0;
+      return - weightedNormalStress * fd.mu0 * V;
   }
 
   void updateAlpha(double alpha) override {}