  We describe a computational framework for simulating suspensions of rigid
particles in Newtonian Stokes flow. One central building block is a
collision-resolution algorithm that overcomes the numerical constraints arising
from particle collisions. This algorithm extends the well-known complementarity
method for non-smooth multi-body dynamics to resolve collisions in dense rigid
body suspensions. This approach formulates the collision resolution problem as
a linear complementarity problem with geometric `non-overlapping' constraints
imposed at each timestep. It is then reformulated as a constrained quadratic
programming problem and the Barzilai-Borwein projected gradient descent method
is applied for its solution. This framework is designed to be applicable for
any convex particle shape, e.g., spheres and spherocylinders, and applicable to
any Stokes mobility solver, including the Rotne-Prager-Yamakawa approximation,
Stokesian Dynamics, and PDE solvers (e.g., boundary integral and immersed
boundary methods). In particular, this method imposes Newton's Third Law and
records the entire contact network. Further, we describe a fast, parallel, and
spectrally-accurate boundary integral method tailored for spherical particles,
capable of resolving lubrication effects. We show weak and strong parallel
scalings up to $8\times 10^4$ particles with approximately $4\times 10^7$
degrees of freedom on $1792$ cores. We demonstrate the versatility of this
framework with several examples, including sedimentation of particle clusters,
and active matter systems composed of ensembles of particles driven to rotate.
