Inertia + damping

Velocity Verlet

// position, velocity, acceleration (and history: last increment)
xt::xtensor<double,2> u   = xt::zeros<double>({nnode, ndim});
xt::xtensor<double,2> v   = xt::zeros<double>({nnode, ndim});
xt::xtensor<double,2> a   = xt::zeros<double>({nnode, ndim});
xt::xtensor<double,2> v_n = xt::zeros<double>({nnode, ndim});
xt::xtensor<double,2> a_n = xt::zeros<double>({nnode, ndim});

// residual force
xt::xtensor<double,2> fr = xt::zeros<double>({nnode, ndim});

// compute mass matrix
// (often assumed constant & diagonal, remove either assumption if needed)
GooseFEM::MatrixDiagonal M(...);

...

// time increments
for ( ... )
{
  // store history
  xt::noalias(v_n) = v;
  xt::noalias(a_n) = a;

  // new displacement
  xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt,2.) * a;

  // update residual force (and mass matrix if needed)
  ...

  // estimate 1: new velocity
  xt::noalias(v) = v_n + dt * a_n;

  // estimate 1: new residual force (and mass matrix if needed)
  ...

  // estimate 1: new acceleration
  M.solve(fr, a);

  // estimate 2: new velocity
  xt::noalias(v) = v_n + .5 * dt * ( a_n + a );

  // estimate 2: new residual force (and mass matrix if needed)
  ...

  // estimate 2: new acceleration
  M.solve(fr, a);

  // new velocity
  xt::noalias(v) = v_n + .5 * dt * ( a_n + a );

  // new residual force (and mass matrix if needed)
  ...

  // new acceleration
  M.solve(fr, a);
}