From 2eb3a6b0db81c801ff620e54de797a10ed32f553 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 27 Oct 2014 20:22:09 +0000
Subject: [PATCH] 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 | 9 ++++++---
 dune/gfe/trustregionsolver.hh | 6 +++++-
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/dune/gfe/trustregionsolver.cc b/dune/gfe/trustregionsolver.cc
index 7bd14731..df3dfd33 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 101d1ddf..064052a7 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_;
 };
-- 
GitLab