Predicting rolling bearing fatigue life requires knowledge of the three-dimensional (3D) stress fields in the roller and raceway near the lubricated contact. Owing to the increasingly severe operating conditions, the effect of localized features such as surface roughness, subsurface inclusions, and even the crystallographic structure of the material becomes important. Achieving such detail requires (locally) extremely dense gridding in simulations, which in 3D is a major challenge. Multigrid techniques have been demonstrated to be capable of solving such problems. In this study, multigrid techniques are shown to further increase the efficiency of the solution by exploiting local grid refinement while maintaining the simplicity of a uniform discretization. This is achieved by employing increasingly finer grids only locally, where the highest resolution is required. Results are presented for dry contact and elastohydrodynamically lubricated contact cases, circular as well as elliptic, with varying crystallographic structure, and with surface roughness. The results show that the developed algorithm is very well suited for detailed analysis, with also excellent prospects for computational diagnostics involving actual material crystallographic structure from electron backscatter diffraction measurements.