  A stable numerical solution of the steady Stokes problem requires
compatibility between the choice of velocity and pressure approximation that
has traditionally proven problematic for meshless methods. In this work, we
present a discretization that couples a staggered scheme for pressure
approximation with a divergence-free velocity reconstruction to obtain an
adaptive, high-order, finite difference-like discretization that can be
efficiently solved with conventional algebraic multigrid techniques. We use
analytic benchmarks to demonstrate equal-order convergence for both velocity
and pressure when solving problems with curvilinear geometries. In order to
study problems in dense suspensions, we couple the solution for the flow to the
equations of motion for freely suspended particles in an implicit monolithic
scheme. The combination of high-order accuracy with fully-implicit schemes
allows the accurate resolution of stiff lubrication forces directly from the
solution of the Stokes problem without the need to introduce sub-grid
lubrication models.
