<p>The gravitational effect of tesseroid, expressed as a Newtonian volume integral, is generally not analytically solvable across most of the globe, in spherical coordinates. A closed-form solution exists only when the computation point is on the polar axis. In this case, spherical and spherical polar coordinates coincide, meaning that a tesseroid regularly defined in spherical coordinates remains regular in spherical polar coordinates. Thus, the analytical solvability of the tesseroid’s gravitational effect depends on whether its domain remains regular in spherical polar coordinates. This study decomposes a tesseroid with an irregular domain in spherical polar coordinates into a set of regular triangular prism elements and derives an analytical solution for the gravitational effect of a triangular prism. This allows us to obtain an analytical solution for a tesseroid at any computation point. To validate this decomposition, we examine the linear relationships along each tesseroid edge and demonstrate that it can be well approximated by four regular triangular prism elements. To assess accuracy, the relative errors of the analytical solution are compared globally against two reference methods: the 3D Gauss–Legendre quadrature (3D-GLQ) and the closed-form solution at the polar axis. Results show that, except near the computation point and its symmetric point located on the opposite side of the earth, our analytical solution surpasses the 3D-GLQ in accuracy and achieves 26.2% of computational efficiency compared to 3D-GLQ with order 7. When the computation point is located on the polar axis, the linear relationships hold exactly, and our solution yields the true value. Finally, the hybrid model where 3D-GLQ replaces our analytical solution in lower-accuracy regions (near the computation point and symmetric point) surpasses the pure 3D-GLQ method, and the superposition error elimination effect occurs in the hybrid model.</p>