Harmonic density interpolation methods for high-order evaluation of Laplace layer potentials in 2D and 3D