A new singularity-free analytical formula has been developed for the gravity field of arbitrary 3D polyhedral mass bodies with horizontally and vertically varying density contrast using third-order polynomial functions. First, the observation sites are moved to the origin of the coordinate system. Then, the volume and surface integral theorems are invoked successively to transform the volume integrals into surface integrals over polygonal faces and into line integrals over the edges of the polyhedral mass bodies. Furthermore, singularity-free closed-form solutions are derived for these line integrals over the edges. Thus, the observation sites can be located inside, on, or outside the 3D distributions. A synthetic prismatic mass body is adopted to verify the accuracy and singularity-free property of our newly developed analytical expressions. Excellent agreements are obtained between our solutions and other published closed-form solutions with relative errors in the order of ${10}^{-12}\%$ to ${10}^{-13}\%$. In addition, an octahedral model and a near-Earth asteroid model are used to verify the accuracy of the presented method for complicated target structures by comparing the results with those from a high-order Gaussian quadrature approach.