[Algorit] Do away with the horrible near-zero mess

x/|x| isn't that big a deal and the logic was a gigantic mess

dune/tectonic/frictionpotential.hh

@@ -24,9 +24,6 @@ class FrictionPotentialWrapper {
     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;
   void virtual updateState(double) = 0;
 };
 
@@ -76,8 +73,6 @@ class FrictionPotential : public FrictionPotentialWrapper {
     return std::abs(second_deriv(V));
   }
 
-  double smallestPositivePoint() const { return V_m; }
-
   //    V_m = V_0 exp(-K/a),
   // with K = -a log(V_m / V_0) = mu_0 + b [ alpha + log(V_0 / L) ]
   void updateState(double state) {
@@ -102,10 +97,6 @@ class TrivialFunction : public FrictionPotentialWrapper {
 
   double regularity(double) const { return 0; }
 
-  double smallestPositivePoint() const {
-    return std::numeric_limits<double>::infinity();
-  }
-
   void updateState(double) {}
 };
 }

dune/tectonic/localfriction.hh

@@ -13,6 +13,23 @@
 
 #include "frictionpotential.hh"
 
+// In order to compute (x * y) / |x|, compute x/|x| first
+template <class VectorType>
+double dotFirstNormalised(VectorType const &x, VectorType const &y) {
+  double const xnorm = x.two_norm();
+  if (xnorm <= 0.0)
+    return 0.0;
+
+  size_t const xsize = x.size();
+  assert(xsize == y.size());
+
+  double sum = 0;
+  for (size_t i = 0; i < xsize; ++i)
+    sum += x[i] / xnorm * y[i];
+
+  return sum;
+}
+
 namespace Dune {
 template <int dimension> class LocalFriction {
 public:
@@ -20,20 +37,17 @@ template <int dimension> class LocalFriction {
   using MatrixType = FieldMatrix<double, dimension, dimension>;
 
   explicit LocalFriction(shared_ptr<FrictionPotentialWrapper> func)
-      : func(func), smp(func->smallestPositivePoint()) {}
+      : func(func) {}
 
   double operator()(VectorType const &x) const {
     return func->evaluate(x.two_norm());
   }
 
-  void updateState(double state) {
-    func->updateState(state);
-    smp = func->smallestPositivePoint();
-  }
+  void updateState(double state) { func->updateState(state); }
 
   double regularity(VectorType const &x) const {
     double const xnorm = x.two_norm();
-    if (xnorm < 1e-14 && xnorm >= smp)
+    if (xnorm <= 0.0)
       return std::numeric_limits<double>::infinity();
 
     return func->regularity(xnorm);
@@ -44,10 +58,10 @@ template <int dimension> class LocalFriction {
   void directionalSubDiff(VectorType const &x, VectorType const &v,
                           Interval<double> &D) const {
     double const xnorm = x.two_norm();
-    if (xnorm < smp)
-      D[0] = D[1] = 0;
+    if (xnorm <= 0.0)
+      D[0] = D[1] = 0.0;
     else
-      D[0] = D[1] = func->differential(xnorm) / xnorm * (x * v);
+      D[0] = D[1] = func->differential(xnorm) * dotFirstNormalised(x, v);
   }
 
   /** Formula for the derivative:
@@ -69,38 +83,14 @@ 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)
+    if (xnorm2 <= 0.0)
       return;
-    double const xnorm3 = xnorm * xnorm2;
 
     double const H1 = func->differential(xnorm);
     double const H2 = func->second_deriv(xnorm);
 
-    // TODO: potential optimisation: factor out (H1 / xnorm), get rid of xnorm3
-    double const weight1 = H2 / xnorm2;
-    double const weight2 = -H1 / xnorm3;
-    double const weight3 = H1 / xnorm;
-
-    // {{{ In what follows, we handle the case 0 * (1/x) = 0 with x
-    // close to 0.
-    double tensorweight = 0;
-
-    if (std::isnan(weight1))
-      assert(H2 == 0);
-    else
-      tensorweight += weight1;
-
-    if (std::isnan(weight2))
-      assert(H1 == 0);
-    else
-      tensorweight += weight2;
-
-    double idweight = 0;
-    if (std::isnan(weight3))
-      assert(H1 == 0);
-    else
-      idweight += weight3;
-    // }}}
+    double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
+    double const idweight = H1 / xnorm;
 
     for (int i = 0; i < dimension; ++i)
       for (int j = 0; j < i; ++j) {
@@ -120,7 +110,7 @@ template <int dimension> class LocalFriction {
     if (std::isinf(func->regularity(xnorm)))
       return;
 
-    if (xnorm >= smp)
+    if (xnorm > 0.0)
       Arithmetic::addProduct(y, func->differential(xnorm) / xnorm, x);
   }
 
@@ -132,7 +122,6 @@ template <int dimension> class LocalFriction {
 
 private:
   shared_ptr<FrictionPotentialWrapper> const func;
-  double smp;
 };
 }
 #endif