Skip to content
Snippets Groups Projects
Commit 445acbcb authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

u_old_old -> u_old_old_ptr

parent 35bd0f6d
No related branches found
No related tags found
No related merge requests found
...@@ -2,14 +2,14 @@ template <class VectorType, class MatrixType, class FunctionType, int dim> ...@@ -2,14 +2,14 @@ template <class VectorType, class MatrixType, class FunctionType, int dim>
class TimeSteppingScheme { class TimeSteppingScheme {
public: public:
TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A, TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A,
VectorType const &_u_old, VectorType const *_u_old_old, VectorType const &_u_old, VectorType const *_u_old_old_ptr,
Dune::BitSetVector<dim> const &_dirichletNodes, Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction, double _time, FunctionType const &_dirichletFunction, double _time,
double _tau) double _tau)
: ell(_ell), : ell(_ell),
A(_A), A(_A),
u_old(_u_old), u_old(_u_old),
u_old_old(_u_old_old), u_old_old_ptr(_u_old_old_ptr),
dirichletNodes(_dirichletNodes), dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction), dirichletFunction(_dirichletFunction),
time(_time), time(_time),
...@@ -30,7 +30,7 @@ class TimeSteppingScheme { ...@@ -30,7 +30,7 @@ class TimeSteppingScheme {
VectorType const &ell; VectorType const &ell;
MatrixType const &A; MatrixType const &A;
VectorType const &u_old; VectorType const &u_old;
VectorType const *u_old_old; VectorType const *u_old_old_ptr;
Dune::BitSetVector<dim> const &dirichletNodes; Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction; FunctionType const &dirichletFunction;
double const time; double const time;
...@@ -61,8 +61,8 @@ class ImplicitEuler ...@@ -61,8 +61,8 @@ class ImplicitEuler
// Treat as Delta u_n for now // Treat as Delta u_n for now
problem_iterate = this->u_old; problem_iterate = this->u_old;
if (this->u_old_old) if (this->u_old_old_ptr)
problem_iterate -= *this->u_old_old; problem_iterate -= *this->u_old_old_ptr;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i) for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) { if (this->dirichletNodes[i].count() == dim) {
...@@ -105,14 +105,14 @@ class ImplicitTwoStep ...@@ -105,14 +105,14 @@ class ImplicitTwoStep
MatrixType &problem_A) const { MatrixType &problem_A) const {
problem_rhs = this->ell; problem_rhs = this->ell;
this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs); this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs);
this->A.usmv(+1.0 / 3.0, *this->u_old_old, problem_rhs); this->A.usmv(+1.0 / 3.0, *this->u_old_old_ptr, problem_rhs);
problem_A = this->A; problem_A = this->A;
problem_A *= 2.0 / 3.0 * this->tau; problem_A *= 2.0 / 3.0 * this->tau;
// The finite difference makes a good start (treat it as such for now) // The finite difference makes a good start (treat it as such for now)
problem_iterate = this->u_old; problem_iterate = this->u_old;
problem_iterate -= *this->u_old_old; problem_iterate -= *this->u_old_old_ptr;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i) for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) { if (this->dirichletNodes[i].count() == dim) {
...@@ -120,13 +120,13 @@ class ImplicitTwoStep ...@@ -120,13 +120,13 @@ class ImplicitTwoStep
this->dirichletFunction.evaluate(this->time, val); this->dirichletFunction.evaluate(this->time, val);
problem_iterate[i] = 0; problem_iterate[i] = 0;
problem_iterate[i].axpy(-2, this->u_old[i]); problem_iterate[i].axpy(-2, this->u_old[i]);
problem_iterate[i].axpy(.5, (*this->u_old_old)[i]); problem_iterate[i].axpy(.5, (*this->u_old_old_ptr)[i]);
problem_iterate[i][0] = problem_iterate[i][0] = 1.5 * val - 2 * this->u_old[i][0] +
1.5 * val - 2 * this->u_old[i][0] + .5 * (*this->u_old_old)[i][0]; .5 * (*this->u_old_old_ptr)[i][0];
} else if (this->dirichletNodes[i][1]) } else if (this->dirichletNodes[i][1])
// Y direction described // Y direction described
problem_iterate[i][1] = problem_iterate[i][1] =
-2 * this->u_old[i][1] + .5 * (*this->u_old_old)[i][1]; -2 * this->u_old[i][1] + .5 * (*this->u_old_old_ptr)[i][1];
// Now treat it as a velocity // Now treat it as a velocity
problem_iterate /= this->tau; problem_iterate /= this->tau;
} }
...@@ -136,7 +136,7 @@ class ImplicitTwoStep ...@@ -136,7 +136,7 @@ class ImplicitTwoStep
solution = problem_iterate; solution = problem_iterate;
solution *= this->tau; solution *= this->tau;
solution.axpy(2, this->u_old); solution.axpy(2, this->u_old);
solution.axpy(-.5, *this->u_old_old); solution.axpy(-.5, *this->u_old_old_ptr);
solution *= 2.0 / 3.0; solution *= 2.0 / 3.0;
// Check if we split correctly // Check if we split correctly
...@@ -145,7 +145,7 @@ class ImplicitTwoStep ...@@ -145,7 +145,7 @@ class ImplicitTwoStep
test *= this->tau; test *= this->tau;
test.axpy(-1.5, solution); test.axpy(-1.5, solution);
test.axpy(+2, this->u_old); test.axpy(+2, this->u_old);
test.axpy(-.5, *this->u_old_old); test.axpy(-.5, *this->u_old_old_ptr);
assert(test.two_norm() < 1e-10); assert(test.two_norm() < 1e-10);
} }
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment