Use computeEnergy()

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

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;

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; }