SmoothesNorm -> SmallestPositivePoint

dune/tectonic/frictionpotential.hh

@@ -18,15 +18,15 @@ class FrictionPotentialWrapper {
 
   double virtual differential(double s) const = 0;
   double virtual second_deriv(double x) const = 0;
-
   double virtual regularity(double s) const = 0;
 
-  // Whether H(|.|) is smooth at zero
-  bool virtual smoothesNorm() const { return false; }
-
   double virtual evaluate(double x) const {
     DUNE_THROW(NotImplemented, "evaluation not implemented");
   }
+
+  // Between 0 and this point, the function is constantly zero (and
+  // thus so are all derivatives)
+  double virtual smallestPositivePoint() const = 0;
 };
 
 // phi(V) = V log(V/V_m) - V + V_m  if V >= V_m
@@ -87,7 +87,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
     return std::abs(second_deriv(V));
   }
 
-  bool virtual smoothesNorm() const { return true; }
+  double virtual smallestPositivePoint() const { return V_m; }
 
 private:
   double const coefficientProduct;
@@ -104,7 +104,9 @@ class TrivialFunction : public FrictionPotentialWrapper {
 
   double virtual regularity(double) const { return 0; }
 
-  bool virtual smoothesNorm() const { return true; }
+  double virtual smallestPositivePoint() const {
+    return std::numeric_limits<double>::infinity();
+  }
 };
 }
 #endif

dune/tectonic/localfriction.hh

@@ -18,16 +18,16 @@ template <int dimension> class LocalFriction {
   using VectorType = FieldVector<double, dimension>;
   using MatrixType = FieldMatrix<double, dimension, dimension>;
 
-  LocalFriction(shared_ptr<FrictionPotentialWrapper const> func) : func(func) {}
+  LocalFriction(shared_ptr<FrictionPotentialWrapper const> func)
+      : func(func), smp(func->smallestPositivePoint()) {}
 
   double operator()(VectorType const &x) const {
     return func->evaluate(x.two_norm());
   }
 
   double regularity(VectorType const &x) const {
-    if (!func->smoothesNorm() && x.two_norm() < 1e-14) // TODO: Make this
-                                                       // controllable
-                                                       // (truncationRadius?)
+    double const xnorm = x.two_norm();
+    if (xnorm < 1e-14 && xnorm >= smp)
       return std::numeric_limits<double>::infinity();
 
     return func->regularity(x.two_norm());
@@ -37,12 +37,11 @@ template <int dimension> class LocalFriction {
   // u and v are assumed to be non-zero
   void directionalSubDiff(VectorType const &x, VectorType const &v,
                           Interval<double> &D) const {
-    if (x.two_norm() == 0) {
-      D[0] = D[1] = func->differential(0) * v.two_norm();
-      return;
-    }
-    double const un = x.two_norm();
-    D[0] = D[1] = func->differential(un) * (x * v) / un;
+    double const xnorm = x.two_norm();
+    if (xnorm < smp)
+      D[0] = D[1] = 0;
+    else
+      D[0] = D[1] = func->differential(xnorm) / xnorm * (x * v);
   }
 
   /** Formula for the derivative:
@@ -64,6 +63,8 @@ template <int dimension> class LocalFriction {
   void addHessian(VectorType const &x, MatrixType &A) const {
     double const xnorm2 = x.two_norm2();
     double const xnorm = std::sqrt(xnorm2);
+    if (xnorm < smp)
+      return;
     double const xnorm3 = xnorm * xnorm2;
 
     double const H1 = func->differential(xnorm);
@@ -110,7 +111,7 @@ template <int dimension> class LocalFriction {
 
   void addGradient(VectorType const &x, VectorType &y) const {
     double const xnorm = x.two_norm();
-    if (xnorm == 0 or std::isinf(func->regularity(xnorm)))
+    if (std::isinf(func->regularity(xnorm)))
       return;
 
     y.axpy(func->differential(xnorm) / xnorm, x);
@@ -129,6 +130,7 @@ template <int dimension> class LocalFriction {
 
 private:
   shared_ptr<FrictionPotentialWrapper const> const func;
+  double const smp;
 };
 }
 #endif