diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 51a7c5f1a4003c99a4d2cb470b347a232c20b19d..afd98c30f4f167ddde4bfa0349f620d192c61190 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -15,6 +15,22 @@
 #include "mydirectionalconvexfunction.hh"
 #include "samplefunctional.hh"
 
+/* Just for debugging */
+template <int dim, class VectorType, class MatrixType>
+double computeEnergy(
+    MatrixType const &A, VectorType const &b,
+    Dune::GlobalNonlinearity<dim, VectorType, MatrixType> const &phi,
+    VectorType const &x) {
+  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;
+}
+
 /** \brief Base class for problems where each block can be solved with a
  * modified gradient method */
 template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
diff --git a/dune/tectonic/mydirectionalconvexfunction.hh b/dune/tectonic/mydirectionalconvexfunction.hh
index 91b278da5efe6c2dc14824d3d5507eea99bd9c74..70667fee41400884418706ed38c8cfebf89b4122 100644
--- a/dune/tectonic/mydirectionalconvexfunction.hh
+++ b/dune/tectonic/mydirectionalconvexfunction.hh
@@ -17,6 +17,13 @@ template <class NonlinearityType> class MyDirectionalConvexFunction {
     phi_.directionalDomain(u_, v_, dom_);
   }
 
+  /* Just for debugging */
+  double operator()(double x) const {
+    VectorType tmp = v_;
+    tmp *= x;
+    return (0.5 * A * x * x) - (b * x) + phi_(tmp);
+  }
+
   double quadraticPart() const { return A; }
 
   double linearPart() const { return b; }