BEM–FEM coupling is desirable for three-dimensional problems involving specific features such as (i) large or unbounded media with linear constitutive properties, (ii) cracks, (iii) critical parts of complex geometry requiring accurate stress analyses. However, for cases with a BEM discretization involving a large number NBEM of degrees of freedom, setting up the BEM contribution to the coupled problem using conventional techniques is an expensive $O(NBEM^{2})$ task. Moreover, the fully-populated BEM block entails a $O(NBEM^{2})$ storage requirement and a $O(NBEM^{3})$ contribution to the solution time via usual direct solvers. To overcome these pitfalls, the BEM contribution is formulated using the fast multipole method (FMM) and the coupl...