Elastoplastic contact problems with hardening are ubiquitous in industrial metal
forming processes as well as many other areas. From a mathematical perspective, they are characterized
by the diculties of variational inequalities for both the plastic behavior as well as the
contact problem. Computationally, they also often lead to very large problems. In this paper, we
present and evaluate a set of methods that allows us to eciently solve such problems. In particular,
we use adaptive nite element meshes with linear and quadratic elements, a Newton linearization
of the plasticity, active set methods for the contact problem, and multigrid-preconditioned linear
solvers. Through a sequence of numerical experiments, we show the performance of these methods.
This includes highly accurate solutions of a benchmark problem and scaling our methods to 1,024
cores and more than a billion unknowns.