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 &_v_initial, Dune::BitSetVector<dim> const &_dirichletNodes, Function const &_dirichletFunction) : A(_A), M(_M), C(_C), 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; } 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; 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(A.N(), A.M()); indices.import(A); indices.import(M); indices.import(C); indices.exportIdx(AM); } AM = 0.0; Arithmetic::addProduct(AM, 1.0 / tau, M); Arithmetic::addProduct(AM, 1.0, C); Arithmetic::addProduct(AM, tau, A); // set up RHS { rhs = ell; Arithmetic::addProduct(rhs, 1.0 / tau, M, v_o); Arithmetic::subtractProduct(rhs, A, u_o); } // v_o makes a good initial iterate; we could use anything, though iterate = 0.0; for (size_t i = 0; i < dirichletNodes.size(); ++i) for (size_t j = 0; j < dim; ++j) if (dirichletNodes[i][j]) { if (j == 0) dirichletFunction.evaluate(relativeTime, iterate[i][j]); else iterate[i][j] = 0; } } template <class Vector, class Matrix, class Function, size_t dim> void BackwardEuler<Vector, Matrix, Function, dim>::postProcess( Vector const &iterate) { postProcessCalled = true; v = iterate; u = u_o; Arithmetic::addProduct(u, tau, 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; }