From b508d875f66bbea1697e4e7d6c4d7f42fe4b76df Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 16 Apr 2013 18:03:57 +0200
Subject: [PATCH] Use computeEnergy()

---
 dune/tectonic/myblockproblem.hh              | 14 +++++---------
 dune/tectonic/myconvexproblem.hh             |  8 ++++----
 dune/tectonic/mydirectionalconvexfunction.hh |  5 ++++-
 3 files changed, 13 insertions(+), 14 deletions(-)

diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index bd7f96ed..707f6aed 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -7,6 +7,9 @@
 #include <dune/common/nullptr.hh>
 #include <dune/common/parametertree.hh>
 
+// Just for debugging
+#include "dune/solvers/computeenergy.hh"
+
 #include <dune/fufem/arithmetic.hh>
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
@@ -17,18 +20,11 @@
 #include "ellipticenergy.hh"
 
 /* Just for debugging */
-template <int dim, class MatrixType, class VectorType>
+template <class MatrixType, class VectorType>
 double computeEnergy(
     MatrixType const &A, VectorType const &x, VectorType const &b,
     Dune::GlobalNonlinearity<MatrixType, VectorType> const &phi) {
-  double ret;
-  VectorType tmp(x.size());
-  A.mv(x, tmp);          //            Ax
-  ret = 0.5 * (tmp * x); // ret = 1/2 <Ax,x>
-
-  ret -= b * x;  // ret = 1/2 <Ax,x> - <b,x>
-  ret += phi(x); // ret = 1/2 <Ax,x> - <b,x> + phi(x)
-  return ret;
+  return computeEnergy(A, x, b) + phi(x);
 }
 
 /** \brief Base class for problems where each block can be solved with a
diff --git a/dune/tectonic/myconvexproblem.hh b/dune/tectonic/myconvexproblem.hh
index b8917229..542ea00a 100644
--- a/dune/tectonic/myconvexproblem.hh
+++ b/dune/tectonic/myconvexproblem.hh
@@ -4,6 +4,9 @@
 #ifndef MY_CONVEX_PROBLEM_HH
 #define MY_CONVEX_PROBLEM_HH
 
+// Just for debugging
+#include "dune/solvers/computeenergy.hh"
+
 #include "globalnonlinearity.hh"
 
 template <class MatrixTypeTEMPLATE, class VectorTypeTEMPLATE>
@@ -29,10 +32,7 @@ class MyConvexProblem {
 
   /* Just for debugging */
   double operator()(VectorType const &x) const {
-    double ret = phi(x) - (f * x);
-    VectorType tmp(x.size());
-    A.mv(x, tmp);
-    return ret + 0.5 * (tmp * x);
+    return phi(x) + computeEnergy(A, x, f);
   }
 
   MatrixType const &A;
diff --git a/dune/tectonic/mydirectionalconvexfunction.hh b/dune/tectonic/mydirectionalconvexfunction.hh
index 4aa33f0f..154aa7a5 100644
--- a/dune/tectonic/mydirectionalconvexfunction.hh
+++ b/dune/tectonic/mydirectionalconvexfunction.hh
@@ -4,6 +4,9 @@
 #ifndef MY_DIRECTIONAL_CONVEX_FUNCTION_HH
 #define MY_DIRECTIONAL_CONVEX_FUNCTION_HH
 
+// just for debugging
+#include <dune/solvers/computeenergy.hh>
+
 #include <dune/fufem/interval.hh>
 
 template <class NonlinearityType> class MyDirectionalConvexFunction {
@@ -21,7 +24,7 @@ template <class NonlinearityType> class MyDirectionalConvexFunction {
   double operator()(double x) const {
     VectorType tmp = v;
     tmp *= x;
-    return (0.5 * A * x * x) - (b * x) + phi(tmp);
+    return computeEnergy(A, x, b) + phi(tmp);
   }
 
   double quadraticPart() const { return A; }
-- 
GitLab