From bcf73446fdcd39a9205a371ad4cc25ce0d3c35b4 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 10 Feb 2013 19:49:23 +0100
Subject: [PATCH] Implement eulerPair scheme (implicitEuler for u,v)
---
src/enum_scheme.cc | 3 ++
src/enums.hh | 3 +-
src/one-body-sample.cc | 5 ++
src/timestepping.cc | 106 +++++++++++++++++++++++++++++++++++++++
src/timestepping.hh | 37 ++++++++++++++
src/timestepping_tmpl.cc | 1 +
6 files changed, 154 insertions(+), 1 deletion(-)
diff --git a/src/enum_scheme.cc b/src/enum_scheme.cc
index 9866d4d5..615d8933 100644
--- a/src/enum_scheme.cc
+++ b/src/enum_scheme.cc
@@ -11,6 +11,9 @@ template <> struct StringToEnum<Config::scheme> {
if (s == "newmark")
return Config::Newmark;
+ if (s == "eulerPair")
+ return Config::EulerPair;
+
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
};
diff --git a/src/enums.hh b/src/enums.hh
index f8d05d21..b6ffb733 100644
--- a/src/enums.hh
+++ b/src/enums.hh
@@ -9,7 +9,8 @@ struct Config {
enum scheme {
ImplicitTwoStep,
ImplicitEuler,
- Newmark
+ Newmark,
+ EulerPair
};
};
#endif
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index c978a555..4525e6bb 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -108,6 +108,11 @@ initTimeStepper(Config::scheme scheme, FunctionType const &dirichletFunction,
Newmark<VectorType, MatrixType, FunctionType, dims>>(
stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
ignoreNodes, dirichletFunction);
+ case Config::EulerPair:
+ return Dune::make_shared<
+ EulerPair<VectorType, MatrixType, FunctionType, dims>>(
+ stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
+ ignoreNodes, dirichletFunction);
}
}
template <class SingletonVectorType, class VectorType>
diff --git a/src/timestepping.cc b/src/timestepping.cc
index 07d1bf41..614ed718 100644
--- a/src/timestepping.cc
+++ b/src/timestepping.cc
@@ -337,5 +337,111 @@ TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
Newmark<VectorType, MatrixType, FunctionType, dim>::clone() {
return new Newmark<VectorType, MatrixType, FunctionType, dim>(*this);
}
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+EulerPair<VectorType, MatrixType, FunctionType, dim>::EulerPair(
+ MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
+ VectorType const &_ud_initial, VectorType const &_udd_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction)
+ : A(_A),
+ B(_B),
+ u(_u_initial),
+ ud(_ud_initial),
+ udd(_udd_initial),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+ udd_old = udd;
+ ud_old = ud;
+ u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::setup(
+ VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+ VectorType &problem_iterate, MatrixType &problem_A) {
+ postProcessCalled = false;
+
+ tau = _tau;
+
+ problem_rhs = ell;
+ B.usmv(1.0 / tau, ud_old, problem_rhs);
+ A.usmv(-1.0, u_old, problem_rhs);
+
+ // For fixed tau, we'd only really have to do this once
+ Dune::MatrixIndexSet indices(A.N(), A.M());
+ indices.import(A);
+ indices.import(B);
+ indices.exportIdx(problem_A);
+ problem_A = 0.0;
+ Arithmetic::addProduct(problem_A, tau, A);
+ Arithmetic::addProduct(problem_A, 1.0 / tau, B);
+
+ // ud_old makes a good initial iterate; we could use anything, though
+ problem_iterate = ud_old;
+
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ switch (dirichletNodes[i].count()) {
+ case 0:
+ continue;
+ case dim:
+ problem_iterate[i] = 0;
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ case 1:
+ if (dirichletNodes[i][0]) {
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ }
+ if (dirichletNodes[i][1]) {
+ problem_iterate[i][1] = 0;
+ break;
+ }
+ assert(false);
+ default:
+ assert(false);
+ }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::postProcess(
+ VectorType const &problem_iterate) {
+ postProcessCalled = true;
+
+ ud = problem_iterate;
+
+ u = u_old;
+ Arithmetic::addProduct(u, tau, ud);
+
+ udd = 0;
+ Arithmetic::addProduct(udd, 1.0 / tau, ud);
+ Arithmetic::addProduct(udd, -1.0 / tau, ud_old);
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractDisplacement(
+ VectorType &displacement) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
+ VectorType &velocity) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ velocity = ud;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+EulerPair<VectorType, MatrixType, FunctionType, dim>::clone() {
+ return new EulerPair<VectorType, MatrixType, FunctionType, dim>(*this);
+}
#include "timestepping_tmpl.cc"
diff --git a/src/timestepping.hh b/src/timestepping.hh
index 84da56b7..f57d67ff 100644
--- a/src/timestepping.hh
+++ b/src/timestepping.hh
@@ -112,6 +112,43 @@ class Newmark
virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
clone();
+private:
+ MatrixType const &A;
+ MatrixType const &B;
+ VectorType u;
+ VectorType ud;
+ VectorType udd;
+ Dune::BitSetVector<dim> const &dirichletNodes;
+ FunctionType const &dirichletFunction;
+
+ VectorType u_old;
+ VectorType ud_old;
+ VectorType udd_old;
+
+ double tau;
+
+ bool postProcessCalled = false;
+};
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class EulerPair
+ : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+ EulerPair(MatrixType const &_A, MatrixType const &_B,
+ VectorType const &_u_initial, VectorType const &_ud_initial,
+ VectorType const &_udd_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction);
+
+ void virtual nextTimeStep();
+ void virtual setup(VectorType const &, double, double, VectorType &,
+ VectorType &, MatrixType &);
+ void virtual postProcess(VectorType const &);
+ void virtual extractDisplacement(VectorType &) const;
+ void virtual extractVelocity(VectorType &) const;
+
+ virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+ clone();
+
private:
MatrixType const &A;
MatrixType const &B;
diff --git a/src/timestepping_tmpl.cc b/src/timestepping_tmpl.cc
index 3db2ba60..5feceaf6 100644
--- a/src/timestepping_tmpl.cc
+++ b/src/timestepping_tmpl.cc
@@ -17,3 +17,4 @@ typedef Dune::VirtualFunction<double, double> FunctionType;
template class ImplicitEuler<VectorType, MatrixType, FunctionType, DIM>;
template class ImplicitTwoStep<VectorType, MatrixType, FunctionType, DIM>;
template class Newmark<VectorType, MatrixType, FunctionType, DIM>;
+template class EulerPair<VectorType, MatrixType, FunctionType, DIM>;
--
GitLab