  In this paper, we couple regularization techniques with the adaptive
$hp$-version of the boundary element method ($hp$-BEM) for the efficient
numerical solution of linear elastic problems with nonmonotone contact boundary
conditions. As a model example we treat the delamination of composite
structures with a contaminated interface layer. This problem has a weak
formulation in terms of a nonsmooth variational inequality. The resulting
hemivariational inequality (HVI) is first regularized and then, discretized by
an adaptive $hp$-BEM. We give conditions for the uniqueness of the solution and
provide an a-priori error estimate. Furthermore, we derive an a-posteriori
error estimate for the nonsmooth variational problem based on a novel
regularized mixed formulation, thus enabling $hp$-adaptivity. Various numerical
experiments illustrate the behavior, strengths and weaknesses of the proposed
high-order approximation scheme.
