Efficient algorithms for the simulation of strained heteroepitaxial growth with intermixing in 2+1 dimensions are presented. The first of these algorithms is an extension of the energy localization method [T. P. Schulze and P. Smereka, An energy localization principle and its application to fast kinetic Monte Carlo simulation of heteroepitaxial growth, J. Mech. Phys. Sol., 3 (2009), 521–538] from 1+1 to 2+1 dimensions. Two approximations of this basic algorithm are then introduced, one of which treats adatoms in a more efficient manner, while the other makes use of an approximation of the change in elastic energy in terms of local elastic energy density. In both cases, it is demonstrated that a reasonable level of fidelity is achieved. Results are presented showing how the film morphology is affected by misfit and deposition rate. In addition, simulations of stacked quantum dots are also presented.