Bug fix: We're computation increments here!

Also, time-independent neumann conditions

src/one-body-sample.cc

@@ -96,8 +96,8 @@ void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
                       int run) { // constant sample function on neumann boundary
   BoundaryPatch<GridView> neumannBoundary(gridView, neumannNodes);
   LocalVectorType SampleVector(0);
-  // FIXME: random values
-  SampleVector[0] = 1 + run / 200;
+  // FIXME: random values (time-independent)
+  SampleVector[0] = 1;
   SampleVector[1] = 0;
   ConstantFunction<LocalVectorType, LocalVectorType> fNeumann(SampleVector);
   NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
@@ -194,6 +194,15 @@ int main(int argc, char *argv[]) {
     VectorType u1(grid.size(grid.maxLevel(), dim));
     u1 = 0;
     VectorType u3 = u1;
+
+    VectorType u1_diff_old(grid.size(grid.maxLevel(), dim));
+    u1_diff_old = 0;
+    VectorType u3_diff_old = u1_diff_old;
+
+    VectorType u1_diff_new(grid.size(grid.maxLevel(), dim));
+    u1_diff_new = 0;
+    VectorType u3_diff_new = u1_diff_new;
+
     CellVectorType vonMisesStress;
 
     VectorType b1;
@@ -213,8 +222,8 @@ int main(int argc, char *argv[]) {
       b3 = b1;
 
       // b += linear update
-      stiffnessMatrix.umv(u1, b1);
-      stiffnessMatrix.umv(u3, b3);
+      stiffnessMatrix.umv(u1_diff_old, b1);
+      stiffnessMatrix.umv(u3_diff_old, b3);
 
       // {{{ Assemble terms for the nonlinearity
 
@@ -258,10 +267,10 @@ int main(int argc, char *argv[]) {
       // }}}
 
       {
-        MyConvexProblemType myConvexProblem(stiffnessMatrix,
-                                            *myGlobalNonlinearity, b1, u1);
+        MyConvexProblemType myConvexProblem(
+            stiffnessMatrix, *myGlobalNonlinearity, b1, u1_diff_new);
         MyBlockProblemType myBlockProblem(parset, myConvexProblem);
-        nonlinearGSStep.setProblem(u1, myBlockProblem);
+        nonlinearGSStep.setProblem(u1_diff_new, myBlockProblem);
 
         LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
                                       solver_tolerance, &energyNorm,
@@ -269,6 +278,9 @@ int main(int argc, char *argv[]) {
         solver.solve();
       }
 
+      u1 += u1_diff_new;
+      u1_diff_old = u1_diff_new;
+
       auto *displacement =
           new BasisGridFunction<P1Basis, VectorType>(p1Basis, u1);
       VonMisesStressAssembler<GridType> localStressAssembler(E, nu,
@@ -294,7 +306,7 @@ int main(int argc, char *argv[]) {
       // same results if phi vanishes (e.g. because the normalstress is zero)
       {
         TruncatedBlockGSStep<OperatorType, VectorType> blockGSStep(
-            stiffnessMatrix, u3, b3);
+            stiffnessMatrix, u3_diff_new, b3);
         blockGSStep.ignoreNodes_ = &ignoreNodes;
 
         LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
@@ -302,6 +314,9 @@ int main(int argc, char *argv[]) {
                                       Solver::QUIET);
         solver.solve();
       }
+
+      u3 += u3_diff_new;
+      u3_diff_old = u3_diff_new;
     }
     std::cout << std::endl;