[Cleanup] Add a convenience function restrict()

dune/tectonic/myblockproblem.hh

@@ -104,11 +104,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
 
     v /= vnorm; // Rescale for numerical stability
 
-    MyDirectionalConvexFunction<GlobalFriction<MatrixType, VectorType>> const
-        psi(computeDirectionalA(problem_.A, v),
-            computeDirectionalb(problem_.A, problem_.f, u, v), problem_.phi, u,
-            v);
-
+    auto const psi = restrict(problem_.A, problem_.f, u, v, problem_.phi);
     Dune::Solvers::Interval<double> D;
     psi.subDiff(0, D);
     // NOTE: Numerical instability can actually get us here

dune/tectonic/mydirectionalconvexfunction.hh

@@ -7,25 +7,6 @@
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/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 Matrix, class Vector>
-double computeDirectionalA(Matrix const &A, Vector const &v) {
-  return Arithmetic::Axy(A, v, v);
-}
-
-template <class Matrix, class Vector>
-double computeDirectionalb(Matrix const &A, Vector const &b, Vector const &u,
-                           Vector const &v) {
-  return Arithmetic::bmAxy(A, b, u, v);
-}
-
 template <class Nonlinearity> class MyDirectionalConvexFunction {
 public:
   using Vector = typename Nonlinearity::VectorType;
@@ -65,4 +46,23 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
 
   Dune::Solvers::Interval<double> dom;
 };
+
+/*
+  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 Matrix, class Vector, class Nonlinearity>
+MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
+                                                   Vector const &b,
+                                                   Vector const &u,
+                                                   Vector const &v,
+                                                   Nonlinearity const &phi) {
+  return MyDirectionalConvexFunction<Nonlinearity>(
+      Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
+}
+
 #endif