Bug fix: Use u instead of delta_u in linear update

src/one-body-sample.cc

@@ -206,11 +206,6 @@ int main(int argc, char *argv[]) {
     VectorType u2 = u1;
     VectorType u3 = u1;
 
-    VectorType u1_diff_old(grid.size(grid.maxLevel(), dim));
-    u1_diff_old = 0;
-    VectorType u2_diff_old = u1_diff_old;
-    VectorType u3_diff_old = u1_diff_old;
-
     VectorType u1_diff_new(grid.size(grid.maxLevel(), dim));
     u1_diff_new = 0;
     VectorType u2_diff_new = u1_diff_new;
@@ -237,9 +232,9 @@ int main(int argc, char *argv[]) {
       b3 = b1;
 
       // b -= linear update
-      stiffnessMatrix.mmv(u1_diff_old, b1);
-      stiffnessMatrix.mmv(u2_diff_old, b2);
-      stiffnessMatrix.mmv(u3_diff_old, b3);
+      stiffnessMatrix.mmv(u1, b1);
+      stiffnessMatrix.mmv(u2, b2);
+      stiffnessMatrix.mmv(u3, b3);
 
       // {{{ Assemble terms for the nonlinearity
 
@@ -295,7 +290,6 @@ int main(int argc, char *argv[]) {
       }
 
       u1 += u1_diff_new;
-      u1_diff_old = u1_diff_new;
 
       auto *displacement =
           new BasisGridFunction<P1Basis, VectorType>(p1Basis, u1);
@@ -329,7 +323,6 @@ int main(int argc, char *argv[]) {
       }
 
       u2 += u2_diff_new;
-      u2_diff_old = u2_diff_new;
 
       // Use a linear solver for comparison; should return roughly the
       // same results if phi vanishes (e.g. because the normalstress is zero)
@@ -344,7 +337,6 @@ int main(int argc, char *argv[]) {
       }
 
       u3 += u3_diff_new;
-      u3_diff_old = u3_diff_new;
     }
     std::cout << std::endl;