Newer
Older
template <class Vector, class Matrix, class Function, size_t dim>
BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
Matrix const &_A, Matrix const &_M, Matrix const &_C,
Vector const &_u_initial, Vector const &_ur_initial,
Vector const &_v_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::nextTimeStep() {
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;
postProcessRelativeQuantitiesCalled = false;
dirichletFunction.evaluate(relativeTime, dirichletValue);
/* We start out with the formulation
M a + C v + A u = ell
Backward Euler means
a1 = 1.0/tau ( v1 - v0 )
u1 = tau v1 + u0
in summary, we get at time t=1
M [1.0/tau ( v1 - v0 )] + C v1
+ A [tau v1 + u0] = ell
or
1.0/tau M v1 + C v1 + tau A v1
= [1.0/tau M + C + tau A] v1
= ell + 1.0/tau M v0 - A u0
*/
// set up LHS (for fixed tau, we'd only really have to do this once)
{
Dune::MatrixIndexSet indices(A.N(), A.M());
indices.import(A);
indices.import(M);
indices.import(C);
Arithmetic::addProduct(AM, 1.0 / tau, M);
Arithmetic::addProduct(AM, 1.0, C);
Arithmetic::addProduct(rhs, 1.0 / tau, M, v_o);
Arithmetic::subtractProduct(rhs, A, u_o);
for (size_t i = 0; i < dirichletNodes.size(); ++i)
for (size_t j = 0; j < dim; ++j)
if (dirichletNodes[i][j])
iterate[i][j] = (j == 0) ? dirichletValue : 0;
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
Vector const &iterate) {
postProcessCalled = true;
Arithmetic::addProduct(u, tau, v);
}
template <class Vector, class Matrix, class Function, size_t dim>
void
BackwardEuler<Vector, Matrix, Function, dim>::postProcessRelativeQuantities() {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
postProcessRelativeQuantitiesCalled = true;
vr = v;
for (size_t i = 0; i < dirichletNodes.size(); ++i)
vr[i][0] -= dirichletValue;
ur = ur_o;
Arithmetic::addProduct(ur, tau, vr);
}
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>::extractRelativeDisplacement(
Vector &displacement) const {
if (!postProcessRelativeQuantitiesCalled)
DUNE_THROW(Dune::Exception,
"It seems you forgot to call postProcessRelativeQuantities!");
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>::extractRelativeVelocity(
Vector &velocity) const {
if (!postProcessRelativeQuantitiesCalled)
DUNE_THROW(Dune::Exception,
"It seems you forgot to call postProcessRelativeQuantities!");
velocity = vr;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractOldVelocity(
Vector &velocity) const {
velocity = v_o;
}