Integral equation methods are extensively used for computational electromagnetism, and accelerated with fast multipole methods or, for problems in layered media, with the fast Fourier transform (FFT). Unfortunately, the efficiency of FFT-based acceleration can be dramatically reduced by the presence of multiscale features. We propose an efficient and robust algorithm to overcome this barrier. A hierarchy of grids of different resolution is used to simultaneously resolve sub-wavelength details and propagate fields efficiently across large distances with the FFT. Integration and pre-correction costs are minimized by adapting projection stencils to the size of each triangle and enabling the use of the quasi-static Green’s function for short distances. Finally, a clever implementation based on sparse matrices exploits empty areas to reduce computational cost and memory consumption. The method is fully automated, and was tested on several structures including layouts of commercial products. Compared to existing algorithms, we demonstrate a speed-up between 7.1 and 19.5 times and a reduction in memory consumption by up to 2.9 times.