Minimum-structure inversion is one of the most effective tools for the inversion of gravity data. However, the standard Gauss-Newton algorithms that are commonly used for the minimization procedure and that employ forward solvers based on analytic formulas require large memory storage for the formation and inversion of the involved matrices. An alternative to the analytical solvers are numerical ones that result in sparse matrices. This sparsity suits gradient-based minimization methods that avoid the explicit formation of the inversion matrices and that solve the system of equations using memory-efficient iterative techniques. We have developed several numerical schemes for the forward modeling of gravity data using the finite-element and finite-volume methods on unstructured grids. In the finite-volume method, a Delaunay tetrahedral grid and its dual Voronoï grid are used to find the primary solution (i.e., gravitational potential) at the centers and vertices of the tetrahedra, respectively (cell-centered and vertex-centered schemes). In the finite-element method, Delaunay tetrahedral grids are used to develop linear and quadratic finite-element schemes. Different techniques are used to recover the vertical component of gravitational acceleration from the gravitational potential. In the finite-volume scheme, a differencing method is used; in the finite-element method, basis functions are used. The capabilities of the finite-volume and finite-element schemes were tested on simple and realistic synthetic examples. The results showed that the quadratic finite-element scheme is the most accurate but also the most computationally demanding scheme. The best trade-offs between accuracy and computational resource requirement were achieved by the linear finite-element and vertex-centered finite-volume schemes.