Modern geoelectrical data acquisition systems allow large amounts of data to be collected in a short time. Inversions of such data sets require powerful forward solvers for predicting the electrical potentials. State-of-the-art solvers are typically based on finite elements. Recent developments in numerical mathematics led to direct matrix solvers that allow the equation systems arising from such finite element problems to be solved very efficiently. They are particularly useful for 3D geoelectrical problems, where many electrodes are involved. Although modern direct matrix solvers include optimized memory saving strategies, their application to realistic, large-scale 3D problems is still somewhat limited. Therefore, we present two novel techniques that allow the number of grid points to be reduced considerably, while maintaining a high solution accuracy. In the areas surrounding an electrode array we attach infinite elements that continue the electrical potentials to infinity. This does not only reduces the number of grid points, but also avoid the artificial Dirichlet or mixed boundary conditions that are well known to be the cause of numerical inaccuracies. Our second development concerns the singularity removal in the presence of significant surface topography. We employ a fast multipole boundary element method for computing the singular potentials. This renders unnecessary mesh refinements near the electrodes, which results in substantial savings of grid points of up to more than 50%. By means of extensive numerical tests we demonstrate that combined application of infinite elements and singularity removal allows the number of grid points to be reduced by a factor of $\approx$ 6 -- 10 compared with traditional finite element methods. This will be key for applying finite elements and direct matrix solver techniques to realistic 3D inversion problems.