From ca2153a21af1c8d69f20bff60fb601a3e48ad384 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 26 Mar 2015 15:45:05 +0100
Subject: [PATCH] [Cleanup] Add a convenience function restrict()

---
 dune/tectonic/myblockproblem.hh              |  6 +---
 dune/tectonic/mydirectionalconvexfunction.hh | 38 ++++++++++----------
 2 files changed, 20 insertions(+), 24 deletions(-)

diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 2091f025..6de650c0 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/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
diff --git a/dune/tectonic/mydirectionalconvexfunction.hh b/dune/tectonic/mydirectionalconvexfunction.hh
index 38c73d46..11f20c52 100644
--- a/dune/tectonic/mydirectionalconvexfunction.hh
+++ b/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
-- 
GitLab