Handle 0 * 1/x for x close to 0 in addHessian

Since 0 * nan = nan

dune/tectonic/localnonlinearity.hh

@@ -54,39 +54,67 @@ template <int dimension> class LocalNonlinearity {
     }
   }
 
-  /*
-    H''(|x|) x \otimes x / |x|^2
-    + H'(|x|) [|x|^2 id - x \otimes x] / |x|^3
-
-    For i == j, this is (writing k and using the notation from below)
-
-    h2 * x[k] * x[k] / normX2 + h1 [normX2 - x[k]^2]/normX^3
-    = h2 * x[k] * x[k] / normX2 + h1ox - h1 x[k]^2/normX^3
-    = (h2-h1ox) * x[k] * x[k] / normX2 + h1ox
-
-    For i != j, this is (again, using the notation from below)
-
-    h2 * x[i] * x[j] / normX2 - h1 * x[i] * x[j]/normX^3
-    = (h2 - h1ox) * x[i] * x[j] / normX2
+  /** Formula for the derivative:
+
+    \f{align*}{
+      \frac {d^2}{dz^2} H(|z|)
+      &= \frac d{dz} \left( H'(|z|) \otimes \frac z{|z|} \right)\\
+      &= H''(|z|) \frac z{|z|} \otimes \frac z{|z|}
+      + H'(|z|) \left( \frac {|z| \operatorname{id} - z \otimes z/|z|}{|z|^2}
+    \right)\\
+      &= \frac {H''(|z|)}{|z|^2} z \otimes z
+      + \frac {H'(|z|)}{|z|} \operatorname{id}
+      - \frac {H'(|z|)}{|z|^3} z \otimes z\\
+      &= \left( \frac {H''(|z|)}{|z|^2} - \frac {H'(|z|)}{|z|^3} \right) z
+    \otimes z
+      + \frac {H'(|z|)}{|z|} \operatorname{id}
+    \f}
   */
   void addHessian(VectorType const &x, MatrixType &A) const {
-    if (x == VectorType(0))
-      return;
-
     double const normX2 = x.two_norm2();
-    double const normX = sqrt(normX2);
-    double const h1ox = func->rightDifferential(normX) / normX;
-    double const h2 = func->second_deriv(normX);
+    double const normX = std::sqrt(normX2);
+    double const normX3 = normX * normX2;
+
+    double const H1 = func->rightDifferential(normX);
+    double const H2 = func->second_deriv(normX);
+
+    // TODO: potential optimisation: factor out (H1 / normX), get rid of normX3
+    double const weight1 = H2 / normX2;
+    double const weight2 = -H1 / normX3;
+    double const weight3 = H1 / normX;
+
+    // {{{ 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;
+    // }}}
 
     for (int i = 0; i < dimension; ++i)
       for (int j = 0; j < i; ++j) {
-        double const entry = (h2 - h1ox) * x[i] * x[j] / normX2;
+        double const entry = tensorweight * x[i] * x[j];
         A[i][j] += entry;
         A[j][i] += entry;
       }
 
-    for (int k = 0; k < dimension; ++k)
-      A[k][k] += (h2 - h1ox) * x[k] * x[k] / normX2 + h1ox;
+    for (int k = 0; k < dimension; ++k) {
+      double const entry = tensorweight * x[k] * x[k];
+      A[k][k] += entry + idweight;
+    }
   }
 
   void addGradient(VectorType const &x, VectorType &y) const {