<p>Gravity anomalies resulting from continuous density media are governed by a three-dimensional Poisson equation. Current modeling techniques primarily depend on integral solutions to address these anomalies. This approach results in three-dimensional gravity anomalies, including vector and tensor components of gravitational fields, which inevitably involve integral singularities. To address this challenge, we present a finite-element algorithm based on partial differential equations that is designed to solve the boundary value problem of the three-dimensional gravitational potential during the forward modeling calculation of gravity anomalies in continuous density media. In order to enhance computational efficiency, non-homogeneous Dirichlet boundary conditions (DBCs) on the boundary surface are implemented. By formulating a variational problem related to the three-dimensional gravitational potential, we derive the associated finite-element equation. Then, the resulting linear system of equations from the finite-element algorithm is solved using a BiCGStab-Jacobi iterative solver. A model with depth-dependent residual density contrasts is utilized to evaluate the computational accuracy and performance of the proposed algorithm. Additionally, we examine a real sedimentary basin to illustrate the application of the finite-element forward modeling algorithm. All numerical results indicate that our forward algorithm can produce gravity responses for continuous density models efficiently and accurately.</p>