Add (temporary) support for point loads

To really make sure I am doing the infamous L-shape example
like the people before me did I need point loads.  Being a
mathematician I am no particular fan of point loads, so they
may go out again eventually.

[[Imported from SVN: r9942]]

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>

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