diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index caac4622c6b4a268f40e341e5af7d826f545ba93..66ec69f47802daae98f53b6b1845d57d169950ca 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -101,118 +101,149 @@ void setup_boundary(GridView const &gridView,
}
}
+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,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction, double _time)
+ : ell(_ell),
+ A(_A),
+ u_old(_u_old),
+ u_old_old(_u_old_old),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction),
+ time(_time) {}
+
+ virtual ~TimeSteppingScheme() {}
+
+ void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
+ MatrixType &problem_A) const = 0;
+
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const = 0;
+
+ void virtual extractDifference(VectorType const &problem_iterate,
+ VectorType &velocity) const = 0;
+
+protected:
+ VectorType const ℓ
+ MatrixType const &A;
+ VectorType const &u_old;
+ VectorType const *u_old_old;
+ Dune::BitSetVector<dim> const &dirichletNodes;
+ FunctionType const &dirichletFunction;
+ double time;
+};
+
// Implicit Euler: Solve the problem
//
// a(Delta u_new, v - Delta u_new) + j_tau(v) - j_tau(Delta u_new)
// >= l(w - Delta u_new) - a(u_new, v - Delta u_new)
-//
-// Setup: Substract a(u_new, .) from rhs
template <class VectorType, class MatrixType, class FunctionType, int dim>
-void implicitEulerSetup(VectorType const &ell, MatrixType const &A,
- VectorType const &u_old, VectorType const *u_old_old,
- VectorType &problem_rhs, VectorType &problem_iterate,
- MatrixType &problem_A,
- Dune::BitSetVector<dim> const &dirichletNodes,
- FunctionType const &dirichletFunction, double time) {
- problem_A = A;
- problem_rhs = ell;
- problem_A.mmv(u_old, problem_rhs);
-
- problem_iterate = u_old;
- if (u_old_old)
- problem_iterate -= *u_old_old;
-
- for (size_t i = 0; i < dirichletNodes.size(); ++i)
- if (dirichletNodes[i].count() == dim) {
- double val;
- dirichletFunction.evaluate(time, val);
- problem_iterate[i] = 0; // Everything prescribed
- problem_iterate[i][0] = val - u_old[i][0]; // Time-dependent X direction
- } else if (dirichletNodes[i][1])
- problem_iterate[i][1] = 0; // Y direction described
-}
+class ImplicitEuler
+ : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+ // Work arond the fact that nobody implements constructor inheritance
+ template <class... Args>
+ ImplicitEuler(Args &&... args)
+ : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(args...) {
+ }
-// Extraction: Add previous solution
-template <class VectorType>
-void implicitEulerExtract(VectorType const &u_old, VectorType const *u_old_old,
- VectorType const &problem_iterate,
- VectorType &solution) {
- solution = u_old;
- solution += problem_iterate;
-}
+ void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
+ MatrixType &problem_A) const {
+ problem_A = this->A;
+ problem_rhs = this->ell;
+ problem_A.mmv(this->u_old, problem_rhs);
+
+ problem_iterate = this->u_old;
+ if (this->u_old_old)
+ problem_iterate -= *this->u_old_old;
+
+ 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
+ } else if (this->dirichletNodes[i][1])
+ problem_iterate[i][1] = 0; // Y direction described
+ }
-template <class VectorType>
-void implicitEulerExtractVelocity(VectorType const &u_old,
- VectorType const *u_old_old,
- VectorType const &problem_iterate,
- VectorType &diff) {
- diff = problem_iterate;
-}
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const {
+ solution = this->u_old;
+ solution += problem_iterate;
+ }
-// two-Stage implicit algorithm: Solve the problem
-//
-// a(Delta u_new, v - Delta u_new) + j_tau(v) - j_tau(Delta u_new)
-// >= l(w - Delta u_new) - a(u_new, v - Delta u_new)
-//
-// Setup: Substract a(u_new, .) from rhs
+ void virtual extractDifference(VectorType const &problem_iterate,
+ VectorType &velocity) const {
+ velocity = problem_iterate;
+ }
+};
+
+// two-Stage implicit algorithm
template <class VectorType, class MatrixType, class FunctionType, int dim>
-void ImplicitTwoStepSetup(VectorType const &ell, MatrixType const &A,
- VectorType const &u_old, VectorType const *u_old_old,
- VectorType &problem_rhs, VectorType &problem_iterate,
- MatrixType &problem_A,
- Dune::BitSetVector<dim> const &dirichletNodes,
- FunctionType const &dirichletFunction, double time) {
- problem_A = A;
- problem_A /= 1.5;
- problem_rhs = ell;
- problem_A.usmv(-2, u_old, problem_rhs);
- problem_A.usmv(.5, *u_old_old, problem_rhs);
-
- // The finite difference makes a good start
- problem_iterate = u_old;
- problem_iterate -= *u_old_old;
-
- for (size_t i = 0; i < dirichletNodes.size(); ++i)
- if (dirichletNodes[i].count() == dim) {
- double val;
- dirichletFunction.evaluate(time, val);
- problem_iterate[i] = 0;
- problem_iterate[i].axpy(-2, u_old[i]);
- problem_iterate[i].axpy(.5, (*u_old_old)[i]);
- problem_iterate[i][0] =
- 1.5 * val - 2 * u_old[i][0] + .5 * (*u_old_old)[i][0];
- } else if (dirichletNodes[i][1])
- // Y direction described
- problem_iterate[i][1] = -2 * u_old[i][1] + .5 * (*u_old_old)[i][1];
-}
+class ImplicitTwoStep
+ : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+ // Work arond the fact that nobody implements constructor inheritance
+ template <class... Args>
+ ImplicitTwoStep(Args &&... args)
+ : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(args...) {
+ }
-template <class VectorType>
-void ImplicitTwoStepExtract(VectorType const &u_old,
- VectorType const *u_old_old,
- VectorType const &problem_iterate,
- VectorType &solution) {
- solution = problem_iterate;
- solution.axpy(2, u_old);
- solution.axpy(-.5, *u_old_old);
- solution *= 2.0 / 3.0;
-
- // Check if we split correctly
- {
- VectorType test = problem_iterate;
- test.axpy(-1.5, solution);
- test.axpy(+2, u_old);
- test.axpy(-.5, *u_old_old);
- assert(test.two_norm() < 1e-10);
+ void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
+ MatrixType &problem_A) const {
+ problem_A = this->A;
+ problem_A /= 1.5;
+ problem_rhs = this->ell;
+ problem_A.usmv(-2, this->u_old, problem_rhs);
+ problem_A.usmv(.5, *this->u_old_old, problem_rhs);
+
+ // The finite difference makes a good start
+ problem_iterate = this->u_old;
+ problem_iterate -= *this->u_old_old;
+
+ 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, (*this->u_old_old)[i]);
+ problem_iterate[i][0] =
+ 1.5 * val - 2 * this->u_old[i][0] + .5 * (*this->u_old_old)[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)[i][1];
}
-}
-template <class VectorType>
-void ImplicitTwoStepExtractVelocity(VectorType const &u_old,
- VectorType const *u_old_old,
- VectorType const &problem_iterate,
- VectorType &diff) {
- diff = problem_iterate;
-}
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const {
+ solution = problem_iterate;
+ solution.axpy(2, this->u_old);
+ solution.axpy(-.5, *this->u_old_old);
+ solution *= 2.0 / 3.0;
+
+ // Check if we split correctly
+ {
+ VectorType test = problem_iterate;
+ test.axpy(-1.5, solution);
+ test.axpy(+2, this->u_old);
+ test.axpy(-.5, *this->u_old_old);
+ assert(test.two_norm() < 1e-10);
+ }
+ }
+
+ void virtual extractDifference(VectorType const &problem_iterate,
+ VectorType &velocity) const {
+ velocity = problem_iterate;
+ }
+};
int main(int argc, char *argv[]) {
try {
@@ -368,17 +399,26 @@ int main(int argc, char *argv[]) {
VectorType problem_rhs(finestSize);
VectorType problem_iterate(finestSize);
- auto setupFunc =
- (run == 1 || !parset.get<bool>("ImplicitTwoStep"))
- ? &implicitEulerSetup<VectorType, MatrixType,
- decltype(dirichletFunction), dim>
- : &ImplicitTwoStepSetup<VectorType, MatrixType,
- decltype(dirichletFunction), dim>;
+ VectorType *u_old_old = (run == 1) ? nullptr : &u_previous;
+
+ typedef TimeSteppingScheme<VectorType, MatrixType,
+ decltype(dirichletFunction), dim> TS;
+ typedef ImplicitEuler<VectorType, MatrixType,
+ decltype(dirichletFunction), dim> IE;
+ typedef ImplicitTwoStep<VectorType, MatrixType,
+ decltype(dirichletFunction), dim> ITS;
+
+ Dune::shared_ptr<TS> timeSteppingScheme;
+ if (run == 1 || !parset.get<bool>("implicitTwoStep"))
+ timeSteppingScheme =
+ Dune::make_shared<IE>(ell, stiffnessMatrix, u, u_old_old,
+ ignoreNodes, dirichletFunction, time);
+ else
+ timeSteppingScheme =
+ Dune::make_shared<ITS>(ell, stiffnessMatrix, u, u_old_old,
+ ignoreNodes, dirichletFunction, time);
- VectorType *u_old_old_ptr = (run == 1) ? nullptr : &u_previous;
- setupFunc(ell, stiffnessMatrix, u, u_old_old_ptr, problem_rhs,
- problem_iterate, problem_A, ignoreNodes, dirichletFunction,
- time);
+ timeSteppingScheme->setup(problem_rhs, problem_iterate, problem_A);
VectorType solution_saved = u;
auto const state_fpi_max =
@@ -400,20 +440,9 @@ int main(int argc, char *argv[]) {
false); // absolute error
overallSolver.solve();
- auto extractFunc = (run == 1 || !parset.get<bool>("ImplicitTwoStep"))
- ? implicitEulerExtract<VectorType>
- : ImplicitTwoStepExtract<VectorType>;
-
- // Extract solution from solver
- extractFunc(u, u_old_old_ptr, problem_iterate, solution);
-
- auto extractDiffFunc =
- (run == 1 || !parset.get<bool>("implicitTwoStep"))
- ? implicitEulerExtractVelocity<VectorType>
- : ImplicitTwoStepExtractVelocity<VectorType>;
-
- // Extract difference from solver
- extractDiffFunc(u, u_old_old_ptr, problem_iterate, u_diff);
+ timeSteppingScheme->extractSolution(problem_iterate, solution);
+ timeSteppingScheme->extractDifference(problem_iterate, u_diff);
+ // FIXME: memory leak
// Update the state
for (size_t i = 0; i < frictionalNodes.size(); ++i) {