Solve velocity problem w/ corresponding Dirichlet

src/one-body-sample.py

@@ -13,7 +13,7 @@ class neumannCondition:
 
 class dirichletCondition:
     def __call__(self, x):
-        return x * .005
+        return .005
 
 Functions = {
     'neumannCondition' : neumannCondition(),

src/timestepping.cc

@@ -56,23 +56,20 @@ class ImplicitEuler
     problem_A = this->A;
     problem_A *= this->tau;
 
-    // Treat as Delta u_n for now
-    problem_iterate = this->u_old;
-    if (u_old_old_ptr)
+    // Use the old velocity as an initial iterate if possible
+    if (u_old_old_ptr) {
+      problem_iterate = this->u_old;
       problem_iterate -= *u_old_old_ptr;
+      problem_iterate /= this->tau;
+    } else
+      problem_iterate = 0.0;
 
     for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
       if (this->dirichletNodes[i].count() == dim) {
-        double val;
-        this->dirichletFunction.evaluate(this->time, val);
-        problem_iterate[i] = 0; // Everything prescribed
-        problem_iterate[i][0] =
-            val - this->u_old[i][0]; // Time-dependent X direction
+        problem_iterate[i] = 0;
+        this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
       } else if (this->dirichletNodes[i][1])
-        problem_iterate[i][1] = 0; // Y direction described
-
-    // Make it a velocity
-    problem_iterate /= this->tau;
+        problem_iterate[i][1] = 0; // Y direction prescribed
   }
 
   void virtual extractSolution(VectorType const &problem_iterate,
@@ -112,25 +109,17 @@ class ImplicitTwoStep
     problem_A = this->A;
     problem_A *= 2.0 / 3.0 * this->tau;
 
-    // The finite difference makes a good start (treat it as such for now)
+    // Use the old velocity as an initial iterate
     problem_iterate = this->u_old;
     problem_iterate -= *u_old_old_ptr;
+    problem_iterate /= this->tau;
 
     for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
       if (this->dirichletNodes[i].count() == dim) {
-        double val;
-        this->dirichletFunction.evaluate(this->time, val);
         problem_iterate[i] = 0;
-        problem_iterate[i].axpy(-2, this->u_old[i]);
-        problem_iterate[i].axpy(.5, (*u_old_old_ptr)[i]);
-        problem_iterate[i][0] =
-            1.5 * val - 2 * this->u_old[i][0] + .5 * (*u_old_old_ptr)[i][0];
+        this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
       } else if (this->dirichletNodes[i][1])
-        // Y direction described
-        problem_iterate[i][1] =
-            -2 * this->u_old[i][1] + .5 * (*u_old_old_ptr)[i][1];
-    // Now treat it as a velocity
-    problem_iterate /= this->tau;
+        problem_iterate[i][1] = 0; // Y direction prescribed
   }
 
   void virtual extractSolution(VectorType const &problem_iterate,