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

[Cleanup] Split up timestepping into multiple headers

parent 72cee6fa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/timestepping.cc
+ 4
− 287
View file @ 4e6db7b0
...@@ -4,294 +4,11 @@ ...@@ -4,294 +4,11 @@
#include <dune/istl/matrixindexset.hh> #include <dune/istl/matrixindexset.hh>
#include <dune/fufem/arithmetic.hh> #include <dune/fufem/arithmetic.hh>
#include "timestepping.hh"
template <class VectorType, class MatrixType, class FunctionType, int dim>
ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::ImplicitEuler(
MatrixType const &_A, VectorType const &_u_initial,
VectorType const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
u(_u_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
v_old = v;
u_old = u;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::setup(
VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
VectorType &problem_iterate, MatrixType &problem_AB) {
postProcessCalled = false;
tau = _tau;
problem_rhs = ell;
Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
// For fixed tau, we'd only really have to do this once
problem_AB = A;
problem_AB *= tau;
// v_old makes a good initial iterate; we could use anything, though
problem_iterate = v_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 ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::postProcess(
VectorType const &problem_iterate) {
postProcessCalled = true;
v = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau, v);
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<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 ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
VectorType &velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = v;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
VectorType const &_v_initial, VectorType const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
B(_B),
u(_u_initial),
v(_v_initial),
a(_a_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
a_old = a;
v_old = v;
u_old = u;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
VectorType &problem_iterate, MatrixType &problem_AB) {
postProcessCalled = false;
tau = _tau;
problem_rhs = ell;
Arithmetic::addProduct(problem_rhs, 2.0 / tau, B, v_old);
Arithmetic::addProduct(problem_rhs, B, a_old);
Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
Arithmetic::addProduct(problem_rhs, -tau / 2.0, A, v_old);
// 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_AB);
problem_AB = 0.0;
Arithmetic::addProduct(problem_AB, tau / 2.0, A);
Arithmetic::addProduct(problem_AB, 2.0 / tau, B);
// v_old makes a good initial iterate; we could use anything, though
problem_iterate = 0.0;
for (size_t i = 0; i < dirichletNodes.size(); ++i) #include "timestepping.hh"
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 Newmark<VectorType, MatrixType, FunctionType, dim>::postProcess(
VectorType const &problem_iterate) {
postProcessCalled = true;
v = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau / 2.0, v);
Arithmetic::addProduct(u, tau / 2.0, v_old);
a = 0;
Arithmetic::addProduct(a, 2.0 / tau, v);
Arithmetic::addProduct(a, -2.0 / tau, v_old);
Arithmetic::addProduct(a, -1.0, a_old);
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<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 Newmark<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
VectorType &velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = v;
}
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 &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
B(_B),
u(_u_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
v_old = v;
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_AB) {
postProcessCalled = false;
tau = _tau;
problem_rhs = ell;
Arithmetic::addProduct(problem_rhs, 1.0 / tau, B, v_old);
Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
// 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_AB);
problem_AB = 0.0;
Arithmetic::addProduct(problem_AB, tau, A);
Arithmetic::addProduct(problem_AB, 1.0 / tau, B);
// v_old makes a good initial iterate; we could use anything, though
problem_iterate = v_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;
v = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau, v);
}
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 = v; #include "timestepping/eulerpair.cc"
} #include "timestepping/impliciteuler.cc"
#include "timestepping/newmark.cc"
#include "timestepping_tmpl.cc" #include "timestepping_tmpl.cc"
src/timestepping.hh
+ 3
− 97
View file @ 4e6db7b0
...@@ -16,102 +16,8 @@ class TimeSteppingScheme { ...@@ -16,102 +16,8 @@ class TimeSteppingScheme {
void virtual extractVelocity(VectorType &velocity) const = 0; void virtual extractVelocity(VectorType &velocity) const = 0;
}; };
template <class VectorType, class MatrixType, class FunctionType, int dim> #include "timestepping/impliciteuler.hh"
class ImplicitEuler #include "timestepping/newmark.hh"
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> { #include "timestepping/eulerpair.hh"
public:
ImplicitEuler(MatrixType const &_A, VectorType const &_u_initial,
VectorType const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
void virtual nextTimeStep() override;
void virtual setup(VectorType const &, double, double, VectorType &,
VectorType &, MatrixType &) override;
void virtual postProcess(VectorType const &) override;
void virtual extractDisplacement(VectorType &) const override;
void virtual extractVelocity(VectorType &) const override;
private:
MatrixType const &A;
VectorType u;
VectorType v;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType v_old;
double tau;
bool postProcessCalled = false;
};
template <class VectorType, class MatrixType, class FunctionType, int dim>
class Newmark
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
Newmark(MatrixType const &_A, MatrixType const &_B,
VectorType const &_u_initial, VectorType const &_v_initial,
VectorType const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
void virtual nextTimeStep() override;
void virtual setup(VectorType const &, double, double, VectorType &,
VectorType &, MatrixType &) override;
void virtual postProcess(VectorType const &) override;
void virtual extractDisplacement(VectorType &) const override;
void virtual extractVelocity(VectorType &) const override;
private:
MatrixType const &A;
MatrixType const &B;
VectorType u;
VectorType v;
VectorType a;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType v_old;
VectorType a_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 &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
void virtual nextTimeStep() override;
void virtual setup(VectorType const &, double, double, VectorType &,
VectorType &, MatrixType &) override;
void virtual postProcess(VectorType const &) override;
void virtual extractDisplacement(VectorType &) const override;
void virtual extractVelocity(VectorType &) const override;
private:
MatrixType const &A;
MatrixType const &B;
VectorType u;
VectorType v;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType v_old;
double tau;
bool postProcessCalled = false;
};
#endif #endif
src/timestepping/eulerpair.cc 0 → 100644
+ 95
− 0
View file @ 4e6db7b0
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 &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
B(_B),
u(_u_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
v_old = v;
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_AB) {
postProcessCalled = false;
tau = _tau;
problem_rhs = ell;
Arithmetic::addProduct(problem_rhs, 1.0 / tau, B, v_old);
Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
// 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_AB);
problem_AB = 0.0;
Arithmetic::addProduct(problem_AB, tau, A);
Arithmetic::addProduct(problem_AB, 1.0 / tau, B);
// v_old makes a good initial iterate; we could use anything, though
problem_iterate = v_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;
v = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau, v);
}
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 = v;
}
src/timestepping/eulerpair.hh 0 → 100644
+ 35
− 0
View file @ 4e6db7b0
#ifndef DUNE_TECTONIC_TIMESTEPPING_EULERPAIR_HH
#define DUNE_TECTONIC_TIMESTEPPING_EULERPAIR_HH
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 &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
void virtual nextTimeStep() override;
void virtual setup(VectorType const &, double, double, VectorType &,
VectorType &, MatrixType &) override;
void virtual postProcess(VectorType const &) override;
void virtual extractDisplacement(VectorType &) const override;
void virtual extractVelocity(VectorType &) const override;
private:
MatrixType const &A;
MatrixType const &B;
VectorType u;
VectorType v;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType v_old;
double tau;
bool postProcessCalled = false;
};
#endif
src/timestepping/impliciteuler.cc 0 → 100644
+ 87
− 0
View file @ 4e6db7b0
template <class VectorType, class MatrixType, class FunctionType, int dim>
ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::ImplicitEuler(
MatrixType const &_A, VectorType const &_u_initial,
VectorType const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
u(_u_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
v_old = v;
u_old = u;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::setup(
VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
VectorType &problem_iterate, MatrixType &problem_AB) {
postProcessCalled = false;
tau = _tau;
problem_rhs = ell;
Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
// For fixed tau, we'd only really have to do this once
problem_AB = A;
problem_AB *= tau;
// v_old makes a good initial iterate; we could use anything, though
problem_iterate = v_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 ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::postProcess(
VectorType const &problem_iterate) {
postProcessCalled = true;
v = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau, v);
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void ImplicitEuler<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 ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
VectorType &velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = v;
}
src/timestepping/impliciteuler.hh 0 → 100644
+ 34
− 0
View file @ 4e6db7b0
#ifndef DUNE_TECTONIC_TIMESTEPPING_IMPLICITEULER_HH
#define DUNE_TECTONIC_TIMESTEPPING_IMPLICITEULER_HH
template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitEuler
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
ImplicitEuler(MatrixType const &_A, VectorType const &_u_initial,
VectorType const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
void virtual nextTimeStep() override;
void virtual setup(VectorType const &, double, double, VectorType &,
VectorType &, MatrixType &) override;
void virtual postProcess(VectorType const &) override;
void virtual extractDisplacement(VectorType &) const override;
void virtual extractVelocity(VectorType &) const override;
private:
MatrixType const &A;
VectorType u;
VectorType v;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType v_old;
double tau;
bool postProcessCalled = false;
};
#endif
src/timestepping/newmark.cc 0 → 100644
+ 104
− 0
View file @ 4e6db7b0
template <class VectorType, class MatrixType, class FunctionType, int dim>
Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
VectorType const &_v_initial, VectorType const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction)
: A(_A),
B(_B),
u(_u_initial),
v(_v_initial),
a(_a_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
a_old = a;
v_old = v;
u_old = u;
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
VectorType &problem_iterate, MatrixType &problem_AB) {
postProcessCalled = false;
tau = _tau;
problem_rhs = ell;
Arithmetic::addProduct(problem_rhs, 2.0 / tau, B, v_old);
Arithmetic::addProduct(problem_rhs, B, a_old);
Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
Arithmetic::addProduct(problem_rhs, -tau / 2.0, A, v_old);
// 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_AB);
problem_AB = 0.0;
Arithmetic::addProduct(problem_AB, tau / 2.0, A);
Arithmetic::addProduct(problem_AB, 2.0 / tau, B);
// v_old makes a good initial iterate; we could use anything, though
problem_iterate = 0.0;
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 Newmark<VectorType, MatrixType, FunctionType, dim>::postProcess(
VectorType const &problem_iterate) {
postProcessCalled = true;
v = problem_iterate;
u = u_old;
Arithmetic::addProduct(u, tau / 2.0, v);
Arithmetic::addProduct(u, tau / 2.0, v_old);
a = 0;
Arithmetic::addProduct(a, 2.0 / tau, v);
Arithmetic::addProduct(a, -2.0 / tau, v_old);
Arithmetic::addProduct(a, -1.0, a_old);
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
void Newmark<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 Newmark<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
VectorType &velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
velocity = v;
}
src/timestepping/newmark.hh 0 → 100644
+ 38
− 0
View file @ 4e6db7b0
#ifndef DUNE_TECTONIC_TIMESTEPPING_NEWMARK_HH
#define DUNE_TECTONIC_TIMESTEPPING_NEWMARK_HH
template <class VectorType, class MatrixType, class FunctionType, int dim>
class Newmark
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
Newmark(MatrixType const &_A, MatrixType const &_B,
VectorType const &_u_initial, VectorType const &_v_initial,
VectorType const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
void virtual nextTimeStep() override;
void virtual setup(VectorType const &, double, double, VectorType &,
VectorType &, MatrixType &) override;
void virtual postProcess(VectorType const &) override;
void virtual extractDisplacement(VectorType &) const override;
void virtual extractVelocity(VectorType &) const override;
private:
MatrixType const &A;
MatrixType const &B;
VectorType u;
VectorType v;
VectorType a;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
VectorType u_old;
VectorType v_old;
VectorType a_old;
double tau;
bool postProcessCalled = false;
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment