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; }