From c6387cb64a93eacbbf22dc9853ad891c2a6c5053 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 11 Sep 2012 14:30:52 +0200
Subject: [PATCH] Move u_old_old_ptr from TimeSteppingScheme
---
src/timestepping.cc | 48 ++++++++++++++++++++++++---------------------
1 file changed, 26 insertions(+), 22 deletions(-)
diff --git a/src/timestepping.cc b/src/timestepping.cc
index dee746d0..caf7f1c0 100644
--- a/src/timestepping.cc
+++ b/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 ℓ
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;
};
--
GitLab