The numerical simulation of contact mechanics problems is computationally challenging, as these problems are locally highly non-linear and non-regular. Numerical solutions of such problems usually rely on adaptive mesh refinement (AMR). Even if efficient parallelizations of standard AMR techniques as h-adaptive methods begin to appear, their combination with contact problems remains a challenging task. Current developments on algorithms for contact problems are focusing either on non-parallelized new AMR methods or on parallelization methods for uniform refinement meshes. This work introduces a High Performance Computing strategy for solving 3D contact elastostatics problems with AMR on hexahedral elements. The contact is treated by a node-to-node algorithm with a penalization technique to deal with primal variables only. It presents the advantages of well modelling the studied phenomenon while not increasing the number of unknowns and not modifying the formulation in an intrusive manner. Concerning the AMR strategy, we rely on a non-conforming h-adaptive refinement solution. This method has shown to be well scalable. Regarding the detection of the refinement zones, a Zienkiewicz-Zhu (ZZ) type error estimator is used to select the elements to be refined through a local detection criterion. A stopping criterion is applied to automatically stop the refinement process. This combined strategy has proven its efficiency. In this contribution, we extend the combination of these contact mechanics and AMR strategies to a parallel framework. To carry out simulations, we place ourselves in the MFEM software environment. The proposed scalable algorithm is first based on a mesh elements partitioning that guarantees the contact paired nodes to be on the same processors. The contact matrix is locally built. The combined AMR-contact algorithm is ruled by two nested iterative loops. The external loop concerns the AMR process while the internal one deals with the contact solution. The penalized contact problem is solved thanks to a penalty constrained solver. The contact solution process is performed until the active contact nodes do not vary. Once this loop converged, the AMR strategy is locally applied and the mesh decomposition is rebalanced with the previously discussed partitioning. The process ends once the AMR stopping criterion is satisfied. To perform the AMR on MFEM, the implemented h-adaptive method is enriched with our own detection and stopping criteria.