Capturing the behavior of faults and fractures is critical to accurately simulate the geomechanical response of a complex subsurface system, such as an aquifer or a hydrocarbon reservoir. Several phenomena, such as micro-seismicity and fracture propagation, depend on the nature of these discontinuities. As a result, it is important to explicitly account for the complex behavior across these fractures and, from a modeling viewpoint, to simulate the influence of the fractures on the mechanical deformation. Dealing with frictional contact problems lies at the core of the challenge [1]. This is one of the most challenging problems in computational mechanics, since it usually produces a stiff non-linear problem associated with a series of linear systems, that is hard to solve efficiently. We use Lagrange multipliers to enforce the constraint and we focus on two different discretization techniques, only one of which is intrinsically stable. To solve the saddle-point Jacobian matrices arising from the linearization, we propose a constraint preconditioner where the primal Schur complement is obtained by eliminating the Lagrange multipliers unknowns. Suitable augmentation is presented for the intrinsically stable case. Finally, an optimal multigrid method [2] is applied to efficiently solve the Schur complement. We provide numerical evidences of the robustness and efficiency by solving large size problems from various applications. Funding was partially provided by TotalEnergies through the FC-MAELSTROM project. REFERENCES [1] Wriggers, P. Computational Contact Mechanics, Springer Berlin Heidelberg, 2nd Ed., (2006). [2] Isotton, G., and Frigo, M. and Spiezia, N. and Janna, C. Chronos: A general purpose classical AMG solver for High Performance Computing, SIAM J. Sci. Comput. (2021) 43(5):C335-C357.