Once inverting potential-field geophysical data an appropriate discretization of the model is required to accurately construct complicated geometries of the causative sources. Rectangular prisms (structured meshes) have limitations to recover and preserve the edges of real geological sources. Here an isoparametric finite-element (FE) methodology is used to design an unstructured mesh for use in 3D inverse modeling of gravity data. The calculation of the sensitivity kernel of the forward operator uses Gauss-Legendre quadrature rather than the analytic formulation. For the sake of instability of inversion operator in gravity data and to solve the Tikhonov norms term, the least-squares QR-factorization (LSQR) technique is used. This method is accompanied with weighted generalized cross-validation (WGCV) for selecting the optimum regularization parameter value. The depth weighting function is also incorporated in the formulation of the objective function to suppress the impact of shallow features and recover sources at an appropriate depth. The proposed algorithm was applied to noise-corrupted synthetic data along with a real case study where a gravity survey was used for iron exploration in Yazd province, central Iran. The obtained results of synthetic example indicated that the proposed 3D inversion method recovers the geometry and density contrast values which are similar to the true structures and its application on a real example data set recovers geologically reasonable complex structures. So, this modified algorithm is a type of smooth inversion for boundary detection and whole understanding of the physical property distribution of the subsurface.