[Cleanup] Helpers for MyDirectionalConvexFunctional

dune/tectonic/minimisation.hh

@@ -7,7 +7,6 @@
 
 #include <dune/fufem/arithmetic.hh>
 #include <dune/fufem/interval.hh>
-#include <dune/solvers/computeenergy.hh>
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
 #include "mydirectionalconvexfunction.hh"
@@ -21,20 +20,9 @@ void descentMinimisation(Functional const &J,
   using SmallVector = typename Functional::SmallVector;
   using LocalNonlinearityType = typename Functional::NonlinearityType;
 
-  // {{{ Construct a restriction of J to the line x + t * v
-
-  /* We have
-
-     1/2 <A(x+tv),x+tv>-<b,x+tv> = 1/2 <Av,v> t^2 - <b-Ax,v> t + <1/2 Ax-b,x>
-
-     since A is symmetric.
-  */
-  SmallVector tmp = J.b;
-  Arithmetic::addProduct(tmp, -1.0, J.A, x);
-  double const JRestb = tmp * v;
-
   MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
-      2.0 * computeEnergy(J.A, v), JRestb, *J.phi, x, v);
+      computeDirectionalA(J.A, v), computeDirectionalb(J.A, J.b, x, v), *J.phi,
+      x, v);
   // }}}
 
   int count;

dune/tectonic/myblockproblem.hh

@@ -105,22 +105,10 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
     Arithmetic::addProduct(tmp, -1.0, problem.A, u);
     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<MatrixType, VectorType>> const psi(localA,
-                                                                    localb,
-                                                                    problem.phi,
-                                                                    u, v);
+        Dune::GlobalNonlinearity<MatrixType, VectorType>> const
+    psi(computeDirectionalA(problem.A, v),
+        computeDirectionalb(problem.A, problem.f, u, v), problem.phi, u, v);
 
     Interval<double> D;
     psi.subDiff(0, D);

dune/tectonic/mydirectionalconvexfunction.hh

@@ -9,6 +9,29 @@
 
 #include <dune/fufem/interval.hh>
 
+/*
+  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>
+*/
+
+template <class MatrixType, class VectorType>
+double computeDirectionalA(MatrixType const &A, VectorType const &v) {
+  VectorType tmp(v.size());
+  A.mv(v, tmp);
+  return tmp * v;
+}
+
+template <class MatrixType, class VectorType>
+double computeDirectionalb(MatrixType const &A, VectorType const &b,
+                           VectorType const &u, VectorType const &v) {
+  VectorType tmp = b;
+  Arithmetic::addProduct(tmp, -1.0, A, u);
+  return tmp * v;
+}
+
 template <class NonlinearityType> class MyDirectionalConvexFunction {
 public:
   using VectorType = typename NonlinearityType::VectorType;