Abstract. In a series of two papers, we present numerical, integral-based methods to compute accurately the self-gravitating field and potential induced by tri-dimensional, axially symmetric fluids, with a special regard for tori, discs and rings. This first article is concerned with a fully numerical approach. Complex shapes, small/large aspect ratios, important density gradients and compact/extended systems can be accounted for. Loop singularities in the Poisson integrals are carefully treated from kernel splitting and/or density splitting. Field components are obtained from density splitting: the local density field is separated into a vertically homogeneous contribution for which the integrable singularity is known in a closed form, plus a "residual" contribution which yields a regular integrand. This technique is exact in the vertically homogeneous limit. The potential is computed from double splitting: kernel splitting (to isolate explicitly the singular function), followed by density splitting. In each two directions, numerical quadratures are performed using a dynamical mesh (very efficient for flat and extended systems), combined with a high-order, irregular-spacing scheme. Global performances are demonstrated through a few test-configurations. The accuracy is potentially very high, and can reach the computer precision (for simple geometries and smooth density profiles). In terms of computing time to precision ratio, present methods are asymptotically more efficient than classical finite-difference methods. These can be employed either as a "Poisson-solver" or only to compute a sub-set of Dirichlet/Neumann-type boundary conditions in order to initialize spectral/finite-difference/finite element methods. As an example of an astrophysical application, we determine the structure of an uniformly rotating polytrope belonging to the compressible, Dyson-Wong sequence. In a second paper (Paper II), we show that singularities are integrable analytically for geometrically thin systems where the density is of the form ρ ∝ z n . Improvements and extensions (such as the non-axi-symmetric case) are possible and discussed.