From 11ae5151979913d4ebb9c0b4928662d1f1209661 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 7 Nov 2011 14:01:24 +0100
Subject: [PATCH] Force zero normal displacement on fric. boundary

---
 dune/tectonic/myblockproblem.hh   | 21 +++++++++++++++++----
 dune/tectonic/samplefunctional.hh | 11 +++++++++--
 src/one-body-sample.cc            |  1 +
 3 files changed, 27 insertions(+), 6 deletions(-)

diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index a3cc1706..b8a50d9d 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -84,9 +84,21 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject {
       LocalVectorType &ui, int m,
       const typename Dune::BitSetVector<block_size>::const_reference ignore) {
     {
-      // TODO: Does it make any sense to ignore single spatial dimensions here?
-      if (ignore.test(0))
-        return;
+      int ignore_component =
+          block_size; // Special value that indicates nothing should be ignored
+      switch (ignore.count()) {
+        case 0: // Full problem
+          break;
+        case 1: // 1 Dimension is fixed
+          assert(
+              ignore[1]); // Only the second component is allowed to stay fixed
+          ignore_component = 1;
+          break;
+        case block_size: // Ignore the whole node
+          return;
+        default:
+          assert(false);
+      }
 
       LocalMatrixType const *localA = NULL;
       LocalVectorType localb(problem.f[m]);
@@ -105,7 +117,8 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject {
       FunctionType f;
       problem.newphi.restriction(m, f);
       Dune::LocalNonlinearity<block_size> const phi(f);
-      Dune::SampleFunctional<block_size> localJ(*localA, localb, phi);
+      Dune::SampleFunctional<block_size> localJ(*localA, localb, phi,
+                                                ignore_component);
 
       LocalVectorType correction;
       Dune::minimise(localJ, ui, 10, bisection); // FIXME: hardcoded value
diff --git a/dune/tectonic/samplefunctional.hh b/dune/tectonic/samplefunctional.hh
index 1d510f54..a19d40a5 100644
--- a/dune/tectonic/samplefunctional.hh
+++ b/dune/tectonic/samplefunctional.hh
@@ -23,8 +23,8 @@ template <int dim> class SampleFunctional {
   typedef LocalNonlinearity<dim> NonlinearityType;
 
   SampleFunctional(SmallMatrix const &A, SmallVector const &b,
-                   NonlinearityType const &phi)
-      : A(A), b(b), phi(phi) {}
+                   NonlinearityType const &phi, int ignore = dim)
+      : A(A), b(b), phi(phi), ignore(ignore) {}
 
   double operator()(SmallVector const &v) const {
     SmallVector y;
@@ -89,12 +89,15 @@ template <int dim> class SampleFunctional {
   SmallMatrix const &A;
   SmallVector const &b;
   NonlinearityType const &phi;
+  int const ignore;
 
 private:
   // Gradient of the smooth part
   void smoothGradient(SmallVector const &x, SmallVector &y) const {
     A.mv(x, y); // y = Av
     y -= b;     // y = Av - b
+    if (ignore != dim)
+      y[ignore] = 0;
   }
 
   void upperGradient(SmallVector const &x, SmallVector &y) const {
@@ -102,6 +105,8 @@ template <int dim> class SampleFunctional {
     SmallVector z;
     smoothGradient(x, z);
     y += z;
+    if (ignore != dim)
+      y[ignore] = 0;
   }
 
   void lowerGradient(SmallVector const &x, SmallVector &y) const {
@@ -109,6 +114,8 @@ template <int dim> class SampleFunctional {
     SmallVector z;
     smoothGradient(x, z);
     y += z;
+    if (ignore != dim)
+      y[ignore] = 0;
   }
 
   // y = (id - P)(d) where P is the projection onto the line t*x
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index e7e9af94..2ce6e9e3 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -79,6 +79,7 @@ void setup_boundary(GridView const &gridView,
       ++frictional_nodes;
       size_t const id = myVertexMapper.map(*it);
       frictionalNodes[id] = true;
+      ignoreNodes[id][1] = true; // Zero displacement in direction y
     }
   }
   std::cout << "Number of Neumann nodes: " << neumann_nodes << std::endl;
-- 
GitLab