[Cleanup] Bug fix: Do not mix old and new BC

During the first fixed point iteration, velocityIterate and the old
velocity are identical except for the new BC: velocityIterate has
already had them imposed.

To simplify the control flow, kill v completely.

Note that since we always make two fixed point iterations and this issue
only occurs during the first step, this change only affects
computational results as (a reduction of) noise

src/sand-wedge.cc

@@ -352,7 +352,6 @@ int main(int argc, char *argv[]) {
         parset.get<double>("boundary.friction.L"),
         parset.get<double>("boundary.friction.V0"));
 
-    Vector v = v_initial;
     Vector v_m(leafVertexCount);
     ScalarVector alpha(leafVertexCount);
 
@@ -413,7 +412,7 @@ int main(int argc, char *argv[]) {
       for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
         timeSteppingScheme->extractOldVelocity(v_m);
         v_m *= 1.0 - lambda;
-        Arithmetic::addProduct(v_m, lambda, v);
+        Arithmetic::addProduct(v_m, lambda, velocityIterate);
 
         stateUpdater->solve(v_m);
         stateUpdater->extractAlpha(alpha);
@@ -421,21 +420,21 @@ int main(int argc, char *argv[]) {
         solveVelocityProblem(velocityIterate, alpha);
         timeSteppingScheme->postProcess(velocityIterate);
         timeSteppingScheme->extractDisplacement(u);
-        timeSteppingScheme->extractVelocity(v);
 
         iterationWriter << iterationCounter << " ";
         if (printProgress)
           std::cout << '.' << std::flush;
 
         if (stateFPI > 1) {
-          double const velocityCorrection = AMNorm.diff(v_saved, v);
+          double const velocityCorrection =
+              AMNorm.diff(v_saved, velocityIterate);
           if (velocityCorrection < fixedPointTolerance)
             break;
         }
         if (stateFPI == maximumStateFPI)
           DUNE_THROW(Dune::Exception, "FPI failed to converge");
 
-        v_saved = v;
+        v_saved = velocityIterate;
       }
       if (printProgress)
         std::cout << std::endl;