Formation anisotropy in complicated geophysical environments can have a significant impact on data interpretation of electromagnetic surveys. To facilitate full 3D modeling of arbitrary anisotropy, we have adopted an h-version geometric multigrid preconditioned finite-element method based on vector basis functions. By using a structured mesh, instead of an unstructured one, our method can conveniently construct the restriction and prolongation operators for multigrid implementation, and then recursively coarsen the grid with the F-cycle coarsening scheme. The geometric multigrid method is used as a preconditioner for the biconjugate-gradient stabilized method to efficiently solve the linear system resulting from the finite-element method. Our method avoids the need of interpolation for arbitrary anisotropy modeling as in Yee’s grid-based finite-difference method, and it is also more capable of large-scale modeling with respect to the p-version geometric multigrid preconditioned finite-element method. A numerical example in geophysical well logging is included to demonstrate its numerical performance. Our h-version geometric multigrid preconditioned finite-element method is expected to help formation anisotropy characterization with electromagnetic surveys in complicated geophysical environments.