Inertia

Time discretisation: 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> fres = 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.0) * a;

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

    // new acceleration
    M.solve(fres, a);

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