Implement computeDampingParameter

dune/tectonic/globalnonlinearity.hh

@@ -10,10 +10,14 @@
 #include "localnonlinearity.hh"
 
 namespace Dune {
-template <int dim, class VectorType = BlockVector<FieldVector<double, dim>>,
-          class MatrixType = BCRSMatrix<FieldMatrix<double, dim, dim>>>
+template <int dim,
+          class VectorTypeTEMPLATE = BlockVector<FieldVector<double, dim>>,
+          class MatrixTypeTEMPLATE = BCRSMatrix<FieldMatrix<double, dim, dim>>>
 class GlobalNonlinearity {
 public:
+  typedef VectorTypeTEMPLATE VectorType;
+  typedef MatrixTypeTEMPLATE MatrixType;
+
   /*
     Return a restriction of the outer function to the i'th node.
   */
@@ -27,6 +31,25 @@ class GlobalNonlinearity {
     }
   }
 
+  void directionalDomain(VectorType const &, VectorType const &,
+                         Interval<double> &dom) const {
+    dom[0] = -std::numeric_limits<double>::max();
+    dom[1] = std::numeric_limits<double>::max();
+  }
+
+  void directionalSubDiff(VectorType const &u, VectorType const &v,
+                          Interval<double> &subdifferential) const {
+    subdifferential[0] = subdifferential[1] = 0;
+    for (size_t i = 0; i < u.size(); ++i) {
+      Interval<double> D;
+      auto res = restriction(i);
+      // TODO: is the coordinate system the right one here?
+      res->directionalSubDiff(u[i], v[i], D);
+      subdifferential[0] += D[0];
+      subdifferential[1] += D[1];
+    }
+  }
+
   void addHessianIndices(Dune::MatrixIndexSet &indices) const {
     for (size_t i = 0; i < indices.rows(); ++i)
       indices.add(i, i);

dune/tectonic/myblockproblem.hh

@@ -10,6 +10,7 @@
 #include <dune/tnnmg/problem-classes/onedconvexfunction.hh>
 
 #include "localnonlinearity.hh"
+#include "mydirectionalconvexfunction.hh"
 #include "samplefunctional.hh"
 
 /** \brief Base class for problems where each block can be solved with a
@@ -71,9 +72,29 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
 
   double computeDampingParameter(VectorType const &u,
                                  VectorType const &projected_v) const {
-    // TODO: implement
-    return 1.0;
-  };
+    VectorType const v = projected_v;
+
+    VectorType tmp = problem.f;
+    problem.A.mmv(u, tmp);
+    double const localb = tmp * v;
+
+    problem.A.mv(v, tmp);
+    double const localA = tmp * v;
+
+    /*
+      1/2 <A(u + hv),u + hv> - <b, u + hv>
+      = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
+
+      localA = <Av,v>
+      localb = <b - Au, v>
+    */
+    MyDirectionalConvexFunction<Dune::GlobalNonlinearity<block_size>> psi(
+        localA, localb, problem.phi, u, v);
+
+    int bisectionsteps = 0;
+    Bisection bisection(0.0, 1.0, 1e-12, true, 0);            // TODO
+    return bisection.minimize(psi, 0.0, 0.0, bisectionsteps); // TODO
+  }
 
   void assembleTruncate(VectorType const &u, Linearization &linearization,
                         Dune::BitSetVector<block_size> const &ignore) const {