Implement addHessian()

dune/tectonic/globallaursennonlinearity.hh

@@ -13,8 +13,11 @@
 #include "nicefunction.hh"
 
 namespace Dune {
-template <int dim, class OuterFunctionType>
-class GlobalLaursenNonlinearity : public GlobalNonlinearity<dim> {
+template <int dim, class OuterFunctionType,
+          class VectorType = BlockVector<FieldVector<double, dim>>,
+          class MatrixType = BCRSMatrix<FieldMatrix<double, dim, dim>>>
+class GlobalLaursenNonlinearity
+    : public GlobalNonlinearity<dim, VectorType, MatrixType> {
 public:
   GlobalLaursenNonlinearity(
       shared_ptr<BlockVector<FieldVector<double, 1>> const> mu,

dune/tectonic/globalnonlinearity.hh

@@ -10,12 +10,22 @@
 #include "localnonlinearity.hh"
 
 namespace Dune {
-template <int dim> class GlobalNonlinearity {
+template <int dim, class VectorType = BlockVector<FieldVector<double, dim>>,
+          class MatrixType = BCRSMatrix<FieldMatrix<double, dim, dim>>>
+class GlobalNonlinearity {
 public:
   /*
     Return a restriction of the outer function to the i'th node.
   */
   virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const = 0;
+
+  virtual void addHessian(const VectorType& v, MatrixType& hessian) const {
+    // TODO: is this correct?
+    for (size_t i = 0; i < v.size(); ++i) {
+      auto res = restriction(i);
+      res->addHessian(v[i], hessian[i][i]);
+    }
+  }
 };
 }
 #endif

dune/tectonic/globalruinanonlinearity.hh

@@ -13,8 +13,10 @@
 #include "nicefunction.hh"
 
 namespace Dune {
-template <int dim>
-class GlobalRuinaNonlinearity : public GlobalNonlinearity<dim> {
+template <int dim, class VectorType = BlockVector<FieldVector<double, dim>>,
+          class MatrixType = BCRSMatrix<FieldMatrix<double, dim, dim>>>
+class GlobalRuinaNonlinearity
+    : public GlobalNonlinearity<dim, VectorType, MatrixType> {
 public:
   GlobalRuinaNonlinearity(
       shared_ptr<BlockVector<FieldVector<double, 1>> const> nodalIntegrals,

dune/tectonic/localnonlinearity.hh

 #ifndef DUNE_TECTONIC_LOCAL_NONLINEARITY_HH
 #define DUNE_TECTONIC_LOCAL_NONLINEARITY_HH
 
+#include <cmath>
 #include <limits>
 
 #include <dune/common/fvector.hh>
@@ -45,6 +46,38 @@ 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
+  */
+  void addHessian(VectorType const &x, MatrixType &A) const {
+    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);
+
+    for (int i = 0; i < dimension; ++i)
+      for (int j = 0; j < i; ++j) {
+        double const entry = (h2 - h1ox) * x[i] * x[j] / normX2;
+        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;
+  }
+
   void upperGradient(VectorType const &x, VectorType &ret) const {
     ret = x;
     ret *= func_->rightDifferential(x.two_norm()) / x.two_norm();

dune/tectonic/nicefunction.hh

@@ -14,6 +14,8 @@ class NiceFunction : public VirtualFunction<double, double> {
 
   virtual double leftDifferential(double s) const = 0;
   virtual double rightDifferential(double s) const = 0;
+
+  double virtual second_deriv(double s) const = 0;
 };
 
 class RuinaFunction : public NiceFunction {
@@ -110,6 +112,8 @@ class LinearFunction : public NiceFunction {
 
   double virtual rightDifferential(double s) const { return coefficient; }
 
+  double virtual second_deriv(double) const { return 0; }
+
 private:
   double coefficient;
 };
@@ -145,6 +149,8 @@ class TrivialFunction : public NiceFunction {
   double virtual leftDifferential(double) const { return 0; }
 
   double virtual rightDifferential(double) const { return 0; }
+
+  double virtual second_deriv(double) const { return 0; }
 };
 
 // slope in [n-1,n] is n