On the coupling of regularization techniques and the boundary element method for a hemivariational inequality modelling a delamination problem
NNoname manuscript No. (will be inserted by the editor)
Coupling regularization and adaptive hp -BEM forthe solution of a delamination problem Nina Ovcharova · Lothar Banz
Received: date / Accepted: date
Abstract
In this paper, we couple regularization techniques with the adaptive hp -version of the boundary element method ( hp -BEM) for the efficient numer-ical solution of linear elastic problems with nonmonotone contact boundaryconditions. As a model example we treat the delamination of composite struc-tures with a contaminated interface layer. This problem has a weak formulationin terms of a nonsmooth variational inequality. The resulting hemivariationalinequality (HVI) is first regularized and then, discretized by an adaptive hp -BEM. We give conditions for the uniqueness of the solution and provide ana-priori error estimate. Furthermore, we derive an a-posteriori error estimatefor the nonsmooth variational problem based on a novel regularized mixedformulation, thus enabling hp -adaptivity. Various numerical experiments il-lustrate the behavior, strengths and weaknesses of the proposed high-orderapproximation scheme. Keywords
Hemivariational inequality · Regularization techniques · Delami-nation problem · a-priori and a-posteriori error estimates · hp -BEM Mathematics Subject Classification (2000) · The interest in robust and efficient numerical methods for the solution of nons-mooth problems in nonmonoton contact like adhesion and delamination prob-lems increases constantly. The nonsmoothness comes from the nonsmooth data
N. OvcharovaUniversit¨at der Bundeswehr M¨unchen, D-85577 Neubiberg/Munich, GermanyE-mail: [email protected]. BanzDepartment of Mathematics, University of Salzburg, Hellbrunner Straße 34, 5020 Salzburg,AustriaE-mail: [email protected] a r X i v : . [ m a t h . NA ] J un Nina Ovcharova, Lothar Banz of the problems itself, in particular from nonmonotone multivalued physicallaws involved in the boundary conditions that lead to nonsmooth functionalsin the variational formulation. There are several approaches to treat this non-differentiability. We can first discretize by finite elements, and then solve thenonconvex optimization problems by novel nonsmooth optimization methods[16], or first regularize the nonsmooth functional, then discretize by finite el-ement methods and finally use standard optimization solvers [32]. Since theinteresting nonmonotone behavior takes place only in a boundary part of theinvolved linear elastic bodies, the domain hemivariational inequaliy can also beformulated as a boundary hemivariational inequality via the Poincar´e-Steklovoperator. Therefore, the boundary element methods (BEM) are the methodsof choice for discretization. There are number of works exploiting the h - andthe high-order hp -version of the BEM for monotone contact problems. Ad-vanced hp -BEM discretizations, a-priori and a-posteriori error estimates forunilateral Signorini problems or contact problems with monotone friction areestablished in [3,4,12,29,30,37]. However, successfully applied hp -BEM tech-niques for nonsmooth nonconvex variational inequalities are still missing. Pure h -BEM with lowest polynomial degree has been investigated for the first timeby Nesemann and Stephan [31] for unilateral contact with adhesion. Theystudy the existence and uniqueness, and introduce a residual error indicator.In this paper, we extend the approximation scheme from [33] (dealing withthe h - version of the BEM) to the hp -version of the BEM for the regularizedproblem. More precisely, the nonsmooth variational inequality is first regu-larized and then discretized by hp -adaptive BEM. For the Galerkin solutionof the regularized problem we provide an a-priori error estimate and derive areliable a-posteriori error estimate based on a novel regularized mixed formula-tion and thus, enabling hp -adaptivity. All our theoretical results are illustratedwith various numerical experiments for a contact problem with adhesion. Wealso provide conditions for the uniqueness of the solution that sharpen theresults due to [31]. Notation:
We denote with C i , c i and such alike generic constants, whichcan take different values at different positions. We use bold font to indicatevector-valued variables, e.g. u ∈ H ( Ω ) := [ H ( Ω )] d . If there are too manyindices, one of them is written as a superscript, e.g. u ε,hp = u εhp . Let Ω ⊂ R d ( d = 2 ,
3) be a bounded domain with Lipschitz boundary Γ = ∂Ω .We assume that the boundary is decomposed into three disjoint open parts Γ D , Γ N , and Γ C such that Γ = Γ D ∪ Γ N ∪ Γ C and, moreover, the measuresof Γ C and Γ D are positive. We fix an elastic body occupying Ω . The bodyis subject to volume force f ∈ [ L ( Ω )] d and g ∈ H / ( Γ C ), g ≥
0, is a gapfunction associating every point x ∈ Γ C with its distance to the rigid obstaclemeasured in the direction of the unit outer normal vector n ( x ). The body is p -adaptive BEM for regularized hemivariational inequalities 3(a) Load-displacement curve for different con-tamination concentrations, see [38] g - u n σ n (b) A nonmonotone adhesion law Fig. 1
Nonmonotone adhesion law fixed along Γ D , surface tractions t ∈ [ L ( Γ N )] d act on Γ N , and on the part Γ C a nonmonotone, generally multivalued boundary condition holds.Further, ε ( u ) = ( ∇ u + ∇ u T ) denotes the linearized strain tensor and σ ( u ) = C : ε ( u ) stands for the stress tensor, where C is the Hooke tensor,assumed to be uniformly positive definite with L ∞ coefficients. The boundarystress vector can be decomposed further into the normal, respectively, thetangential stress: σ n = σ ( u ) n · n , σ t = σ ( u ) n − σ n n . By assuming that the structure is symmetric and the forces applied to theupper and lower part are the same, it suffices to consider only the upper half ofthe specimen represented by Ω , see Fig. 2 left for the 2D benchmark problem. t D Fig. 2
Reference configuration for the 2D benchmark under loading. Under applied tractionforce t the crack front propagates to the left Nina Ovcharova, Lothar Banz The delamination problem under consideration is the following: Find u ∈ H ( Ω ) := [ H ( Ω )] d such that − div σ ( u ) = f in Ω (1a) u = 0 on Γ D (1b) σ ( u ) n = t on Γ N (1c) u n ≤ g on Γ C (1d) σ t ( u ) = 0 on Γ C (1e) − σ n ( u ) ∈ ∂f ( u n ) on Γ C . (1f)The contact law (1f), written as a differential inclusion by means of the Clarkesubdifferential ∂f [13] of a locally Lipschitz function f , describes the nonmono-tone, multivalued behavior of the adhesive. More precisely, ∂f is the physicallaw between the normal component σ n of the stress boundary vector and thenormal component u n = u · n of the displacement u on Γ C . A typical zig-zagged nonmonotone adhesion law is shown in Figure 1(b).To give a variational formulation of the above boundary value problem wedefine H D ( Ω ) = (cid:8) v ∈ H ( Ω ) : v | Γ D = 0 (cid:9) , K = { v ∈ H D ( Ω ) : v n ≤ g a.e. on Γ C } and introduce the H D ( Ω )-coercive and continuous bilinear form of linear elas-ticity a ( u , v ) = (cid:90) Ω σ ( u ) : ε ( v ) dx. Multiplying the equilibrium equation (1a) by v − u , integrating over Ω andapplying the divergence theorem yields (cid:90) Ω σ ( u ) : ε ( v − u ) dx = (cid:90) Ω f · ( v − u ) dx + (cid:90) Γ σ ( u ) n · ( v − u ) ds. From the definition of the Clarke subdifferential, the nonmonotone boundarycondition (1f) is equivalent to − σ n ( u )( v n − u n ) ≤ f ( u n ; v n − u n ) . Here, the notation f ( x ; z ) stands for the generalized directional derivativeof f at x in direction z . Substituting σ ( u ) n by t on Γ N , using on Γ C thedecomposition σ ( u ) n · ( v − u ) = σ t ( u ) · ( v t − u t ) + σ n ( u )( v n − u n )and taking into account that on Γ C no tangential stresses are assumed, c.f. (1e),we obtain the hemivariational inequality (HVI): Find u ∈ K such that a ( u , v − u ) + (cid:90) Γ C f ( u n ( s ); v n ( s ) − u n ( s )) ds ≥ (cid:90) Ω f · ( v − u ) dx + (cid:90) Γ N t · ( v − u ) ds ∀ v ∈ K . (2) p -adaptive BEM for regularized hemivariational inequalities 5 Since the main difficulties of the boundary value problem, namely the non-monotone adhesion law (1f) and the unilateral contact condition (1d) ap-pear on the boundary, it might be viewed advantageous to formulate (2) as aboundary integral problem. To this end, we introduce the free boundary part Γ = Γ \ Γ D = Γ N ∪ Γ C and recall the Sobolev spaces [25]: H / ( Γ ) = { v ∈ L ( Γ ) : ∃ v (cid:48) ∈ H ( Ω ) , tr v (cid:48) = v } ,H / ( Γ ) = { v = v (cid:48) | Γ : ∃ v (cid:48) ∈ H / ( Γ ) } , ˜ H / ( Γ ) = { v = v (cid:48) | Γ : ∃ v (cid:48) ∈ H / ( Γ ) , supp v (cid:48) ⊂ Γ } with the standard norms (cid:107) u (cid:107) H / ( Γ ) = inf v ∈ H / ( Γ ) ,v | Γ = u (cid:107) v (cid:107) H / ( Γ ) and (cid:107) u (cid:107) ˜ H / ( Γ ) = (cid:107) u (cid:107) H / ( Γ ) , where u is the extension of u onto Γ by zero. The Sobolev space of negativeorder on Γ are defined by duality as H − / ( Γ ) = ( ˜ H / ( Γ )) ∗ and ˜ H − / ( Γ ) = ( H / ( Γ )) ∗ . Moreover, from [25, Lemma 4.3.1] we have the inclusions˜ H / ( Γ ) ⊂ H / ( Γ ) ⊂ L ( Γ ) ⊂ ˜ H − / ( Γ ) ⊂ H − / ( Γ ) . For the solution u ( x ) of (1a) with x ∈ Ω \ Γ we have the following repre-sentation formula, also known as Somigliana’s identity, see e.g. [26]: u ( x ) = (cid:90) Γ G ( x , y ) ( T y u ( y )) ds y − (cid:90) Γ T y G ( x , y ) u ( y ) ds y + (cid:90) Ω G ( x , y ) f ( y ) d y , (3) where G ( x , y ) is the fundamental solution of the Navier-Lam´e equation de-fined by G ( x , y ) = λ +3 µ πµ ( λ +2 µ ) (cid:16) log | x − y | I + λ + µλ +3 µ ( x − y )( x − y ) (cid:62) | x − y | (cid:17) , if d=2 λ +3 µ πµ ( λ +2 µ ) (cid:16) | x − y | − I + λ + µλ +3 µ ( x − y )( x − y ) (cid:62) | x − y | (cid:17) , if d=3with the Lam´e constants λ, µ > E and the Poisson’s ratio ν : λ = Eν − ν , µ = E ν . Here, T y stands for the traction operator with respect to y defined by T y ( u ) := σ ( u ( y )) · n y , where n y is the unit outer normal vector at y ∈ Γ . Letting Ω \ Γ (cid:51) x → Γ in (3), we obtain the well-known Calder´on operator (cid:32) uT x u (cid:33) = (cid:32) I − K VW I + K (cid:48) (cid:33) (cid:32) uT x u (cid:33) + (cid:32) N f N f (cid:33) , Nina Ovcharova, Lothar Banz with the single layer potential V , the double layer potential K , its formaladjoint K (cid:48) , and the hypersingular integral operator W defined for x ∈ Γ asfollows:( V φ ) ( x ) := (cid:90) Γ G ( x , y ) φ ( y ) ds y , ( Kφ ) ( x ) := (cid:90) Γ T y G (cid:62) ( x , y ) φ ( y ) ds y ( K (cid:48) φ ) ( x ) := T x (cid:90) Γ G ( x , y ) φ ( y ) ds y , ( W φ ) ( x ) := − T x ( Kφ ) ( x ) . The Newton potentials N , N are given for x ∈ Γ by N f = (cid:90) Γ G ( x , y ) f ( y ) ds y , N f = T x (cid:90) Γ G ( x , y ) f ( y ) ds y . From [14] it is known that the linear operators V : H − / σ ( Γ ) → H / σ ( Γ ) , K : H / σ ( Γ ) → H / σ ( Γ ) K (cid:48) : H − / σ ( Γ ) → H − / σ ( Γ ) , W : H / σ ( Γ ) → H − / σ ( Γ )are well-defined and continuous for | σ | ≤
12 . Moreover, V is symmetric andpositive definite (elliptic on H − / ( Γ )) in R and, if the capacity of Γ issmaller than 1, also in R . That can always be achieved by scaling, since thecapacity (or conformal radius or transfinite diameter) of Γ is smaller than 1,if Ω is contained in a disc with radius < W is symmetricand positive semidefinite with kernel R (elliptic on ˜ H / ( Γ )). Hence, since V : H − / ( Γ ) → H / ( Γ ) is invertible, we obtain by taking the Schur complementof the Calderon projector that T x u = P u − N f (4)with the symmetric Poincar´e-Steklov operator P : H / ( Γ ) → H − / ( Γ ) andits Newton potential N given by P u = W u + (cid:18) K (cid:48) + 12 I (cid:19) V − (cid:18) K + 12 I (cid:19) u , N f = (cid:18) K (cid:48) + 12 I (cid:19) V − N f − N f . Note that if f = 0, P maps u to its traction and, therefore, the Poincar´e-Steklov operator P is sometimes called the Dirichlet-to-Neumann mapping.Moreover, P induces a symmetric bilinear form on H / ( Γ ), and is continuousand ˜ H / ( Γ )-elliptic, see e.g. [9], i.e. there exist constants c P , C P > (cid:107) P v (cid:107) H − / ( Γ ) ≤ C P (cid:107) v (cid:107) H / ( Γ ) ∀ v ∈ H / ( Γ ) , (cid:104) P v , v (cid:105) Γ ≥ c P (cid:107) v (cid:107) ˜H / ( Γ ) ∀ v ∈ ˜H / ( Γ ) . Here, (cid:104)· , ·(cid:105) Γ is the L ( Γ )-duality pairing between the involved spaces. Usingthe boundary functional sets V = ˜ H / ( Γ ) and K Γ = { v ∈ V : v n ≤ g a.e. on Γ C } p -adaptive BEM for regularized hemivariational inequalities 7 we obtain as in the domain based cases the boundary integral hemivariationalinequality (BIHVI), Problem ( P ): Find u ∈ K Γ such that (cid:104) P u , v − u (cid:105) Γ + (cid:90) Γ C f ( u n ( s ); v n ( s ) − u n ( s )) ds ≥ (cid:104) t , v − u (cid:105) Γ N + (cid:104) N f , v − u (cid:105) Γ ∀ v ∈ K Γ . (5)To shorten the right hand side we introduce the linear functional (cid:104) F , v (cid:105) Γ = (cid:90) Γ N t · v ds + (cid:104) N f , v (cid:105) Γ . In this section, we review from [32,34] a class of smoothing approximations fornonsmooth functions that can be re-expressed by means of the plus functionp( t ) = t + = max { t, } . The approximation is based on smoothing of the plusfunction via convolution.Let ε > f : R × R ++ → R of a locally Lipschitz function f : R → R is defined viaconvolution by ˆ f ( x, ε ) = (cid:90) R f ( x − εt ) ρ ( t ) dt, where ρ : R → R + is a probability density function such that κ = (cid:90) R | t | ρ ( t ) dt < ∞ and R + = { ε ∈ R : ε ≥ } , R ++ = { ε ∈ R : ε > } . We consider a class of nonsmooth functions that can be reformulated by meansof the plus function. To this class of functions belong maximum, minimum ornested max-min function. If f is a maximum function, for example, f ( x ) =max { g ( x ) , g ( x ) } , then f can be reformulated by using the plus function as f ( x ) = g ( x ) + p( g ( x ) − g ( x )) , (6)The smoothing function of f is then defined by S ( x, ε ) = g ( x ) + ˆ p ( g ( x ) − g ( x ) , ε ) , (7)where ˆ p ( t, ε ) is the smoothing function of p( t ) via convolution.Choose, for example, the Zang probability density function ρ ( t ) = (cid:40) − ≤ t ≤ . Nina Ovcharova, Lothar Banz
Then ˆp( t, ε ) = (cid:90) R p( t − εs ) ρ ( s ) ds = t < − ε ε ( t + ε ) if − ε ≤ t ≤ ε t if t > ε and hence, S ( x, ε ) = g ( x ) if g ( x ) − g ( x ) ≤ − ε ε [ g ( x ) − g ( x )] + ( g ( x ) + g ( x )) + ε if | g ( x ) − g ( x ) | ≤ ε g ( x ) if g ( x ) − g ( x ) ≥ ε . (8)The relation (7) can be extended to the maximum function f : R → R of m continuous functions g , . . . , g m , i.e. f ( x ) = max { g ( x ) , g ( x ) , . . . , g m ( x ) } (9)by S ( x, ε ) = g ( x ) + ˆp ( g ( x ) − g ( x ) + . . . + ˆp ( g m ( x ) − g m − ( x ) , ε ) , ε ) . (10)The major properties of the function S in (10) are listed in the followinglemma: Lemma 1 [32, lemma 10] Let f and S be defined as in (9) , (10) , respectively,with g i ∈ C . Then there holds:(i) For any ε > and for all x ∈ R , | S ( x, ε ) − f ( x ) | ≤ ( m − κε. (ii) The function S is continuously differentiable on R × R ++ and for any x ∈ R and ε > there exist Λ i ∈ [0 , such that m (cid:88) i =1 Λ i = 1 and ∂S ( x, ε ) ∂x = S x ( x, ε ) = m (cid:88) i =1 Λ i g (cid:48) i ( x ) . (11) Moreover, (cid:26) lim z → x,ε → + S x ( z, ε ) (cid:27) ⊆ ∂f ( x ) . (12) Remark 1
Since f ( x ) = min { g ( x ) , . . . , g m ( x ) } = − max {− g ( x ) , . . . , − g m ( x ) } these type of functions can be handled as above. Assumption 1
Assume that there exists positive constants c i , d i such that forall x ∈ R | g (cid:48) i ( x ) | ≤ c i (1 + | x | ) (13a) g (cid:48) i ( x ) · ( − x ) ≤ d i | x | . (13b) p -adaptive BEM for regularized hemivariational inequalities 9 By using (11) - (12) and (13a) - (13b), the following auxiliary result can beeasily deduced.
Lemma 2
For any x, z, ξ ∈ R it holds that | S x ( x, ε ) · z | ≤ c (1 + | x | ) | z | (14a) S x ( x, ε ) · ( − x ) ≤ d | x | (14b)lim sup z → x,ε → + S x ( ε, x ) ξ ≤ f ( x ; ξ ) . (14c)Further, we introduce the functional J ε : H / ( Γ ) → R defined by J ε ( u ) = (cid:90) Γ C S ( u n ( s ) , ε ) ds. The regularized problem ( P ε ) of ( P ) is given by: Find u ε ∈ K Γ such that (cid:104) P u ε , v − u ε (cid:105) Γ + (cid:104) DJ ε ( u ε ) , v − u ε (cid:105) Γ C ≥ (cid:104) F , v − u ε (cid:105) Γ ∀ v ∈ K Γ , (15)where DJ ε : H / ( Γ ) → H / ( Γ ) is the Gˆateaux derivative of the functional J ε and is given by (cid:104) DJ ε ( u ) , v (cid:105) Γ C = (cid:90) Γ C S x ( u n ( s ) , ε ) v n ( s ) ds. To simplify the notations we introduce ϕ : H / ( Γ ) × H / ( Γ ) → R definedby ϕ ( u , v ) = (cid:90) Γ C f ( u n ( s ); v n ( s ) − u n ( s )) ds ∀ u , v ∈ H / ( Γ ) . (16)The regularized version of ϕ ( · , · ) is therewith ϕ ε : H / ( Γ ) × H / ( Γ ) → R , ϕ ε ( u , v ) := (cid:104) DJ ε ( u ) , v − u (cid:105) Γ C . The existence result for (5), resp. (15), is due to [19] and relies in both caseson the pseudomonotonicity of ϕ and ϕ ε , respectively. For details we refer thereader to [32,34].Finally, we recall that the functional ϕ : X × X → R , where X is a realreflexive Banach space, is pseudomonotone if u n (cid:42) u (weakly) in X andlim inf n →∞ ϕ ( u n , u ) ≥ v ∈ X , we have lim sup n →∞ ϕ ( u n , v ) ≤ ϕ ( u, v ) . In this section, we discuss the uniqueness of solution of the boundary hemivari-ational inequality and of the corresponding regularization problem. The mainresults are presented in Theorem 2 which is based on the abstract uniquenessTheorem 1 from [33] that gives a sufficient condition for uniqueness. Similaruniqueness result but for the regularized problem is derived in Theorem 3.We are also dealing with the question which classes of smoothing functionspreserve the property of unique solvability of the original problem as ε → + . Assumption 2
Assume that there exists an α ∈ [0 , c P ) such that for any u , v ∈ V it holds ϕ ( u , v ) + ϕ ( v , u ) ≤ α (cid:107) u − v (cid:107) V . (17) To make the results of uniqueness self-consistent, we introduce from [33] thefollowing theorem.
Theorem 1 [33, Theorem 5.1] Under the assumption (17) with α < c P ,there exists a unique solution to the BIHVI problem ( P ) , which depends Lips-chitz continuously on the right hand side F ∈ V ∗ . Proof
Assume that u , ˜ u are two solutions of ( P ). Then the inequalities belowhold: (cid:104) P u − F , v − u (cid:105) Γ + ϕ ( u , v ) ≥ ∀ v ∈ K Γ (cid:104) P ˜ u − F , v − ˜ u (cid:105) Γ + ϕ (˜ u , v ) ≥ ∀ v ∈ K Γ . Setting v = ˜ u in the first inequality and v = u in the second one, and summingup the resulting inequalities, we get (cid:104) P u − P ˜ u , ˜ u − u (cid:105) Γ + ϕ ( u , ˜ u ) + ϕ (˜ u , u ) ≥ . (18)From the coercivity of the operator P and the assumption (17) we obtain c P (cid:107) u − ˜ u (cid:107) V ≤ ϕ ( u , ˜ u ) + ϕ (˜ u , u ) ≤ α (cid:107) u − ˜ u (cid:107) V . Hence, since α ∈ [0 , c P ), if u (cid:54) = ˜ u we receive a contradiction.Now let F i ∈ V ∗ and denote u i = u F i , i = 1 , . Analogously to (18), we findthat (cid:104) P u − F − P u + F , u − u (cid:105) Γ + ϕ ( u , u ) + ϕ ( u , u ) ≥ . Hence, c P (cid:107) u − u (cid:107) V ≤ ϕ ( u , u ) + ϕ ( u , u ) + (cid:104) F − F , u − u (cid:105) Γ and by (17),( c P − α ) (cid:107) u − u (cid:107) V ≤ (cid:104) F − F , u − u (cid:105) Γ ≤ (cid:107) F − F (cid:107) V ∗ (cid:107) u − u (cid:107) V . Also, since α < c P we deduce that (cid:107) u − u (cid:107) V ≤ c P − α (cid:107) F − F (cid:107) V ∗ , which concludes the proof of the theorem. (cid:117)(cid:116) Further, we present a class of locally Lipschitz functions for which thecrucial assumption (17) is satisfied. Let f : R → R be a function such that( ξ ∗ − η ∗ ) ( ξ − η ) ≥ − α | ξ − η | ∀ ξ ∗ ∈ ∂f ( ξ ) , ∀ η ∗ ∈ ∂f ( η ) (19) p -adaptive BEM for regularized hemivariational inequalities 11 f Fig. 3
An exemplary function b with one positive jump and its anti-derivative f ( x ) = (cid:82) x b ( t ) dt for any ξ, η ∈ R and some α ≥
0. From the definition of the Clarke generalizedderivative [13] we get f ( ξ ; η − ξ ) = max ξ ∗ ∈ ∂f ( ξ ) ξ ∗ ( η − ξ ) . Rewriting (19) as ξ ∗ ( η − ξ ) + η ∗ ( ξ − η ) ≤ α | ξ − η | we find f ( ξ ; η − ξ ) + f ( η ; ξ − η ) ≤ α | ξ − η | . (20)Hence, ϕ ( u , v ) + ϕ ( v , u ) = (cid:90) Γ C f ( u n ; v n − u n ) ds + (cid:90) Γ C f ( v n ; u n − v n ) ds ≤ α (cid:107) u n − v n (cid:107) L ( Γ C ) ≤ α (cid:107) u − v (cid:107) V (21)and consequently the assumption (17) is satisfied provided that α is suffi-ciently small ( α < c P ).Next, we show that if ∂f includes only non-negative jumps, then the con-dition (19) is globally satisfied. Whereas for negative jumps, (19) holds onlylocally. Example 1
Let f : R → R be a continuous, piecewise C , function such thatits first derivative b ( x ) := f (cid:48) ( x ) has finite non-negative jumps at the points x Ji of discontinuity (see Fig. 3). This means that b ( x Ji ) ≤ b ( x Ji ), where b ( x Ji ) := f (cid:48) ( x Ji −
0) = lim h → − f ( x Ji + h ) − f ( x Ji ) h and b ( x Ji ) := f (cid:48) ( x Ji + 0) = lim h → + f ( x Ji + h ) − f ( x Ji ) h . Hence, ∂f ( x Ji ) = (cid:2) b ( x Ji ) , b ( x Ji ) (cid:3) . Setting E = ∪ E i the finite set of open intervals E i , where the function b | E i issmooth with the Lipschitz constant c L i , we define α = max i c L i . (22)Let { x Ji } ki =1 be the set of jags between x and x with x < x J < · · · < x Jk
If the graph of ∂f consists of several decreasing straight line seg-ments and non-negative jumps, then the value − α in (19) is the steepestdecreasing slope. The next example treats the case of negative jumps.
Example 2
Let f : R → R be a continuous, piecewise C , function such thatits first derivative has finite negative jumps at the points x Ji of discontinuity(see Fig. 4). In this case, ∂f ( x Ji ) = (cid:2) b ( x Ji ) , b ( x Ji ) (cid:3) . Then, for any x > x , we have b ( x ) − b ( x ) x − x ≥ − α + k (cid:88) i =1 c Ji x − x , (23) p -adaptive BEM for regularized hemivariational inequalities 13 f Fig. 4
An exemplary function b with one negative jump and its anti-derivative f ( x ) = (cid:82) x b ( t ) dt with α ≥ c Ji = b ( x Ji ) − b ( x Ji ) < x Ji .Further, for the sake of simplicity, we consider a function with one negativejump c J and estimate as follows: c J ( x − x ) = c J ( x − x J ) + c J ( x J − x ) = c J ( x − x J ) + + c J ( x J − x ) + . (24)Using x ≥ ( x − a ) + for any a ≥ , we get( x − x ) ≥ ( x − x −
14 ) + ≥ ( x − x − ( x J − x )) + = ( x − x J ) + ∀ x ≤ x J − x − x ) ≥ ( x − x −
14 ) + ≥ ( x − x − ( x − x J )) + = ( x J − x ) + ∀ x ≥ x J + 14 . Hence, the condition (19) is fulfilled with α − c J >
0, but only for those x , x satisfying x ≥ x J + and x ≤ x J − .Let h ( t ) be a Lipschitz continuous function with Lipschitz constant c L . Weinclude the negative jump c J at the point t J through the Heaviside function Θ by b ( t ) = h ( t ) + c J Θ ( t − t J )and compute f ( t ; d ) = (cid:40) h ( t ) d, if t < t J h ( t ) d + c J d, if t > t J , f ( t J ; d ) = (cid:40) h ( t J ) d, if d > h ( t J ) d + c J d, if d ≤ . Hence, for any t , t < t J , we have f ( t ; t − t ) + f ( t ; t − t ) = h ( t )( t − t ) + h ( t )( t − t )=( h ( t ) − h ( t ))( t − t ) ≤ c L ( t − t ) . (25) Analogously, for any t , t > t J , we get f ( t ; t − t ) + f ( t ; t − t ) = h ( t )( t − t ) + c J ( t − t ) + h ( t )( t − t )+ c J ( t − t ) ≤ c L ( t − t ) . In both cases, the condition (20) is satisfied.Let now t < t J ≤ t , then we have f ( t ; t − t ) + f ( t ; t − t ) = h ( t )( t − t ) + h ( t )( t − t ) + c J ( t − t ) ≤ c L ( t − t ) + c J ( t − t ) ≤ α ( t − t ) (26)for some α if and only if c J t − t + c L ≤ α . (27)Thus, we assume c L < α and t ≤ t J − c J c L − α , (28)which, since t J ≤ t , implies (27) immediately.Based on the results presented in the Examples 1 and 2 we derive thefollowing uniqueness result. Theorem 2
Let u be a solution of the BIHVI problem ( P ) . Then u is uniqueif one of the following conditions holds:(a) The jumps c iJ are non-negative and the Lipschitz constant c L satisfies c L We note that our uniqueness result is more general then the resultin [31, Theorem 7], which demands ˜ c = 4 . p -adaptive BEM for regularized hemivariational inequalities 15 Assumption 3 Assume that the assumption (19) is satisfied for the regular-ization function, i.e. there exists a constant α ≥ (in general depending on ε > ) such that ( S x ( x , ε ) − S x ( x , ε )) ( x − x ) ≥ − α | x − x | ∀ x , x ∈ R . (30)Hence, for any u , v ∈ V , we have (cid:104) DJ ε ( u ) − DJ ε ( v ) , v − u (cid:105) Γ C = (cid:90) Γ C ( S x ( u n ( s ) , ε ) − S x ( v n ( s ) , ε )) ( v n ( s ) − u n ( s )) ds ≤ α (cid:107) u n − v n (cid:107) L ( Γ C ) ≤ α (cid:107) u − v (cid:107) V . (31)Due to Theorem 1, we have the following uniqueness result for the regularizedproblem. Theorem 3 Under the assumption (30) with α < c p , there exists a uniquesolution to the regularized problem ( P ε ) , which depends Lipschitz continuouslyon the right hand side F ∈ V ∗ . One of the most interesting question is, which classes of smoothing functionspreserve the property of unique solvability of the original problem as ε → + .As the next example shows, this is the case for the smoothing function (8)relevant to our computation. In particular, we show that the constant α in(30) does not depend on ε . Example 3 In view of our application to contact problems with nonmonotonelaws in the boundary conditions, we consider f ( x ) = max { g ( x ) , g ( x ) } with g i ( x ) = a i x + c i x + d i . We use the smoothing function (8) as a regularizationof f . The corresponding approximation of ∂f is given by S x ( x, ε ) = g (cid:48) ( x ) if g ( x ) − g ( x ) ≤ − ε λ g (cid:48) ( x ) + λ g (cid:48) ( x ) if | g ( x ) − g ( x ) | < ε g (cid:48) ( x ) if g ( x ) − g ( x ) ≥ ε , where λ = 1 ε ( g ( x ) − g ( x )) + 12 , λ = 1 − λ . We need to consider only the region I ε := { x ∈ R : | g ( x ) − g ( x ) | < ε } , where f is actually approximated. Otherwise, the condition (30) is automaticallysatisfied.Let x i ∈ I ε , i = 1 , 2. Setting g ( x ) := g ( x ) − g ( x ), we introduce λ i = 1 ε ( g ( x i ) − g ( x i )) = − ε g ( x i ) and λ i = 1 − λ i , i = 1 , . Hence, λ − λ = λ − λ = 1 ε ( g ( x ) − g ( x )) (32) and thus, S x ( x , ε ) − S x ( x , ε ) = λ ( a x + c ) + λ ( a x + c ) − λ ( a x + c ) − λ ( a x + c )= a λ ( x − x ) + a λ ( x − x ) + ( λ − λ )( a x + c )+ ( λ − λ )( a x + c )= (cid:0) a λ + a λ (cid:1) ( x − x ) + ( λ − λ )( g (cid:48) ( x ) − g (cid:48) ( x ))= (cid:0) a λ + a λ (cid:1) ( x − x ) + ( λ − λ ) g (cid:48) ( x ) . (33)Similarly to (33), we find that S x ( x , ε ) − S x ( x , ε ) = (cid:0) a λ + a λ (cid:1) ( x − x ) + ( λ − λ ) g (cid:48) ( x ) . (34)Summing now (33) and (34), we get2 ( S x ( x , ε ) − S x ( x , ε )) = ( a λ + a λ + a λ + a λ )( x − x )+ ( λ − λ )( g (cid:48) ( x ) + g (cid:48) ( x )) . (35)From (32), using the structure of g and g , we find that λ − λ = 1 ε ( g ( x ) − g ( x )) = 1 ε (cid:18) a − a x − x ) + ( c − c )( x − x ) (cid:19) = 1 ε ( x − x ) (cid:18) a − a x + c − c a − a x + c − c (cid:19) = 12 ε ( x − x )( g (cid:48) ( x ) + g (cid:48) ( x )) . (36)Hence,2 ( S x ( x , ε ) − S x ( x , ε )) ( x − x ) ≥ ( a λ + a λ + a λ + a λ )( x − x ) = ( a (1 − λ ) + a λ + a (1 − λ ) + a λ )( x − x ) ≥ ( −| a || (1 − λ ) | − | a || λ | − | a || (1 − λ ) | − | a || λ | )( x − x ) ≥ − max {| a | , | a |} ( | (1 − λ ) | + | λ | + | (1 − λ ) | + | λ | )( x − x ) ≥ − max {| a | , | a |} (3 / / / / x − x ) = − {| a | , | a |} ( x − x ) and consequently, (30) is satisfied with α = 2 max {| a | , | a |} . Remark 4 If we use (8) as a smoothing approximation for the maximum(resp. minimum) function of several quadratic functions a i x + c i x + d i , then(30) is satisfied and, therefore, the regularized problem is always uniquely solv-able if all | a i | are sufficiently small. Since the constant α does not depend on ε , the property of uniqueness is preserved as ε → + . p -adaptive BEM for regularized hemivariational inequalities 17 To avoid domain approximation, let Ω ⊂ R d , d = 2 , 3, also be a polygonaldomain. Let T h be a sufficiently fine finite element mesh of the boundary Γ respecting the decomposition of Γ into Γ D , Γ N and Γ C , p = ( p T ) T ∈T h apolynomial degree distribution over T h , P p T ( ˆ T ) the space of polynomials oforder p T on the reference element ˆ T , and Ψ T : ˆ T → T ∈ T h a bijective, (bi)-linear transformation. In 2D, ˆ T is the interval [ − , − , . Let Σ T,hp be the set of all affinely transformed (tensorproduct based) Gauss-Lobatto nodes on the element T of the partition T h of Γ , corresponding to the polynomial degree p T , and set Σ hp := (cid:91) T ∈T h | ΓC Σ T,hp ,see [27,29,20]. Furthermore, we assume in this section that g ∈ C ( Γ C ) toallow point evaluation.For the discretization of the displacement u we use V hp = { v hp ∈ C ( Γ ) : v hp | T ◦ Ψ T ∈ [ P p T ( ˆ T )] d ∀ T ∈ T h , v hp = 0 on Γ D } , K Γhp = { v hp ∈ V hp : ( v hp · n )( P i ) ≤ g ( P i ) ∀ P i ∈ Σ hp } . In general K Γhp (cid:54)⊆ K Γ . For the approximation of the Poincar´e-Steklov operator,namely V − , we need the space W hp = { ψ hp ∈ L ( Γ ) : ψ hp | T ◦ Ψ T ∈ [ P p T − ( ˆ T )] d ∀ T ∈ T h } ⊂ H − / ( Γ ) . For more details on the approximation techniques based on boundary elementmethod see e.g. [10,14,15,18,22,23,29,28,35]. Let { φ i } N D i =1 and { ψ j } N N j =1 be abases of V hp and W hp , respectively. Then, the boundary matrices are given by( V hp ) j,i = (cid:104) V ψ i , ψ j (cid:105) , ( K hp ) j,i = (cid:104) Kφ i , ψ j (cid:105) , ( W hp ) j,i = (cid:104) W φ i , φ j (cid:105) , ( I hp ) j,i = (cid:104) φ i , ψ j (cid:105) and, therewith, the approximation of the Galerkin matrix is P hp = W hp + (cid:18) K hp + 12 I hp (cid:19) (cid:62) V − hp (cid:18) K hp + 12 I hp (cid:19) . With the canonical embeddings i hp : W hp (cid:44) → H − / ( Γ ) , j hp : V hp (cid:44) → H / ( Γ )and their duals i ∗ hp and j ∗ hp , the discrete Poincar´e-Steklov operator P hp : V hp →V ∗ hp can be also represented by P hp = j ∗ hp W j hp + j ∗ hp (cid:18) K (cid:48) + 12 I (cid:19) i hp ( i ∗ hp V i hp ) − i ∗ hp (cid:18) K + 12 I (cid:19) j hp . According to [7], see e.g. [2] for the hp -version, this functional is well-definedand satisfies ∃ c > (cid:104) P hp v hp , v hp (cid:105) Γ ≥ c (cid:107) v hp (cid:107) ˜ H / ( Γ ) ∀ v hp ∈ V hp (37) if min { h, p − } is sufficiently small. Further, we consider the operator E hp : H / ( Γ ) → H − / ( Γ ), reflecting the consistency error in the discretization ofthe Poincar´e-Steklov operator P , defined by E hp := P − P hp = (cid:18) I + K (cid:48) (cid:19) (cid:0) V − − i hp ( i ∗ hp V i hp ) − i ∗ hp (cid:1) (cid:18) I + K (cid:19) . From [29] (see also [7]), the operator E hp is bounded and there exists a constant c > (cid:107) E hp v (cid:107) H − / ( Γ ) ≤ c inf ψ hp ∈W hp (cid:107) V − ( I + K ) v − ψ hp (cid:107) H − / ( Γ ) ∀ v ∈ H / ( Γ ) . Now, we turn to the discretization of the regularized problem ( P ε ). The dis-cretized regularized problem ( P ε,hp ) is: Find u εhp ∈ K Γhp such that for all ∀ v hp ∈ K Γhp (cid:104) P hp u εhp , v hp − u εhp (cid:105) Γ + (cid:104) DJ ε ( u εhp ) , v hp − u εhp (cid:105) Γ C ≥ (cid:104) F , v hp − u εhp (cid:105) Γ . (38)The convergence result for the hp -solution u εhp of ( P ε,hp ) is due to the followingabstract approximation result from [21].Let K be a closed, convex nonempty subset of a reflexive Banach space X and V I ( ψ, f, K ) be the pseudo-monotone variational inequality: Find u ∈ K such that ψ ( u, v ) ≥ (cid:104) f, v − u (cid:105) ∀ v ∈ K . Let T be a directed set. We introduce thefamily {K t } t ∈ T of nonempty, closed and convex sets and admit the followinghypotheses:(H1) If { v t (cid:48) } t (cid:48) ∈ T (cid:48) weakly converges to v in X , v t (cid:48) ∈ K t (cid:48) ( t (cid:48) ∈ T (cid:48) ) for a subnet {K t (cid:48) } t (cid:48) ∈ T (cid:48) of the net {K t } t ∈ T , then v ∈ K .(H2) For any v ∈ K and any t ∈ T there exists v t ∈ K t such that v t → v in X .(H3) ψ t is pseudomonotone for any t ∈ T .(H4) f t → f in V ∗ .(H5) For any nets { u t } and { v t } such that u t ∈ K t , v t ∈ K t , u t (cid:42) u , and v t → v in X it follows thatlim inf t ∈ T ψ t ( u t , v t ) ≤ ψ ( u, v ) . (H6) There exist constants c > d , d ∈ R and α > t ∈ T )such that for some w t ∈ K t with w t → w there holds − ψ t ( u t , w t ) ≥ c (cid:107) u t (cid:107) αV + d (cid:107) u t (cid:107) V + d , ∀ u t ∈ K t , ∀ t ∈ T . Theorem 4 Under the hypotheses ( H - ( H , there exists a solution to theproblem V I ( ψ t , f t , K t ) . Moreover, the family { u t } of solutions to the problem V I ( ψ t , f t , K t ) is bounded in X and there exists a subnet of { u t } that convergesweakly in X to a solution of the problem V I ( ψ, f, K ) . p -adaptive BEM for regularized hemivariational inequalities 19 The hypotheses (H1) and (H2) are due to Glowinski [17] and describe theMosco convergence [1] of the family K t to K , whereas f t in (H4) is a standardapproximation of the linear functional f , for example, by numerical integra-tion. The hypotheses (H5)-(H6) have been verified in [33].We apply Theorem 4 with t = ( ε, h, p ) and show convergence for ε → + ,and either p → ∞ or h → 0. We emphasize that in fact the weak convergenceof { u εhp } to u in H / ( Γ ) can be replaced by the strong one. The completeproof of this result for the h -version of BEM can be found in [33]. In this section we provide an a-priori error estimate for the hp -approximatesolution u εhp of the regularized boundary variational problem ( P ε ) under theregularity assumptions of [29] for Signorini contact, i.e. u ε ∈ H / ( Γ ) and g ∈ H / ( Γ C ). For our more general variational problem ( P ε ) we arrive atthe same convergence order of O ( h / p − / ) as Maischak and Stephan in [29].To this end, we present the following C´ea-Falk approximation lemma for theregularized problem ( P ε ). Lemma 3 Let u ε ∈ K Γ be the solution of the problem ( P ε ) and let u εhp ∈ K Γhp be the solution of the problem ( P ε,hp ) . Assume that α < c p , u ε ∈ H / ( Γ ) , g ∈ H / ( Γ C ) and P u ε − F ∈ L ( Γ ) . Then there exists a constant c = c ( u ε , g, f , t ) > , but independent of h and p such that c (cid:107) u ε − u εhp (cid:107) H / ( Γ ) ≤ (cid:107) E hp ( u ε ) (cid:107) H − / ( Γ ) + inf v ∈K Γ {(cid:107) P u ε − F (cid:107) L ( Γ ) (cid:107) u εhp − v (cid:107) L ( Γ ) + (cid:104) DJ ε ( u ε ) , v − u εhp (cid:105) Γ C } + inf v hp ∈K Γhp (cid:110) (cid:107) u ε − v hp (cid:107) H / ( Γ ) + (cid:107) P u ε − F (cid:107) L ( Γ ) (cid:107) u ε − v hp (cid:107) L ( Γ ) + (cid:104) DJ ε ( u εhp ) , v hp − u ε (cid:105) Γ C (cid:9) . Proof Using the definitions of ( P ε ) and ( P ε,hp ), and estimates similar to [29,Theorem 3], we obtain for all v ∈ K Γ , v hp ∈ K Γhp c P (cid:107) u ε − u εhp (cid:107) H / ( Γ ) ≤ (cid:107) E hp ( u ε ) (cid:107) H − / ( Γ ) + (cid:107) u ε − v hp (cid:107) H / ( Γ ) + (cid:107) P u ε − F (cid:107) L ( Γ ) (cid:0) (cid:107) u ε − v hp (cid:107) L ( Γ ) + (cid:107) u εhp − v (cid:107) L ( Γ ) (cid:1) + D , where we abbreviateD = (cid:104) DJ ε ( u ε ) , v − u ε (cid:105) Γ C + (cid:104) DJ ε ( u εhp ) , v hp − u εhp (cid:105) Γ C . To bound the term D, we use (31) and estimate as follows:D = (cid:104) DJ ε ( u ε ) , v − u εhp (cid:105) Γ C + (cid:104) DJ ε ( u εhp ) , v hp − u ε (cid:105) Γ C + (cid:104) DJ ε ( u ε ) − DJ ε ( u εhp ) , u εhp − u ε (cid:105) Γ C ≤ (cid:104) DJ ε ( u ε ) , v − u εhp (cid:105) Γ C + (cid:104) DJ ε ( u εhp ) , v hp − u ε (cid:105) Γ C + α (cid:107) u ε − u εhp (cid:107) V . Therefore, since α < c P , see Theorem 3, we obtain the assertion. (cid:117)(cid:116) Theorem 5 Let u ε ∈ K Γ be the solution of the problem ( P ε ) and let u εhp ∈K Γhp be the solution of the problem ( P ε,hp ) . Assume that α < c p , u ε ∈ H / ( Γ ) , g ∈ H / ( Γ C ) and P u ε − F ∈ L ( Γ ) . Then there exists a constant c = c ( u ε , g, f , t ) > , but independent of h and p such that (cid:107) u ε − u εhp (cid:107) H / ( Γ ) ≤ ch / p − / . (39) Proof Taking into account the estimates obtained by Maischak and Stephanin [29, Theorem 3] for the consistency error, the approximation error, and for (cid:107) E hp u (cid:107) H − / ( Γ ) , we only need to estimateinf {(cid:104) DJ ε ( u ε ) , v − u εhp (cid:105) Γ C : v ∈ K Γ } (40)and inf {(cid:104) DJ ε ( u εhp ) , v hp − u ε (cid:105) Γ C : v hp ∈ K Γhp } . (41)To estimate (41) we take v hp := I hp u ε ∈ K Γhp , the interpolate of u ε ∈ H / ( Γ ) ⊂ C ( Γ ). From (14a), (cid:104) DJ ε ( u εhp ) , v hp − u ε (cid:105) Γ C = (cid:90) Γ C S x ( u εhp,n , ε )( v hp,n − u ε,n ) ds ≤ c (cid:0) (cid:107) u εhp,n (cid:107) L ( Γ C ) (cid:1) (cid:107) u ε,n − v hp,n (cid:107) L ( Γ C ) . (42)By [5, Theorems 4.2 and 4.5] and by the real interpolation between H ( Γ )and L ( Γ ) there exists a constant C > (cid:107) u ε,n − v hp,n (cid:107) H / ( Γ ) ≤ C h p − (cid:107) u ε (cid:107) H / ( Γ ) . (43)To estimate (40) similarly as [29] we define v ∗ ∈ K Γ ∩ H ( Γ ) by v ∗ := u εhp,t + [ g + inf { u εhp,n − g hp , } ] n on Γ C Γ D γ N u εhp on Γ N , where g hp := I hp g is the interpolate of the gap function g , and γ N is the tracemap onto Γ N .Analogously to (42), we have (cid:104) DJ ε ( u ε ) , v ∗ − u εhp (cid:105) Γ C = (cid:90) Γ C S x ( u ε,n , ε )( v ∗ n − u εhp,n ) ds ≤ c (cid:0) (cid:107) u ε,n (cid:107) L ( Γ C ) (cid:1) (cid:107) v ∗ n − u εhp,n (cid:107) L ( Γ C ) . (44)The elaborate analysis in [29], see the proof of Theorem 3, estimates (31) -(33), gives (cid:107) v ∗ n − u εhp,n (cid:107) L ( Γ C ) ≤ C h / p − / (cid:0) (cid:107) g (cid:107) H / ( Γ C ) + (cid:107) u εhp (cid:107) H / ( Γ C ) (cid:1) . (45)Finally, combining the error estimates for the interpolation (43) and the con-sistency (45) with (42) and (44), respectively, and taking into account theboundedness of (cid:107) u εhp (cid:107) in H / ( Γ C ) (see Theorem 4), we prove the wishedbound for (40) and (41). (cid:117)(cid:116) p -adaptive BEM for regularized hemivariational inequalities 21 Remark 5 Let u ∈ K Γ be the solution of the problem ( P ) and u ε ∈ K Γ bethe solution of the problem ( P ε ) . Taking v = u ε in (5) and v = u in (15), andadding the two resulting inequalities yields (cid:104) P u − P u ε , u ε − u (cid:105) Γ + ϕ ( u , u ε ) + ϕ ε ( u ε , u ) ≥ . Hence, c P (cid:107) u ε − u (cid:107) V ≤ (cid:104) P u ε − P u , u ε − u (cid:105) Γ ≤ ϕ ( u , u ε ) + ϕ ε ( u ε , u ) . We emphasize that, in general, the functional ϕ ε does not approximate thefunctional ϕ arbitrarily close as ε → + . Nevertheless, we have convergenceof the sequence { u ε } of the solutions of the regularized problem ( P ε ) to thesolution u of the boundary hemivariational problem ( P ). In particular, forany v ε (cid:42) v it holds lim sup ε → ϕ ε ( v ε , v ) ≤ ϕ ( v , v ) ∀ v ∈ K Γ , which according to Theorem 4 is a sufficient condition for convergence.However, for convergence rates in ε , ϕ ( u , u ε ) + ϕ ε ( u ε , u ) ≤ c ε β (cid:107) u ε − u (cid:107) V for some positive constants c and β is needed. An estimate of this form can befound in [32] for the approximation of the absolute value function. To be able to split the approximation error into the discretization error ofa simpler variational equation and contributions arising from the constraintson Γ C we introduce the regularized mixed formulation (46a)-(46b), which isequivalent to the regularized problem ( P ε ).Find ( u ε , λ ε ) ∈ V × M ( u ε ) such that (cid:104) P u ε , v (cid:105) Γ + (cid:104) λ ε , v n (cid:105) Γ C = (cid:104) F , v (cid:105) Γ ∀ v ∈ V (46a) (cid:104) µ − λ ε , u εn − g (cid:105) Γ C ≤ ∀ µ ∈ M ( u ε ) (46b)with the set of admissible Lagrange multipliers M ( u ε ) := (cid:8) µ ∈ X ∗ : (cid:104) µ, η (cid:105) Γ C ≥ (cid:104) DJ ε ( u ε ) , η (cid:105) Γ C ∀ η ∈ X, η ≥ Γ C (cid:9) where X = { w | ∃ v ∈ V , v n | Γ C = w } ⊂ H / ( Γ C ) and X ∗ its dual space. Lemma 4 1. Let u ε solve the regularized problem ( P ε ) , then there exists a λ ε ∈ M ( u ε ) such that ( u ε , λ ε ) solves (46) .2. Let ( u ε , λ ε ) solve (46) , then u ε solves ( P ε ) . Proof 1. Define λ ε ∈ X ∗ by (cid:104) λ ε , v n (cid:105) Γ C = (cid:104) F , v (cid:105) Γ − (cid:104) P u ε , v (cid:105) Γ ∀ v ∈ V . (47)Then, by (15), we obtain for all v ∈ K Γ with v n | Γ C = u εn | Γ C − η for0 ≤ η ∈ X (cid:104) λ ε , v n − u εn (cid:105) Γ C = (cid:104) F , v − u ε (cid:105) Γ − (cid:104) P u ε , v − u ε (cid:105) Γ ≤ (cid:104) DJ ε ( u ε ) , v n − u εn (cid:105) Γ C ⇒ (cid:104) λ ε , − η (cid:105) Γ C ≤ (cid:104) DJ ε ( u ε ) , − η (cid:105) Γ C , i.e. λ ε ∈ M ( u ε ). Analogously, we obtain for v n | Γ C = g , v n | Γ C = 2 u εn | Γ C − g (cid:104) λ ε , g − u εn (cid:105) Γ C ≤ (cid:104) DJ ε ( u ε ) , g − u εn (cid:105) Γ C and (cid:104) λ ε , u εn − g (cid:105) Γ C ≤ (cid:104) DJ ε ( u ε ) , u εn − g (cid:105) Γ C , which implies (cid:104) λ ε , u εn − g (cid:105) Γ C = (cid:104) DJ ε ( u ε ) , u εn − g (cid:105) Γ C . It remains to show (46b). From the previous equation we obtain for µ ∈ M ( u ε ) (cid:104) µ − λ ε , u εn − g (cid:105) Γ C = (cid:104) µ, u εn − g (cid:105) Γ C − (cid:104) DJ ε ( u ε ) , u εn − g (cid:105) Γ C ≤ (cid:104) DJ ε ( u ε ) , u εn − g (cid:105) Γ C − (cid:104) DJ ε ( u ε ) , u εn − g (cid:105) Γ C = 0 . 2. There holds trivially DJ ε ( u ε ), 2 λ − DJ ε ( u ε ) ∈ M ( u ε ), and thus (46b)yields (cid:104)− λ ε , u εn − g (cid:105) Γ C ≤ (cid:104)− DJ ε ( u ε ) , u εn − g (cid:105) Γ C , (cid:104) λ ε , g − u εn (cid:105) Γ C ≤ (cid:104) DJ ε ( u ε ) , g − u εn (cid:105) Γ C ⇒ (cid:104) λ ε , u εn − g (cid:105) Γ C = (cid:104) DJ ε ( u ε ) , u εn − g (cid:105) Γ C . For v ∈ K Γ , i.e. v n | Γ C = g − η with 0 ≤ η ∈ X , we therefore obtain (cid:104) λ ε , v n − u εn (cid:105) Γ C = (cid:104) λ ε , g − u εn (cid:105) Γ C + (cid:104) λ ε , − η (cid:105) Γ C = (cid:104) DJ ε ( u ε ) , g − u εn (cid:105) Γ C + (cid:104) λ ε , − η (cid:105) Γ C ≤ (cid:104) DJ ε ( u ε ) , g − u εn (cid:105) Γ C + (cid:104) DJ ε ( u ε ) , − η (cid:105) Γ C = (cid:104) DJ ε ( u ε ) , v n − u εn (cid:105) Γ C . It remains to show that u εn − g ≤ 0. Assume there exists a measurable set inwhich u εn − g > 0. Then, from (46b) we obtain for µ = DJ ε ( u ε )+ χ u εn − g> ∈ M ( u ε ) with indicator function χ ≥ (cid:104) µ − λ ε , u εn − g (cid:105) Γ C = (cid:104) µ − DJ ε ( u ε ) , u εn − g (cid:105) Γ C = (cid:10) χ u εn − g> , u εn − g (cid:11) Γ C > , i.e. u ε ∈ K Γ . (cid:117)(cid:116) Given the discrete solution u εhp ∈ K Γhp to ( P ε,hp ), we reconstruct λ εhp ∈ span { ψ i } Mi =1 such that (cid:10) λ εhp , v n (cid:11) Γ C = (cid:104) F , v (cid:105) Γ − (cid:10) P hp u εhp , v (cid:11) Γ ∀ v ∈ V hp (48) p -adaptive BEM for regularized hemivariational inequalities 23 by solving a potentially over-constrained system of linear equations for anarbitrary choice of basis { ψ } . Following the Braess trick [6] as e.g. in [4], wedefine the auxiliary problem z ∈ V : (cid:104) P z , v (cid:105) Γ = (cid:104) F , v (cid:105) Γ − (cid:10) λ εhp , v n (cid:11) Γ C ∀ v ∈ V . (49)Subtracting (46a) and (49) yields (cid:104) P ( u ε − z ) , v (cid:105) Γ = (cid:10) λ εhp − λ ε , v n (cid:11) Γ C ∀ v ∈ V (50)and additionally with the continuous inf-sup condition [11, Theorem 3.2.1] thisyields (see [4]) (cid:13)(cid:13) λ εhp − λ ε (cid:13)(cid:13) X ∗ ≤ Cβ (cid:107) u ε − z (cid:107) V ≤ Cβ (cid:13)(cid:13) u ε − u εhp (cid:13)(cid:13) V + Cβ (cid:13)(cid:13) u εhp − z (cid:13)(cid:13) V (51)with inf-sup constant β > 0. See [11, Theorem 3.2.1] for a proof of the inf-supcondition for the difficult case when ¯ Γ C ∩ ¯ Γ D = ∅ , i.e. X ∗ = ˜ H − / ( Γ C ). Theorem 6 Under the assumption (30) and if S x ( · , ε ) is Lipschitz continu-ous, then there exists a constant C independent of h and p such that for ς > arbitrary ( c P − α − ς ) (cid:107) u ε − u εhp (cid:107) V ≤ (cid:18) Cς + 1 (cid:19) (cid:107) z − u εhp (cid:107) V + 14 ς (cid:13)(cid:13) ( λ εhp − DJ ε ( u εhp )) − (cid:13)(cid:13) X ∗ + C (cid:18) ς + 1 β + 1 ςβ (cid:19) (cid:13)(cid:13) ( u εhp,n − g ) + (cid:13)(cid:13) X − (cid:10) ( λ εhp − DJ ε ( u εhp )) + , ( u εhp,n − g ) − (cid:11) Γ C with λ εhp satisfying (48) . Proof From the Lipschitz continuity it follows that (cid:10) DJ ε ( u εhp ) − DJ ε ( u ε ) , v n (cid:11) Γ C = (cid:90) Γ C (cid:0) S x ( u εhp,n ) − S x ( u εn ) (cid:1) v n ds ≤ Lip (cid:107) u εhp,n − u εn (cid:107) L ( Γ C ) (cid:107) v n (cid:107) L ( Γ C ) and thus, by Young’s inequality ( ς > (cid:10) DJ ε ( u εhp ) − DJ ε ( u ε ) , ( u εhp,n − g ) + (cid:11) Γ C ≤ ς (cid:107) u εhp − u ε (cid:107) L ( Γ C ) + C ς (cid:107) ( u εhp,n − g ) + (cid:107) L ( Γ C ) . (52)From the conformity V hp ⊂ V we obtain with (50) c P (cid:107) u ε − u εhp (cid:107) V ≤ (cid:10) P ( u ε − u εhp ) , u ε − u εhp (cid:11) Γ = (cid:10) P ( z − u εhp ) , u ε − u εhp (cid:11) Γ + (cid:10) λ εhp − λ ε , u εn − u εhp,n (cid:11) Γ C ≤ C P (cid:107) z − u εhp (cid:107) V (cid:107) u ε − u εhp (cid:107) V + (cid:10) λ εhp − λ ε , u εn − u εhp,n (cid:11) Γ C . From complementarity (cid:104) λ ε , u εn − g (cid:105) Γ C = (cid:104) DJ ε ( u ε ) , u εn − g (cid:105) Γ C , λ ε ∈ M ( u ε ), v = v + + v − = max { , v } + min { , v } , u εn − g ≤ Γ C we obtain (cid:10) λ εhp − λ ε , u εn − u εhp,n (cid:11) Γ C = (cid:10) λ εhp − DJ ε ( u εhp ) , u εn − u εhp,n (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − λ ε , u εn − g (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − λ ε , g − u εhp,n (cid:11) Γ C = (cid:10) λ εhp − DJ ε ( u εhp ) , u εn − u εhp,n (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − DJ ε ( u ε ) , u εn − g (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − λ ε , g − u εhp,n (cid:11) Γ C = (cid:10) λ εhp − DJ ε ( u εhp ) , u εn − u εhp,n (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − DJ ε ( u ε ) , u εn − u εhp,n (cid:11) Γ C + (cid:10) λ ε − DJ ε ( u ε ) , u εhp,n − g (cid:11) Γ C ≤ (cid:10) ( λ εhp − DJ ε ( u εhp )) + , u εn − g + g − u εhp,n (cid:11) Γ C + (cid:10) ( λ εhp − DJ ε ( u εhp )) − , u εn − u εhp,n (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − DJ ε ( u ε ) , u εn − u εhp,n (cid:11) Γ C + (cid:10) λ ε − DJ ε ( u ε ) , ( u εhp,n − g ) + (cid:11) Γ C ≤ (cid:10) ( λ εhp − DJ ε ( u εhp )) + , g − u εhp,n (cid:11) Γ C + (cid:10) λ ε − DJ ε ( u ε ) , ( u εhp,n − g ) + (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − DJ ε ( u ε ) , u εn − u εhp,n (cid:11) Γ C + (cid:10) ( λ εhp − DJ ε ( u εhp )) − , u εn − u εhp,n (cid:11) Γ C . Since (cid:10) ( λ εhp − DJ ε ( u εhp )) + , g − u εhp,n (cid:11) Γ C + (cid:10) λ ε − DJ ε ( u ε ) , ( u εhp,n − g ) + (cid:11) Γ C = − (cid:10) ( λ εhp − DJ ε ( u εhp )) + , ( u εhp,n − g ) + + ( u εhp,n − g ) − (cid:11) Γ C + (cid:10) λ ε − DJ ε ( u ε ) − ( λ εhp − DJ ε ( u εhp )) + ( λ εhp − DJ ε ( u εhp )) + , ( u εhp,n − g ) + (cid:11) Γ C + (cid:10) ( λ εhp − DJ ε ( u εhp )) − , ( u εhp,n − g ) + (cid:11) Γ C ≤ (cid:10) λ ε − λ εhp , ( u εhp,n − g ) + (cid:11) Γ C + (cid:10) DJ ε ( u εhp ) − DJ ε ( u ε ) , ( u εhp,n − g ) + (cid:11) Γ C − (cid:10) ( λ εhp − DJ ε ( u εhp )) + , ( u εhp,n − g ) − (cid:11) Γ C , application of the Cauchy-Schwarz inequality, Young’s inequality with ς > c P − α − ς ) (cid:107) u ε − u εhp (cid:107) V ≤ (cid:18) Cς + 1 (cid:19) (cid:107) z − u εhp (cid:107) V + 14 ς (cid:13)(cid:13) ( λ εhp − DJ ε ( u εhp )) − (cid:13)(cid:13) X ∗ + C (cid:18) ς + 1 β + 1 ςβ (cid:19) (cid:13)(cid:13) ( u εhp,n − g ) + (cid:13)(cid:13) X − (cid:10) ( λ εhp − DJ ε ( u εhp )) + , ( u εhp,n − g ) − (cid:11) Γ C . (cid:117)(cid:116) The a-posteriori error estimate decomposes into the discretization error of avariational equality (cid:107) z − u εhp (cid:107) V , which can be further estimated by e.g. residualerror estimates [8] or bubble error estimates, e.g. [4], and violation of the con-sistency condition (cid:13)(cid:13)(cid:13) ( λ εhp − DJ ε ( u εhp )) − (cid:13)(cid:13)(cid:13) X ∗ , violation of the non-penetrationcondition (cid:13)(cid:13)(cid:13) ( u εhp,n − g ) + (cid:13)(cid:13)(cid:13) X and violation of the complementarity condition − (cid:68) ( λ εhp − DJ ε ( u εhp )) + , ( u εhp,n − g ) − (cid:69) Γ C . Localizing an approximation of the p -adaptive BEM for regularized hemivariational inequalities 25 global a-posteriori error estimate gives rise to the following solve-mark-refinealgorithm for hp -adaptivity. Algorithm 1 (Solve-mark-refine algorithm for hp -adaptivity)1. Choose initial discretization T h and p , steering parameters θ ∈ (0 , 1) and δ ∈ (0 , k = 0 , , , . . . do(a) solve regularized discrete problem ( P varepsilon,hp ).(b) compute discrete Lagrange multiplier by (48)(c) compute local error indicators Ξ to current solution.(d) mark all elements T ∈ NN := argmin {|{ (cid:98) N ⊂ T h : (cid:88) T ∈ (cid:98) N Ξ ( T ) ≥ θ (cid:88) T ∈T h Ξ ( T ) }|} (e) estimate local analyticity, i.e. compute Legendre coefficients of v hp | T ( Ψ T ( x )) = p T (cid:88) j =0 a i L i ( x ) , a i = 2 i + 12 (cid:90) − v hp | T ( Ψ T ( x )) L i ( x ) dx and use a least square approach to compute the slop m of | log | a i || = mi + b , for each direction of u hp on Γ Σ , of ψ hp on Γ D , respectively. If e − m ≤ δ for all directions then p -refine, else h -refine marked element E . If p E = 0 always p -refine to have a decision basis next time.(f) refine marked elements based on the decision in 2(e). For the adaptivity algorithm 1 we use D¨orfler marking with marking parameter θ = 0 . 3, and the Legendre expansion strategy [24] with δ = 0 . h - and p -refinement. The error of the auxiliary problem is estimatedby a bubble error estimate as in e.g. [4] and the non-localized norms in the de-lamination related error contributions are approximated by scaled L -norms,namely hp − (cid:107) ( u εhp,n − g ) + (cid:107) L ( Γ C ) and h − p (cid:107) ( λ εhp − DJ ε ( u εhp )) − (cid:107) L ( Γ C ) . The in-tegral (cid:104) DJ ε ( u εhp ) , v hp − u εhp (cid:105) Γ C is computed by a composite Gauss-Quadraturewith (cid:100) p T (cid:101) + 17 quadrature points per element. If not mentioned otherwise theregularization parameter is ε = 10 − .For the numerical experiments we choose Ω = (0 , / , Γ D = { }× [0 , / Γ C = (0 , / ×{ } , Γ N = ∂Ω \ ( Γ D ∪ Γ C ). The material parameters are E = 5, ν = 0 . f ≡ t = 0 . 25 on [1 / , / × { } and zero elsewhere, g = 0. Thedelamination law is given via f ( u n ( x )) = min { g ( g ( x ) − u n ( x )) , g ( g ( x ) − u n ( x )) , g ( g ( x ) − u n ( x )) } = − max {− g ( − u n ( x )) , − g ( − u n ( x )) , − g ( − u n ( x )) } g-u n - S x ( u n ) Fig. 5 Regularized delamination law S x for ε = 10 − with g ( y ) = A t y , g ( y ) = b ( y − t ) + d , g ( y ) = d and parameters A = 0 . , A = 0 . , t = 0 . , t = 0 . ,b = A t , d = A t , d = b ( t − t ) + d . The regularized delamination law S x with regularization parameter ε = 10 − is plotted in Figure 5. The characteristic saw tooth shape is already present,but the absolute value in the tips and the slope approximating the jump arestill noticeable coarse approximated.The discrete Lagrange multiplier λ εhp is obtain by solving (48) where ψ i are discontinuous, piecewise polynomials on Γ C on a one time coarsened mesh( H = 2 h ) with polynomial degree reduced by one ( q = p − 1) compared tothe mesh and polynomial degree distribution of u εhp . Figure 6 displays thedeformation of the rectangle and the normal stresses on Γ C obtained fromthe lowest order uniform h -method with 16384 elements and regularizationparameter ε = 10 − . The normal stress on Γ C , Figure 6 (b), reflects thedelamination law from Figure 5 well.Figure 7 displays the reduction of the error in ( u , λ ) and of the error es-timate. Since the exact solution is not known, we compute the error approx-imately by (cid:107) u fine − u hp (cid:107) P and (cid:107) λ fine − λ hp (cid:107) V , with norms induced by thePoincar´e-Steklov operator, single layer potential acting on Γ C , respectively.The pair ( u fine , λ fine ) is a very fine (last) approximation for each sequence ofdiscretization. These error expressions are only meaningful if the distance issufficiently large, which at the end of the different discretization sequences isnot true and, therefore, are then omitted. Hence, the error curves are ”shorter”then the error estimate curves.The uniform h -version with p = 1 exhibits an experimental order of con-vergence (eoc) of 0.65 for the error and of 0.53 for the error estimate, all p -adaptive BEM for regularized hemivariational inequalities 27 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.600.10.20.30.40.50.6 (a) Deformation Γ C no r m a l c on t a c t s t r e ss -0.0100.010.020.030.040.05 (b) Contact stresses Fig. 6 Uniform h -version, p=1, 16384 elements with 4096 elements on Γ C Degrees of Freedom -4 -3 -2 -1 error approx., uniform h, p=1error estimate, uniform h, p=1error approx., h-adaptiveerror estimate, h-adaptiveerror approx., hp-adaptiveerror estimate, hp-adaptive Fig. 7 (cid:113) (cid:107) u fine − u hp (cid:107) P + (cid:107) λ fine − λ hp (cid:107) V and error estimate for different families ofdiscrete solutions, ε = 10 − w.r.t. degrees of freedom (dof). Hence, the error estimate is only reliable butnot efficient. For the h -adaptive scheme with p = 1, the eoc is 0.84 for the errorand only 0.56 for the error estimate over the last 20, 10 iterations, respectively.In case of hp -adaptivity the eoc appears to be exponential for the error and1.52 for the error estimate also over the last 10 iterations.Looking at the single contributions of the error estimate plotted in Figure 8we find that in all cases the complementarity contribution is the dominant er-ror contribution with the lowest order of convergence and that the consistencycontribution has always a significantly higher order of convergence. Given thenature of this numerical experiment, the violation of the non-penetration con-dition is always zero and, therefore, is omitted in these plots. For the uniform h -version, Figure 8 (a), we have that the Neumann and V − contributions con-verge at 0 . 62 almost as fast as the error approximation. Given the nature, that λ h, is piecewise constant with double mesh size, − u h, · n is linearly increasing Degrees of Freedom -7 -6 -5 -4 -3 -2 -1 Error estimateNeumann contr.V -1 contr.Consistency contr.Complementarity contr. (a) uniform h -version with p = 1 Degrees of Freedom -7 -6 -5 -4 -3 -2 -1 Error estimateNeumann contr.V -1 contr.Consistency contr.Complementarity contr. (b) h -adaptive with p = 1 Degrees of Freedom -6 -5 -4 -3 -2 -1 Error estimateNeumann contr.V -1 contr.Consistency contr.Complementarity contr. (c) hp -adaptive Fig. 8 Error contributions of the error estimates and thus − S x ( u h, · n ) in the simplest case as well, the complementarity con-tribution is zero on the left element and non-zero on the right element. Thisoscillatory case occurs over a large fraction of Γ C and thus the h -adaptivescheme performs almost uniform mesh refinements on Γ C , c.f. Figure 9(a).Hence, the little gain in h -adaptivity, Figure 8 (b), for the error estimate. Thisoscillatory effect can be reduced by increasing the polynomial degree as in hp -adaptivity, which balances out the single error contributions, Figure 8 (c)and is thus capable of producing more localized refinements, c.f. Figure 9(b).Finally, Figure 10 shows the influence of the regularization parameter onthe errors (cid:107) u fine − u εh, (cid:107) P and (cid:107) λ fine − λ εh, (cid:107) V where ( u fine , λ fine ) is theGalerkin solution to a uniform mesh with 16384 elements, p = 1 and regular-ization parameter ε = 2 . · − . Here, the pair ( u εh, , λ εh, ) is computed on auniform mesh with 8192 elements and the regularization parameter is varied.We find that for very large ε , the regularized delamination law has lost itssaw tooth characteristic and small variations of ε have no noticeable influence.For very small ε the discretization error is dominant which sets in for u much p -adaptive BEM for regularized hemivariational inequalities 29 (a) h -adaptive, mesh nr. 10 (inner), nr. 20(outer) (b) hp -adaptive, mesh nr. 15 (inner), nr. 25(outer) Fig. 9 Adaptively generated meshes Regularization parameter -6 -5 -4 -3 -2 -1 -6 -5 -4 -3 -2 111 -0.03 error uerror λ error estimate0.155 · ǫ Fig. 10 Error and error estimate for different regularization parameters ε , uniform h -version8192 elements, p = 1 earlier than for λ . In the intermediate range the eoc w.r.t. ε is 0.75 for u ,whereas λ displays an (asymptotic) eoc of 1.0 w.r.t. ε . The error estimate forthe regularized problem increases in ε slightly. Acknowledgment The authors are grateful to J. Gwinner for the helpfuldiscussions and comments. References 1. Attouch, H.: Variational convergence for functions and operators. Pitman, Boston (1984)0 Nina Ovcharova, Lothar Banz2. Banz, L., Gimperlein, H., Issaoui, A., Stephan, E.P.: Stabilized mixed hp -BEM for fric-tional contact problems in linear elasticity. Numer. Math. (2016). DOI: 10.1007/s00211-016-0797-y3. Banz, L., Stephan, E.P.: A posteriori error estimates of hp -adaptive IPDG-FEM forelliptic obstacle problems. Appl. Numer. Math. , 76–92 (2014)4. Banz, L., Stephan, E.P.: On hp-adaptive BEM for frictional contact problems in linearelasticity. Comput. Math. Appl. (7), 559–581 (2015)5. Bernardi, C., Maday, Y.: Spectral methods. In: Handb. Numer. Anal., vol. V, pp.209–485. North-Holland, Amsterdam (1997)6. Braess, D.: A posteriori error estimators for obstacle problems–another look. Numer.Math. (3), 415–421 (2005)7. Carstensen, C.: Interface problem in holonomic elastoplasticity. Math. Methods Appl.Sci. (11), 819–835 (1993)8. Carstensen, C.: A posteriori error estimate for the symmetric coupling of finite elementsand boundary elements. Computing (4), 301–322 (1996)9. Carstensen, C., Funken, S.A., Stephan, E.P.: On the adaptive coupling of FEM andBEM in 2–d–elasticity. Numer. Math. (2), 187–221 (1997)10. Carstensen, C., Gwinner, J.: FEM and BEM coupling for a nonlinear transmissionproblem with Signorini contact. SIAM J. Numer. Anal. (5), 1845–1864 (1997)11. Chernov, A.: Nonconforming boundary elements and finite elements for interface andcontact problems with friction - hp-version for mortar, penalty and Nitsche’s methods.Ph.D. thesis, Universit¨at Hannover (2006)12. Chernov, A., Maischak, M., Stephan, E.P.: hp-Mortar boundary element method fortwo-body contact problems with friction. Math. Methods Appl. Sci. (17), 2029–2054(2008)13. Clarke, F.H.: Optimization and nonsmooth analysis. John Wiley and Sons (1983)14. Costabel, M.: Boundary integral operators on Lipschitz domains: Elementary results.SIAM J. Math. Anal. (3), 613–626 (1988)15. Costabel, M., Stephan, E.P.: Coupling of finite and boundary element methods for anelastoplastic interface problem. SIAM J. Numer. Anal. (5), 1212–1226 (1990)16. Dao, M., Gwinner, J., Noll, D., Ovcharova, N.: Nonconvex bundle method with applica-tion to a delamination problem. Computational Optimization and Applications (2016).DOI: 10.1007/s10589-016-9834-017. Glowinski, R.: Numerical methods for nonlinear variational problems. Springer Seriesin Computational Physics (1984)18. Guediri, H.: On a boundary variational inequality of the second kind modelling a frictionproblem. Math. Methods Appl. Sci. (2), 93–114 (2002)19. Gwinner, J.: Nichtlineare Variationsungleichungen mit Anwendungen. Ph.D. thesis,Universit¨at Mannheim (1978)20. Gwinner, J.: hp-FEM convergence for unilateral contact problems with Tresca frictionin plane linear elastostatics. J. Comput. Appl. Math. , 175–184 (2013)21. Gwinner, J., Ovcharova, N.: From solvability and approximation of variational inequal-ities to solution of nondifferentiable optimization problems in contact mechanics. Opti-mization (ahead-of-print), 1–20 (2015)22. Gwinner, J., Stephan, E.P.: A boundary element procedure for contact problems inlinear elastostatics. ESAIM Math. Model. Numer. Anal. (4), 457 – 480 (1993)23. Han, H.D.: A direct boundary element method for Signorini problems. Math. Comp. (191), 115–128 (1990)24. Houston, P., S¨uli, E.: A note on the design of hp-adaptive finite element methods forelliptic partial differential equations. Comput. Methods Appl. Mech. Engrg. (2-5),229–243 (2005)25. Hsiao, G.C., Wendland, W.L.: Boundary integral equations. Springer (2008)26. Kleiber, M., Borkowski, A.: Handbook of computational solid mechanics: survey andcomparison of contemporary methods. Springer (1998)27. Krebs, A., Stephan, E.: A p-version finite element method for nonlinear elliptic varia-tional inequalities in 2D. Numer. Math. (3), 457–480 (2007)28. Maischak, M., Stephan, E.P.: A FEM–BEM coupling method for a nonlinear trans-mission problem modelling Coulomb friction contact. Comput. Methods Appl. Mech.Engrg. (2), 453–466 (2005) p -adaptive BEM for regularized hemivariational inequalities 3129. Maischak, M., Stephan, E.P.: Adaptive hp-versions of BEM for Signorini problems.Appl. Numer. Math. (3), 425–449 (2005)30. Maischak, M., Stephan, E.P.: Adaptive hp-versions of boundary element methods forelastic contact problems. Comput. Mech. (5), 597–607 (2007)31. Nesemann, L., Stephan, E.P.: Numerical solution of an adhesion problem with FEMand BEM. Appl. Numer. Math. (5), 606–619 (2012)32. Ovcharova, N.: Regularization Methods and Finite Element Approximation of Hemi-variational Inequalities with Applications to Nonmonotone Contact Problems. Ph.D.thesis, Universit¨at der Bundeswehr M¨unchen (2012)33. Ovcharova, N.: On the coupling of regularization techniques and the boundary elementmethod for a hemivariational inequality modelling a delamination problem. Math. Meth-ods Appl. Sci. (2015). ArXiv 1603.0509134. Ovcharova, N., Gwinner, J.: A study of regularization techniques of nondifferentiableoptimization in view of application to hemivariational inequalities. J. Optim. TheoryAppl. (3), 754–778 (2014)35. Sloan, I.H., Spence, A.: The Galerkin method for integral equations of the first kindwith logarithmic kernel: Theory. IMA J. Numer. Anal. (1), 105–122 (1988)36. Steinbach, O.: Numerische N¨aherungsverfahren f¨ur elliptische Randwertprobleme: FiniteElemente und Randelemente. Vieweg + Teubner (2003)37. Stephan, E.P.: The hp-version of BEMfast convergence, adaptivity and efficient precon-ditioning. J. Comput. Math. (2-3), 348–359 (2009)38. Wetzel, M., Holtmannsp¨otter, J., Gudladt, H.J., Czarnecki, J.: Sensitivity of doublecantilever beam test to surface contamination and surface pretreatment. Int. J. Adhes.Adhes.46