template <class VectorType, class MatrixType, class FunctionType, int dim>
Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
    MatrixType const &_A, MatrixType const &_M, VectorType const &_u_initial,
    VectorType const &_v_initial, VectorType const &_a_initial,
    Dune::BitSetVector<dim> const &_dirichletNodes,
    FunctionType const &_dirichletFunction)
    : A(_A),
      M(_M),
      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_o = a;
  v_o = v;
  u_o = 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_AM) {
  postProcessCalled = false;

  tau = _tau;

  MatrixType const &C = A;
  double const wc = 0.0;
  /* We start out with the formulation

       M a + C v + A u = ell

     Newmark means

       a1 = 2/tau ( v1 - v0 ) - a0
       u1 = tau/2 ( v1 + v0 ) + u0

     in summary, we get at time t=1

       M [2/tau ( u1 - u0 ) - a0] + C v1
       + A [tau/2 ( v1 + v0 ) + u0] = ell

     or

       2/tau M v1 + C v1 + tau/2 A v1
       = [2/tau M + C + tau/2 A] v1
       = ell + 2/tau M v0 + M a0
       - tau/2 A 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(problem_AM);
  }
  problem_AM = 0.0;
  Arithmetic::addProduct(problem_AM, (1.0 - wc) * 2.0 / tau, M);
  Arithmetic::addProduct(problem_AM, wc, C);
  Arithmetic::addProduct(problem_AM, tau / 2.0, A);

  // set up RHS
  {
    problem_rhs = ell;
    Arithmetic::addProduct(problem_rhs, (1.0 - wc) * 2.0 / tau, M, v_o);
    Arithmetic::addProduct(problem_rhs, 1.0 - wc, M, a_o);
    Arithmetic::subtractProduct(problem_rhs, tau / 2.0, A, v_o);
    Arithmetic::subtractProduct(problem_rhs, A, u_o);
  }

  // v_o 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;

  // u1 = tau/2 ( v1 + v0 ) + u0
  u = u_o;
  Arithmetic::addProduct(u, tau / 2.0, v);
  Arithmetic::addProduct(u, tau / 2.0, v_o);

  // a1 = 2/tau ( v1 - v0 ) - a0
  a = 0;
  Arithmetic::addProduct(a, 2.0 / tau, v);
  Arithmetic::subtractProduct(a, 2.0 / tau, v_o);
  Arithmetic::subtractProduct(a, 1.0, a_o);
}

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