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

[Cleanup] Move time stepping logic into parent class 

While we're at it, reset postProcessCalled in nextTimeStep()
parent 21a6b769
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/timestepping.cc
+ 45
− 0
View file @ eaf14ab5
......@@ -4,6 +4,51 @@
#include "timestepping.hh"
template <class Vector, class Matrix, class Function, size_t dim>
TimeSteppingScheme<Vector, Matrix, Function, dim>::TimeSteppingScheme(
Matrices<Matrix> const &_matrices, Vector const &_u_initial,
Vector const &_v_initial, Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: matrices(_matrices),
u(_u_initial),
v(_v_initial),
a(_a_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class Vector, class Matrix, class Function, size_t dim>
void TimeSteppingScheme<Vector, Matrix, Function, dim>::nextTimeStep() {
u_o = u;
v_o = v;
a_o = a;
postProcessCalled = false;
}
template <class Vector, class Matrix, class Function, size_t dim>
void TimeSteppingScheme<Vector, Matrix, Function, dim>::extractDisplacement(
Vector &displacement) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
template <class Vector, class Matrix, class Function, size_t dim>
void TimeSteppingScheme<Vector, Matrix, Function, dim>::extractVelocity(
Vector &velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = v;
}
template <class Vector, class Matrix, class Function, size_t dim>
void TimeSteppingScheme<Vector, Matrix, Function, dim>::extractOldVelocity(
Vector &oldVelocity) const {
oldVelocity = v_o;
}
#include "timestepping/backward_euler.cc"
#include "timestepping/newmark.cc"
......
src/timestepping.hh
+ 24
− 5
View file @ eaf14ab5
......@@ -10,18 +10,37 @@
template <class Vector, class Matrix, class Function, size_t dim>
class TimeSteppingScheme {
protected:
TimeSteppingScheme(Matrices<Matrix> const &_matrices,
Vector const &_u_initial, Vector const &_v_initial,
Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
public:
void virtual nextTimeStep() = 0;
void nextTimeStep();
void virtual setup(Vector const &ell, double _tau, double relativeTime,
Vector &rhs, Vector &iterate, Matrix &AB) = 0;
void virtual postProcess(Vector const &iterate) = 0;
void virtual extractDisplacement(Vector &displacement) const = 0;
void virtual extractVelocity(Vector &velocity) const = 0;
void virtual extractOldVelocity(Vector &velocity) const = 0;
void extractDisplacement(Vector &displacement) const;
void extractVelocity(Vector &velocity) const;
void extractOldVelocity(Vector &velocity) const;
std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function,
dim>> virtual clone() const = 0;
protected:
Matrices<Matrix> const &matrices;
Vector u, v, a;
Dune::BitSetVector<dim> const &dirichletNodes;
Function const &dirichletFunction;
double dirichletValue;
Vector u_o, v_o, a_o;
double tau;
bool postProcessCalled = true;
};
#include "timestepping/newmark.hh"
......@@ -42,7 +61,7 @@ initTimeStepper(Config::scheme scheme,
case Config::BackwardEuler:
return std::make_shared<
BackwardEuler<Vector, Matrix, Function, dimension>>(
matrices, u_initial, v_initial, velocityDirichletNodes,
matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
velocityDirichletFunction);
default:
assert(false);
......
src/timestepping/backward_euler.cc
+ 29
− 56
View file @ eaf14ab5
......@@ -3,27 +3,19 @@
template <class Vector, class Matrix, class Function, size_t dim>
BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
Matrices<Matrix> const &_matrices, Vector const &_u_initial,
Vector const &_v_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
Vector const &_v_initial, Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: matrices(_matrices),
u(_u_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::nextTimeStep() {
v_o = v;
u_o = u;
}
: TimeSteppingScheme<Vector, Matrix, Function, dim>(
_matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
_dirichletFunction) {}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::setup(
Vector const &ell, double _tau, double relativeTime, Vector &rhs,
Vector &iterate, Matrix &AM) {
postProcessCalled = false;
dirichletFunction.evaluate(relativeTime, dirichletValue);
tau = _tau;
this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
this->tau = _tau;
/* We start out with the formulation
......@@ -48,66 +40,47 @@ void BackwardEuler<Vector, Matrix, Function, dim>::setup(
// set up LHS (for fixed tau, we'd only really have to do this once)
{
Dune::MatrixIndexSet indices(matrices.elasticity.N(),
matrices.elasticity.M());
indices.import(matrices.elasticity);
indices.import(matrices.mass);
indices.import(matrices.damping);
Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
this->matrices.elasticity.M());
indices.import(this->matrices.elasticity);
indices.import(this->matrices.mass);
indices.import(this->matrices.damping);
indices.exportIdx(AM);
}
AM = 0.0;
Arithmetic::addProduct(AM, 1.0 / tau, matrices.mass);
Arithmetic::addProduct(AM, 1.0, matrices.damping);
Arithmetic::addProduct(AM, tau, matrices.elasticity);
Arithmetic::addProduct(AM, 1.0 / this->tau, this->matrices.mass);
Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
Arithmetic::addProduct(AM, this->tau, this->matrices.elasticity);
// set up RHS
{
rhs = ell;
Arithmetic::addProduct(rhs, 1.0 / tau, matrices.mass, v_o);
Arithmetic::subtractProduct(rhs, matrices.elasticity, u_o);
Arithmetic::addProduct(rhs, 1.0 / this->tau, this->matrices.mass,
this->v_o);
Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
}
iterate = v_o;
iterate = this->v_o;
for (size_t i = 0; i < dirichletNodes.size(); ++i)
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
for (size_t j = 0; j < dim; ++j)
if (dirichletNodes[i][j])
iterate[i][j] = (j == 0) ? dirichletValue : 0;
if (this->dirichletNodes[i][j])
iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
Vector const &iterate) {
postProcessCalled = true;
this->postProcessCalled = true;
v = iterate;
this->v = iterate;
u = u_o;
Arithmetic::addProduct(u, tau, v);
}
this->u = this->u_o;
Arithmetic::addProduct(this->u, this->tau, this->v);
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractDisplacement(
Vector &displacement) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractVelocity(
Vector &velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = v;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractOldVelocity(
Vector &velocity) const {
velocity = v_o;
this->a = this->v;
this->a -= this->v_o;
this->a /= this->tau;
}
template <class Vector, class Matrix, class Function, size_t dim>
......
src/timestepping/backward_euler.hh
+ 1
− 18
View file @ eaf14ab5
......@@ -5,34 +5,17 @@ template <class Vector, class Matrix, class Function, size_t dim>
class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
public:
BackwardEuler(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
Vector const &_v_initial,
Vector const &_v_initial, Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
void nextTimeStep() override;
void setup(Vector const &, double, double, Vector &, Vector &,
Matrix &) override;
void postProcess(Vector const &) override;
void extractDisplacement(Vector &) const override;
void extractVelocity(Vector &) const override;
void extractOldVelocity(Vector &) const override;
std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
const;
private:
Matrices<Matrix> const &matrices;
Vector u;
Vector v;
Dune::BitSetVector<dim> const &dirichletNodes;
Function const &dirichletFunction;
double dirichletValue;
Vector u_o;
Vector v_o;
double tau;
bool postProcessCalled = true;
};
#endif
src/timestepping/newmark.cc
+ 32
− 65
View file @ eaf14ab5
......@@ -6,19 +6,9 @@ Newmark<Vector, Matrix, Function, dim>::Newmark(
Vector const &_v_initial, Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: matrices(_matrices),
u(_u_initial),
v(_v_initial),
a(_a_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::nextTimeStep() {
a_o = a;
v_o = v;
u_o = u;
}
: TimeSteppingScheme<Vector, Matrix, Function, dim>(
_matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
_dirichletFunction) {}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
......@@ -26,9 +16,8 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
double relativeTime,
Vector &rhs, Vector &iterate,
Matrix &AM) {
postProcessCalled = false;
dirichletFunction.evaluate(relativeTime, dirichletValue);
tau = _tau;
this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
this->tau = _tau;
/* We start out with the formulation
......@@ -54,76 +43,54 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
// set up LHS (for fixed tau, we'd only really have to do this once)
{
Dune::MatrixIndexSet indices(matrices.elasticity.N(),
matrices.elasticity.M());
indices.import(matrices.elasticity);
indices.import(matrices.mass);
indices.import(matrices.damping);
Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
this->matrices.elasticity.M());
indices.import(this->matrices.elasticity);
indices.import(this->matrices.mass);
indices.import(this->matrices.damping);
indices.exportIdx(AM);
}
AM = 0.0;
Arithmetic::addProduct(AM, 2.0 / tau, matrices.mass);
Arithmetic::addProduct(AM, 1.0, matrices.damping);
Arithmetic::addProduct(AM, tau / 2.0, matrices.elasticity);
Arithmetic::addProduct(AM, 2.0 / this->tau, this->matrices.mass);
Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
Arithmetic::addProduct(AM, this->tau / 2.0, this->matrices.elasticity);
// set up RHS
{
rhs = ell;
Arithmetic::addProduct(rhs, 2.0 / tau, matrices.mass, v_o);
Arithmetic::addProduct(rhs, matrices.mass, a_o);
Arithmetic::subtractProduct(rhs, tau / 2.0, matrices.elasticity, v_o);
Arithmetic::subtractProduct(rhs, matrices.elasticity, u_o);
Arithmetic::addProduct(rhs, 2.0 / this->tau, this->matrices.mass,
this->v_o);
Arithmetic::addProduct(rhs, this->matrices.mass, this->a_o);
Arithmetic::subtractProduct(rhs, this->tau / 2.0, this->matrices.elasticity,
this->v_o);
Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
}
iterate = v_o;
iterate = this->v_o;
for (size_t i = 0; i < dirichletNodes.size(); ++i)
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
for (size_t j = 0; j < dim; ++j)
if (dirichletNodes[i][j])
iterate[i][j] = (j == 0) ? dirichletValue : 0;
if (this->dirichletNodes[i][j])
iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::postProcess(
Vector const &iterate) {
postProcessCalled = true;
this->postProcessCalled = true;
v = iterate;
this->v = iterate;
// u1 = tau/2 ( v1 + v0 ) + u0
u = u_o;
Arithmetic::addProduct(u, tau / 2.0, v);
Arithmetic::addProduct(u, tau / 2.0, v_o);
this->u = this->u_o;
Arithmetic::addProduct(this->u, this->tau / 2.0, this->v);
Arithmetic::addProduct(this->u, this->tau / 2.0, this->v_o);
// a1 = 2/tau ( v1 - v0 ) - a0
a = 0;
Arithmetic::addProduct(a, 2.0 / tau, v);
Arithmetic::subtractProduct(a, 2.0 / tau, v_o);
Arithmetic::subtractProduct(a, 1.0, a_o);
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractDisplacement(
Vector &displacement) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u;
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractVelocity(
Vector &velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = v;
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractOldVelocity(
Vector &velocity) const {
velocity = v_o;
this->a = 0.0;
Arithmetic::addProduct(this->a, 2.0 / this->tau, this->v);
Arithmetic::subtractProduct(this->a, 2.0 / this->tau, this->v_o);
Arithmetic::subtractProduct(this->a, 1.0, this->a_o);
}
template <class Vector, class Matrix, class Function, size_t dim>
......
src/timestepping/newmark.hh
+ 0
− 21
View file @ eaf14ab5
......@@ -9,32 +9,11 @@ class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
void nextTimeStep() override;
void setup(Vector const &, double, double, Vector &, Vector &,
Matrix &) override;
void postProcess(Vector const &) override;
void extractDisplacement(Vector &) const override;
void extractVelocity(Vector &) const override;
void extractOldVelocity(Vector &) const override;
std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
const;
private:
Matrices<Matrix> const &matrices;
Vector u;
Vector v;
Vector a;
Dune::BitSetVector<dim> const &dirichletNodes;
Function const &dirichletFunction;
double dirichletValue;
Vector u_o;
Vector v_o;
Vector a_o;
double tau;
bool postProcessCalled = true;
};
#endif
src/timestepping_tmpl.cc
+ 1
− 0
View file @ eaf14ab5
......@@ -8,5 +8,6 @@
using Function = Dune::VirtualFunction<double, double>;
template class TimeSteppingScheme<Vector, Matrix, Function, MY_DIM>;
template class Newmark<Vector, Matrix, Function, MY_DIM>;
template class BackwardEuler<Vector, Matrix, Function, MY_DIM>;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment