Abstract This paper reports an indirect boundary element method in elasticity that is most suitable to deal with particulate solids. The method involves a distribution of a double layer potential and, after a suitable completion and deflation, is amenable to iterative solution techniques. It can therefore accommodate a large number of particles with complex geometries. Convergence of the method is significantly improved by the introduction of a simple domain decomposition to solve the system of equations. The method is illustrated by the translating sphere problem, the load transfer problem between two spheres at near contact, and the shear deformation of a cluster of 125 spheres initially in a simple cubic array.