In this work, a fast multipole boundary element method for 3D elasticity problem was developed by the application of the fast multipole algorithm and isoparametric 8-node boundary elements with quadratic shape functions. The problem is described by the boundary integral equation involving the Kelvin solutions. In order to keep the numerical integration error on appropriate level, an adaptive method with subdivision of boundary elements into subelements, described in the literature, was applied. An extension of the neighbour list of boundary element clusters, corresponding to near-field computations, was proposed in order to reduce the truncation error of expansions in problems with high stress concentration. Efficiency of the method is illustrated by numerical examples including a solid with single spherical cavity, solids with two interacting spherical cavities, and numerical homogenization of solids with cubic arrangement of spherical cavities. All results agree with analytical models available in the literature. The examples show that the method can be applied to the analysis of porous structures.