diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 8fd4671a8a5b3a894dc63be9b0dbb2fc3b702e63..02e33783cfc72b35c24da18ae455d91e0e713fa5 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -19,6 +19,7 @@
 template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
 public:
   typedef MyConvexProblemTypeTEMPLATE MyConvexProblemType;
+  typedef typename MyConvexProblemType::FunctionType FunctionType;
   typedef typename MyConvexProblemType::NonlinearityType NonlinearityType;
   typedef typename MyConvexProblemType::VectorType VectorType;
   typedef typename MyConvexProblemType::MatrixType MatrixType;
@@ -98,9 +99,9 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject {
       }
       assert(localA != NULL);
 
-      // FIXME: Hardcoding a fixed function here for now
-      Dune::TrivialFunction func;
-      Dune::MyNonlinearity<block_size> phi(func);
+      FunctionType f;
+      problem.newphi.restriction(m, f);
+      Dune::MyNonlinearity<block_size> const phi(f);
       Dune::SampleFunctional<block_size> localJ(*localA, localb, phi);
 
       LocalVectorType correction;
diff --git a/dune/tectonic/myconvexproblem.hh b/dune/tectonic/myconvexproblem.hh
index ecbd73ce9d58618fc6f81008ab805c513ef6b6b3..9b61b489d45d3901bdf6bb7a0cee6c87b8e8f3eb 100644
--- a/dune/tectonic/myconvexproblem.hh
+++ b/dune/tectonic/myconvexproblem.hh
@@ -12,11 +12,11 @@
    part
 */
 template <class NonlinearityTypeTEMPLATE, class MatrixTypeTEMPLATE,
-          class VectorTypeTEMPLATE>
+          class VectorTypeTEMPLATE, class FunctionTypeTEMPLATE>
 class MyConvexProblem {
 public:
   typedef NonlinearityTypeTEMPLATE NonlinearityType;
-
+  typedef FunctionTypeTEMPLATE FunctionType;
   typedef VectorTypeTEMPLATE VectorType;
   typedef MatrixTypeTEMPLATE MatrixType;
   typedef typename VectorType::block_type LocalVectorType;
@@ -31,12 +31,14 @@ class MyConvexProblem {
       \param f The linear functional
       \param u The solution vector
   */
-  MyConvexProblem(MatrixType const &A, NonlinearityType &phi,
-                  VectorType const &f, VectorType &u)
-      : A(A), phi(phi), f(f), u(u) {};
+  MyConvexProblem(MatrixType const &A,
+                  Dune::MyGlobalNonlinearity<block_size, FunctionType> &newphi,
+                  NonlinearityType &phi, VectorType const &f, VectorType &u)
+      : A(A), newphi(newphi), phi(phi), f(f), u(u) {};
 
   MatrixType const &A;
   NonlinearityType &phi;
+  Dune::MyGlobalNonlinearity<block_size, FunctionType> const &newphi;
 
   VectorType const &f;
 
diff --git a/dune/tectonic/myglobalnonlinearity.hh b/dune/tectonic/myglobalnonlinearity.hh
index dfa47e5ff41d5e004da815290e26bd6c8d3003e5..80f4c69d1f95cd8966f19808081ed4a58be6b406 100644
--- a/dune/tectonic/myglobalnonlinearity.hh
+++ b/dune/tectonic/myglobalnonlinearity.hh
@@ -18,7 +18,7 @@ template <int dim, class OuterFunctionType> class MyGlobalNonlinearity {
         normalStress(normalStress),
         nodalIntegrals(nodalIntegrals) {}
 
-  void restriction(int i, OuterFunctionType &f) {
+  void restriction(int i, OuterFunctionType &f) const {
     double coefficient = normalStress[i] * nodalIntegrals[i];
     // FIXME: Assume Gamma = id and h_{n+1} = 1 for now;
     // We then only have to evaluate (1 + F_xi)(|x|)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 3b648758e9b2dd4129aae170d703b89d48c6ab2d..3c533a902da3bd3c8bf9dbaf626d56823b251ab3 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -36,6 +36,7 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 
+#include <dune/tectonic/myglobalnonlinearity.hh>
 #include <dune/tectonic/myconvexproblem.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
@@ -184,6 +185,10 @@ int main() {
           neumannBoundaryAssembler, nodalIntegrals,
           true); // resize the output vector and zero all of its entries
 
+      // FIXME: We should be using S_F instead of S_N here
+      Dune::MyGlobalNonlinearity<dim, Dune::LinearFunction>
+      myGlobalNonlinearity(coefficientOfFriction, normalStress, nodalIntegrals);
+
       { // constant 2D function
         std::vector<SmallVector> b;
         SmallVector SampleVector;
@@ -205,10 +210,11 @@ int main() {
       typedef ZeroNonlinearity<SmallVector, SmallMatrix> NonlinearityType;
       NonlinearityType phi;
 
-      typedef MyConvexProblem<NonlinearityType, OperatorType, VectorType>
-      MyConvexProblemType;
+      typedef MyConvexProblem<NonlinearityType, OperatorType, VectorType,
+                              Dune::LinearFunction> MyConvexProblemType;
 
-      MyConvexProblemType myConvexProblem(stiffnessMatrix, phi, f, u1);
+      MyConvexProblemType myConvexProblem(stiffnessMatrix, myGlobalNonlinearity,
+                                          phi, f, u1);
 
       typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
       MyBlockProblemType myBlockProblem(myConvexProblem);