Newer
Older
#include <dune/solvers/common/arithmetic.hh>
#include "backward_euler.hh"
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, Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: RateUpdater<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) {
this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
this->tau = _tau;
/* 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(this->matrices.elasticity.N(),
this->matrices.elasticity.M());
indices.import(this->matrices.elasticity);
indices.import(this->matrices.mass);
indices.import(this->matrices.damping);
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);
Arithmetic::addProduct(rhs, 1.0 / this->tau, this->matrices.mass,
this->v_o);
Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
iterate = this->v_o;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
for (size_t j = 0; j < dim; ++j)
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) {
this->postProcessCalled = true;
this->u = this->u_o;
Arithmetic::addProduct(this->u, this->tau, this->v);
this->a = this->v;
this->a -= this->v_o;
this->a /= this->tau;
template <class Vector, class Matrix, class Function, size_t dim>
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>>
BackwardEuler<Vector, Matrix, Function, dim>::clone() const {
return std::make_shared<BackwardEuler<Vector, Matrix, Function, dim>>(*this);
}