We present a 2.5-dimensional (2.5-D) finite element algorithm for direct current (DC) resistivity modelling in anisotropic media with singularity removal. First, we provide the weak form of the integral equation for the boundary value problem and simplify the Euler angles while calculating the primary potential so that the Fourier transform of the background potential with the dip angle can be avoided because it is mathematically difficult. A two-layered model is then simulated when the first covering is anisotropic. The relative error between this numerical solution and the analytical solution is < 1%. We then model a number of more complicated scenarios, using the algorithm developed in this paper. We test the model response to a small body at depth whose resistivity is isotropic, and then test whether the longitudinal or transverse resistivities affect the final results more. Based on this analysis, we found that longitudinal resistivity has more of an effect on the apparent resistivity than transverse resistivity in collinear arrays, such as pole-pole, dipole-dipole and Wenner arrays. Finally, through calculation of the current density and anomalous current density of several arrays, we conclude that the causes of different responses of longitudinal and transverse resistivity by each array is the distribution of current density in the subsurface. We also show that the sensitivity of each array type to variations in longitudinal and transverse resistivity can be understood when looked at from the perspective of current density.