diff --git a/dune/gfe/trustregionsolver.cc b/dune/gfe/trustregionsolver.cc
index 7bd147319bde4ec8cc50a896b05ba0e4b04d203e..df3dfd334fd02f7fa174c642f700c07ee4b80165 100644
--- a/dune/gfe/trustregionsolver.cc
+++ b/dune/gfe/trustregionsolver.cc
@@ -37,7 +37,8 @@ setup(const GridType& grid,
          int nu1,
          int nu2,
          int baseIterations,
-         double baseTolerance)
+         double baseTolerance,
+         const SolutionType& pointLoads)
 {
     grid_                     = &grid;
     assembler_                = assembler;
@@ -48,6 +49,7 @@ setup(const GridType& grid,
     innerIterations_          = multigridIterations;
     innerTolerance_           = mgTolerance;
     ignoreNodes_              = &dirichletNodes;
+    pointLoads_               = pointLoads;
 
     int numLevels = grid_->maxLevel()+1;
 
@@ -204,7 +206,7 @@ void TrustRegionSolver<GridType,VectorType>::solve()
     //   Trust-Region Solver
     // /////////////////////////////////////////////////////
 
-    double oldEnergy = assembler_->computeEnergy(x_);
+    double oldEnergy = assembler_->computeEnergy(x_, pointLoads_);
 
     bool recomputeGradientHessian = true;
     CorrectionType rhs;
@@ -226,6 +228,7 @@ void TrustRegionSolver<GridType,VectorType>::solve()
         if (recomputeGradientHessian) {
 
             assembler_->assembleGradientAndHessian(x_,
+                                                   pointLoads_,
                                                    rhs,
                                                    *hessianMatrix_,
                                                    i==0    // assemble occupation pattern only for the first call
@@ -323,7 +326,7 @@ void TrustRegionSolver<GridType,VectorType>::solve()
         for (size_t j=0; j<newIterate.size(); j++)
             newIterate[j] += corr[j];
 
-        double energy    = assembler_->computeEnergy(newIterate);
+        double energy    = assembler_->computeEnergy(newIterate, pointLoads_);
 
         // compute the model decrease
         // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
diff --git a/dune/gfe/trustregionsolver.hh b/dune/gfe/trustregionsolver.hh
index 101d1ddfa62c30a1921e2144e94bba54ea308552..064052a7b933f591648407b499326cd5dfd3c438 100644
--- a/dune/gfe/trustregionsolver.hh
+++ b/dune/gfe/trustregionsolver.hh
@@ -66,7 +66,8 @@ public:
                int nu1,
                int nu2,
                int baseIterations,
-               double baseTolerance);
+               double baseTolerance,
+               const SolutionType& pointLoads);
 
     void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
     {
@@ -128,6 +129,9 @@ protected:
 
     /** \brief An L2-norm, really.  The H1SemiNorm class is badly named */
     std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
+
+    SolutionType pointLoads_;
+
 public:
     VectorType identity_;
 };