diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index 9a284eaa23a6d70cb7683b5b4b55de143e8ebd0e..5bf0d30f7d2695db0124f9e895c3b913b83a6e07 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -27,6 +27,11 @@ class GlobalNonlinearity {
     }
   }
 
+  void addHessianIndices(Dune::MatrixIndexSet &indices) const {
+    for (size_t i = 0; i < indices.rows(); ++i)
+      indices.add(i, i);
+  }
+
   void addGradient(VectorType const &v, VectorType &gradient) const {
     for (size_t i = 0; i < v.size(); ++i) {
       auto res = restriction(i);
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 7f0bac45c144d1f982d45dca5d134e8cfc6bd1cb..83fbfbc22dbbdc3fead55199a330f7aec7640f7b 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -91,8 +91,7 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
     // construct sparsity pattern for linearization
     Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());
     indices.import(problem.A);
-    // TODO: handle smooth domain
-    // problem.phi.addHessianIndices(indices);
+    problem.phi.addHessianIndices(indices);
 
     // construct matrix from pattern and initialize it
     indices.exportIdx(linearization.A);