Contact mechanics can be addressed numerically by finite elements using either a penalty formulation or the Lagrange multipliers. The penalty approach leads to a linearized symmetric positive definite system which can prove severely ill-conditioned, with the iterative solution to large 3D problems requiring expensive preconditioners to accelerate, or even to allow for, convergence. In the present paper a mixed constraint preconditioner (MCP) is developed and implemented for the iterative solution of contact problems based on the penalty approach. A theoretical bound for the eigenspectrum of the preconditioned matrix is provided depending on the user-specified parameters that control the preconditioner fill-in. The performance of MCP combined with the preconditioned conjugate gradient is studied in two real 3D test cases describing the geomechanics of faulted porous media and compared to a well-established preconditioner such as the incomplete Cholesky decomposition ILLT with optimal fill-in. Numerical results show that typically MCP outperforms ILLT especially in the most challenging problems where very large penalty parameters must be used. Finally, a modification aimed at simplifying the MCP implementation and increasing its robustness is developed and experimented with, showing that the modified MCP can be viewed as an efficient and viable alternative to both MCP and ILLT, especially in problems where the size of the contact surfaces is large with respect to the global domain.