Compare the outcome, not just the last step

Make Neumann boundary t-dependent as well

src/one-body-sample.cc

@@ -89,8 +89,8 @@ void setup_boundary(GridView const &gridView,
 template <class GridType, class GridView, class LocalVectorType, class FEBasis>
 void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
                       Dune::BitSetVector<1> const &neumannNodes,
-                      Dune::BlockVector<LocalVectorType> &
-                          f) { // constant sample function on neumann boundary
+                      Dune::BlockVector<LocalVectorType> &f,
+                      int run) { // constant sample function on neumann boundary
   BoundaryPatch<GridView> neumannBoundary(gridView, neumannNodes);
   LocalVectorType SampleVector(0);
   // FIXME: random values
@@ -188,7 +188,8 @@ int main() {
     u1 = 0;
     VectorType u2 = u1;
 
-    VectorType b(grid.size(grid.maxLevel(), dim));
+    VectorType b1(grid.size(grid.maxLevel(), dim));
+    VectorType b2(grid.size(grid.maxLevel(), dim));
 
     std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
     assemble_frictional<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
@@ -198,10 +199,12 @@ int main() {
       std::cout << "Run: " << run << std::endl;
       // b = neumann
       assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
-          leafView, p1Basis, neumannNodes, b);
+          leafView, p1Basis, neumannNodes, b1, run);
+      b2 = b1;
 
       // b += linear update
-      stiffnessMatrix.umv(u1, b);
+      stiffnessMatrix.umv(u1, b1);
+      stiffnessMatrix.umv(u2, b2);
 
       // {{{ Assemble terms for the nonlinearity
       // TODO: Random value
@@ -221,7 +224,7 @@ int main() {
 
       {
         MyConvexProblemType myConvexProblem(stiffnessMatrix,
-                                            myGlobalNonlinearity, b, u1);
+                                            myGlobalNonlinearity, b1, u1);
         MyBlockProblemType myBlockProblem(myConvexProblem);
         nonlinearGSStep.setProblem(u1, myBlockProblem);
 
@@ -260,7 +263,7 @@ int main() {
       {
         // TODO: Why does blockGSStep even provide a default constructor?
         BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u2,
-                                                          b);
+                                                          b2);
         blockGSStep.ignoreNodes_ = &ignoreNodes;
 
         LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,