As a commonly used method for near-surface surveys, there is a demand for lightweight, fast, and accurate 3D DC geoelectric resistivity modeling that can be execute on portable computers. This study employed a singularity removal technique, solving the primary potential equation analytically and the secondary potential equation numerically using the finite element method. The modeling domain was discretized using orthosceme elements derived from dividing a hexahedron into six elements. The modeling domain (area of interest [AOI] + padding) covered [-300, 300] m × [-250, 250] m × [0, 250] m, divided into 59 × 39 × 20 nodes. The global matrix equation of the secondary potential is solved using the preconditioned conjugate gradients (PCG) with incomplete Cholesky factorization (ICF) as a preconditioner, making it 6.7 times faster than the direct solver. The combination of the singularity removal, orthosceme discretization, and the PCG solver enabled lightweight, fast, and accurate 3D resistivity modeling on a portable computer. Benchmarking this modeling against a layered Earth model showed that the Wenner and Schlumberger arrays achieved apparent resistivity with a relative error of mostly <
5 %. Applying this modeling to the vertical contact and buried block 3D anomaly Earth models demonstrated that the apparent resistivity profiles from Wenner, Schlumberger, and dipole-dipole arrays effectively captured resistivity changes at the anomaly boundaries, aligning with previous studies.