Solutions of variational inequalities often have limited regularity. In particular, the nonsmooth parts are local, while other parts of the solution have higher regularity. To overcome this limitation, we apply hp-adaptivity, which uses a combination of locally finer meshes and varying polynomial degrees to separate the different features of the the solution. For this, we employ Discontinuous Galerkin (DG) methods and show some novel error estimates for the obstacle problem which emphasize the use in hp-adaptive algorithms. Besides this analysis, we present how to efficiently compute numerical solutions using error estimators, fast algebraic solvers which can also be employed in a parallel setup, and discuss implementation details. Finally, some numerical examples and applications to phase field models are presented.
