Force zero normal displacement on fric. boundary

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

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

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;