Move u_old_old_ptr from TimeSteppingScheme

src/timestepping.cc

@@ -2,14 +2,13 @@ template <class VectorType, class MatrixType, class FunctionType, int dim>
 class TimeSteppingScheme {
 public:
   TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A,
-                     VectorType const &_u_old, VectorType const *_u_old_old_ptr,
+                     VectorType const &_u_old,
                      Dune::BitSetVector<dim> const &_dirichletNodes,
                      FunctionType const &_dirichletFunction, double _time,
                      double _tau)
       : ell(_ell),
         A(_A),
         u_old(_u_old),
-        u_old_old_ptr(_u_old_old_ptr),
         dirichletNodes(_dirichletNodes),
         dirichletFunction(_dirichletFunction),
         time(_time),
@@ -30,7 +29,6 @@ class TimeSteppingScheme {
   VectorType const &ell;
   MatrixType const &A;
   VectorType const &u_old;
-  VectorType const *u_old_old_ptr;
   Dune::BitSetVector<dim> const &dirichletNodes;
   FunctionType const &dirichletFunction;
   double const time;
@@ -45,14 +43,14 @@ template <class VectorType, class MatrixType, class FunctionType, int dim>
 class ImplicitEuler
     : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
 public:
-  ImplicitEuler(VectorType const &_ell, MatrixType const &_A,
-                VectorType const &_u_old, VectorType const *_u_old_old_ptr,
-                Dune::BitSetVector<dim> const &_dirichletNodes,
-                FunctionType const &_dirichletFunction, double _time,
-                double _tau)
+  ImplicitEuler(
+      VectorType const &_ell, MatrixType const &_A, VectorType const &_u_old,
+      VectorType const *_u_old_old_ptr, // FIXME: this should not be necessary
+      Dune::BitSetVector<dim> const &_dirichletNodes,
+      FunctionType const &_dirichletFunction, double _time, double _tau)
       : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
-            _ell, _A, _u_old, _u_old_old_ptr, _dirichletNodes,
-            _dirichletFunction, _time, _tau) {}
+            _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
+        u_old_old_ptr(_u_old_old_ptr) {}
 
   void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
                      MatrixType &problem_A) const {
@@ -64,8 +62,8 @@ class ImplicitEuler
 
     // Treat as Delta u_n for now
     problem_iterate = this->u_old;
-    if (this->u_old_old_ptr)
-      problem_iterate -= *this->u_old_old_ptr;
+    if (u_old_old_ptr)
+      problem_iterate -= *u_old_old_ptr;
 
     for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
       if (this->dirichletNodes[i].count() == dim) {
@@ -91,6 +89,9 @@ class ImplicitEuler
                                VectorType &velocity) const {
     velocity = problem_iterate;
   }
+
+private:
+  VectorType const *u_old_old_ptr;
 };
 
 // two-Stage implicit algorithm
@@ -104,21 +105,21 @@ class ImplicitTwoStep
                   FunctionType const &_dirichletFunction, double _time,
                   double _tau)
       : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
-            _ell, _A, _u_old, _u_old_old_ptr, _dirichletNodes,
-            _dirichletFunction, _time, _tau) {}
+            _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
+        u_old_old_ptr(_u_old_old_ptr) {}
 
   void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
                      MatrixType &problem_A) const {
     problem_rhs = this->ell;
     this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs);
-    this->A.usmv(+1.0 / 3.0, *this->u_old_old_ptr, problem_rhs);
+    this->A.usmv(+1.0 / 3.0, *u_old_old_ptr, problem_rhs);
 
     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)
     problem_iterate = this->u_old;
-    problem_iterate -= *this->u_old_old_ptr;
+    problem_iterate -= *u_old_old_ptr;
 
     for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
       if (this->dirichletNodes[i].count() == dim) {
@@ -126,13 +127,13 @@ class ImplicitTwoStep
         this->dirichletFunction.evaluate(this->time, val);
         problem_iterate[i] = 0;
         problem_iterate[i].axpy(-2, this->u_old[i]);
-        problem_iterate[i].axpy(.5, (*this->u_old_old_ptr)[i]);
-        problem_iterate[i][0] = 1.5 * val - 2 * this->u_old[i][0] +
-                                .5 * (*this->u_old_old_ptr)[i][0];
+        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];
       } else if (this->dirichletNodes[i][1])
         // Y direction described
         problem_iterate[i][1] =
-            -2 * this->u_old[i][1] + .5 * (*this->u_old_old_ptr)[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;
   }
@@ -142,7 +143,7 @@ class ImplicitTwoStep
     solution = problem_iterate;
     solution *= this->tau;
     solution.axpy(2, this->u_old);
-    solution.axpy(-.5, *this->u_old_old_ptr);
+    solution.axpy(-.5, *u_old_old_ptr);
     solution *= 2.0 / 3.0;
 
     // Check if we split correctly
@@ -151,7 +152,7 @@ class ImplicitTwoStep
       test *= this->tau;
       test.axpy(-1.5, solution);
       test.axpy(+2, this->u_old);
-      test.axpy(-.5, *this->u_old_old_ptr);
+      test.axpy(-.5, *u_old_old_ptr);
       assert(test.two_norm() < 1e-10);
     }
   }
@@ -160,4 +161,7 @@ class ImplicitTwoStep
                                VectorType &velocity) const {
     velocity = problem_iterate;
   }
+
+private:
+  VectorType const *u_old_old_ptr;
 };