We describe two different algorithms for automatic generation (without any user interaction) of optimal goal-oriented grids within a 3D higher-order Finite Element (FE) package, in context of simulation of electromagnetic borehole resistivity measurements for the assessment of rock formation properties. The first algorithm is based on transferring a two-dimensional optimal grid obtained from a 2D self-adaptive method. The resulting grid enables fast simulations (few seconds per logging position), although the accuracy of the final 3D results may be compromised in the case of highly deviated wells when combined with rapid spatial variations in electrical conductivity. For those cases, we developed a second algorithm based on a fully automatic three-dimensional goal-oriented strategy that, although it is considerably slower, it guarantees accurate results. A proper combination of both algorithms enables highly accurate numerical simulations of a variety of 3D resistivity logging measurements at different frequencies. This is illustrated in the manuscript with a Through Casing Resistivity (TCR) model problem. Numerical results corresponding to laterolog measurements indicate that differential voltage measurements remain sensitive to the 3D effects of various 3D sources and receivers only in deviated wells. In vertical wells, measurements are insensitive to the 3D effects of the sources and receivers. Thus, the resulting 3D combined effect (due to the use of 3D sources and receivers and deviation of the well) is different from the superposition of individual 3D effects.