A generalization of a recently introduced recursive numerical method for the exact evaluation of integrals of regular solid harmonics and their normal derivatives over simplex elements in ℝ3 is presented. The original Quadrature to Expansion (Q2X) method achieves optimal per-element asymptotic complexity, however, it considered only constant density functions over the elements. Here, we generalize this method to support arbitrary degree polynomial density functions, which is achieved in an extended recursive framework while maintaining the optimality of the complexity. The method is derived for 1- and 2- simplex elements in ℝ3 and can be used for the boundary element method and vortex methods coupled with the fast multipole method.