On the SCD semismooth* Newton method for generalized equations with application to a class of static contact problems with Coulomb friction

In the paper, a variant of the semismooth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{*}$$\end{document}∗ Newton method is developed for the numerical solution of generalized equations, in which the multi-valued part is a so-called SCD (subspace containing derivative) mapping. Under a rather mild regularity requirement, the method exhibits (locally) superlinear convergence behavior. From the main conceptual algorithm, two implementable variants are derived whose efficiency is tested via a generalized equation modeling a discretized static contact problem with Coulomb friction.


Introduction
Consider the generalized equation (GE) 0 ∈ H(x) := f (x) + Q(x), (1.1) where f : R n → R n is continuously differentiable and Q : R n ⇒ R n is a set-valued mapping with a closed graph.For the numerical solution of (1.1) various methods are available, including the semismooth * Newton method developed in [13].In this method, the approximation/linearization of the multi-valued term in (1.1) is performed on the basis of either the graph of the respective strict derivative or the limiting coderivative.In each Newton step, one has to solve a linear square system with a non-singular matrix.The method thus differs both from the approach of Josephy ( [20], [21]), where the multi-valued part is not approximated at all, and from the Newton-type methods in [1] and [6], where the multi-valued term is approximated in a different way.
In [14] the authors have shown that for a class of the so-called SCD (subspace containing derivatives) mappings the semismooth * Newton can be improved.In particular, at these mappings, we dispose at each point with linear subspaces belonging to the graphs of the above-mentioned generalized derivatives, which generate the linear systems in the Newton step in a straightforward way.Moreover, the "regularity" requirement, needed to ensure the (locally) superlinear convergence, could have been substantially relaxed.In [15] this so-called SCD semismooth * Newton method has been implemented in a class of variational inequalities (VIs) of the second kind that includes, among various problems of practical importance, also a class of discretized contact problems with Tresca friction ( [16]).The very good performance of the new method, when applied to those problems, has led us to consider more complicated problems, in which Q does not amount to the subdifferential of a convex function.This clearly requires a generalization of the theory of [15].As a concrete representative problem of this type, we have chosen a discretized static contact problem, where the Tresca friction is replaced by the (physically more realistic) Coulomb friction.
Starting with the pioneering paper [27] there are many papers and a comprehensive monograph [11] devoted to static, quasi-static and dynamic contact problems with Coulomb friction for various types of material of the bodies in contact.Concerning the static contact of two elastic bodies or an elastic body with a rigid obstacle, it is known ( [27]) that this problem has a (not necessarily unique) solution whenever the friction coefficient belongs to the interval (0, b], where b > 0 is a bound depending on the Poisson constant.This is a great difference from the corresponding discretized problems where, for small (mesh-dependent) values of the friction coefficient, one has to do with a unique solution.However, when the friction coefficient increases, discretized models may allow multiple solutions, as shown, e.g., in [30], [26], For simplicity, we consider only the contact of an elastic body with a rigid obstacle (Signorini problem with friction).From the algebraic point of view, each two-body problem can be rewritten formally as a one-body problem, see [33].Such problems can be modeled as quasi-variational inequalities (QVIs) expressed in terms of the so-called dual variables (having the physical meaning of stresses).This enables us to employ a variety of methods developed for the numerical solution of QVIs, cf., e.g., [18,Chapter 5], [12].As shown in [23], [28], to these methods we can also count the classical semismooth Newton method ( [25], [29]).Another approach has been used in [2], [3], where a solution is computed as a fixed point of a mapping generated by solving the corresponding contact problems with the Tresca friction.These problems can be solved, for example, by a specially tailored minimization routine; see [24].In this paper, we use a new discrete model formulated in displacements, which is obtained by a modification of the GE used in [3].We thus work with a purely "primal" model, well suited for a direct application of the SCD semismooth * Newton method.
The plan of the paper is as follows: In Section 2 we recall some standard notions from variational analysis, which are extensively used throughout the paper.Section 3 is devoted to the SCD mappings; we list their basic properties and provide some calculus rules, indispensable in the construction of the SCD semismooth * Newton method in Section 4 and its subsequent implementation.In Section 4 we present the main conceptual algorithm which is thereafter, in Section 5, implemented to the numerical solution of GE (1.1).As a result, we obtain an efficient tool, applicable to a wide range of equilibrium models including VIs of the first and second kind, hemivariational inequalities ( [19]) and many others.The first part of Section 6 is devoted to the construction of the new model of the discrete contact problem with Coulomb friction mentioned above (Section 6.1).In Section 6.2, it is shown that the multifunction, which arises in the respective GE, is an SCD mapping and possesses the semismooth * property.This finally enables us to specialize the formulas for the approximation and Newton step, developed in Section 5, for the GE considered.The results of the numerical experiments are presented in Section 7.
Our notation is basically standard.Given a linear subspace L ⊆ R n , L ⊥ denotes its orthogonal complement.For an element u ∈ R n , u denotes its Euclidean norm, B ε (u) denotes the closed ball around u with radius ε and S R n stands for the unit sphere in R n .For a matrix A, rge A signifies its range.To avoid possible confusion, in some situations the dimension of a unit matrix I will be indicated by a subscript (I n ) .Given a set Ω ⊂ R s , we define the distance from a point x to Ω by d Ω (x) := dist(x, Ω) := inf{ y − x y ∈ Ω}, the respective indicator function is denoted by δ Ω and Ω x → x means convergence within Ω.When a mapping F : R n → R m is differentiable at x, we denote by ∇F(x) its Jacobian.

Preliminaries
Throughout the paper, we will frequently use the following basic notions of modern variational analysis.
Definition 2.1.Let A be a set in R s , x ∈ A and A be locally closed around x. Then (i) The tangent (contingent, Bouligand) cone to A at x is given by The set A is geometrically derivable at x if every tangent vector u to A at x is derivable.
is the regular (Fréchet) normal cone to A at x, and is the limiting (Mordukhovich) normal cone to A at x.
In this definition "Lim sup" stands for the Painlevé-Kuratowski outer (upper) set limit, see, e.g., [32].The above listed cones enable us to describe the local behavior of set-valued maps via various generalized derivatives.Let F : R n ⇒ R m be a (set-valued) mapping with the domain and the graph dom Definition 2.2.Consider a (set-valued) mapping F : R n ⇒ R m and let gph F be locally closed around some ( x, ȳ) ∈ gph F.
(ii) The multifunction D * F( x, ȳ) : R m ⇒ R n , defined by is called the limiting coderivative of F at ( x, ȳ).
Let us now recall the following regularity notions.
1. F is said to be metrically subregular at ( x, ȳ) if there exists a real κ ≥ 0 along with some neighborhood X of x such that dist(x, 2. F is said to be strongly metrically subregular at ( x, ȳ) if it is metrically subregular at ( x, ȳ) and there exists a neighborhood X of x such that F −1 ( ȳ) ∩ X = { x}.
3. F is said to be metrically regular around ( x, ȳ) if there is some κ ≥ 0 together with neighborhoods X of x and Y of ȳ such that dist(x, 4. F is said to be strongly metrically regular around ( x, ȳ) if it is metrically regular around ( x, ȳ) and F −1 has a single-valued localization around ( ȳ, ȳ), i.e., there are open neighborhoods Y of ȳ, X of x and a mapping h It is easy to see that the strong metric regularity around ( x, ȳ) implies the strong metric subregularity at ( x, ȳ) and the metric regularity around ( x, ȳ) implies the metric subregularity at ( x, ȳ).To check the metric regularity one often employs the so-called Mordukhovich criterion, see, e.g.[32,Theorem 9.43], according to which this property around ( x, ȳ) is equivalent to the condition 0 ∈ D * F( x, ȳ)(y * ) ⇒ y * = 0. (2.2) Another useful characterization of metric regularity in terms of the graphical derivative is provided by the so-called Aubin-criterion by Dontchev et al [9].For pointwise characterizations of the other stability properties from Definition 2.3, the reader is referred to [14,Theorem 2.7].
In this preparatory section, we end with a definition of the semismooth * property, which paved the way for both the semismooth * Newton method in [13] and the SCD semismooth * Newton method in [14].
3 On SCD mappings

Basic properties
In this section, we wish to recall the basic definitions and characteristics of the SCD property introduced in the recent paper [14].
In what follows, we denote by Z n the metric space of all n-dimensional subspaces of R 2n equipped with the metric where P L i is the symmetric 2n × 2n matrix representing the orthogonal projection onto L i , i = 1, 2. Sometimes, we will also work with bases for the subspaces L ∈ Z n .Let M n denote the collection of all 2n × n matrices with full column rank n and define for L ∈ Z n the set i.e., the columns of Z ∈ M (L) create a basis for L.
We treat every element of R 2n as a column vector.To keep the notation simple, we write (u, v) instead of u v ∈ R 2n when this does not cause confusion.
Let L ∈ Z n and consider Z ∈ M (L).Then we can divide Z into two n × n matrices A and B and write Furthermore, for every L ∈ Z n we can define the adjoint space It can be shown that . Thus, the mapping L → L * defines an isometry on Z n .Definition 3.1.Consider a mapping F : R n ⇒ R n .
1. We call F graphically smooth of dimension n at (x, y) ∈ gph F, if T gph F (x, y) = gph DF(x, y) ∈ Z n .In addition, we denote by O F the set of all points where F is graphically smooth of dimension n.
2. We associate with F the four mappings S F : gph 3. We say that F has the SCD (subspace containing derivative) property at (x, y) ∈ gph F, if S * F(x, y) = / 0. F is said to have the SCD property around (x, y) ∈ gph F, if there is a neighborhood W of (x, y) such that F has the SCD property at every (x , y ) ∈ gph F ∩W .Finally, we call F an SCD mapping if F has the SCD property at every point of its graph.
Since L → L * is an isometry on Z n and (L * ) * = L, the mappings S * F and S F are related through S * F(x, y) = {L * L ∈ S F(x, y)}, S F(x, y) = {L * L ∈ S * F(x, y)}. (3.4) The name SCD property is motivated by the following statement.In the recent paper [14] one can find several calculus rules to work with the SCD property including the next result.Proposition 3.3 (cf.[14,Proposition 3.15]).Let F : R n ⇒ R n have the SCD property at (x, y) ∈ gph F and let h : U → R n be continuously differentiable at x ∈ U where U ⊆ R n is open.Then F + h has the SCD property at (x, y + h(x)) and Note that these sum rules are also valid at the points (x, y) ∈ gph F where F does not have the SCD property: In this case, we simply have S (F + h)(x, y + h(x)) = S * (F + h)(x, y + h(x)) = S F(x, y) = S * F(x, y) = / 0. In addition, we will need some calculus rules for the Cartesian product of mappings.Consider the mapping F : where each multifunction F i : R n i ⇒ R n i , i = 1, . . ., p, has a closed graph.Note that where the first equation follows from the identity gph F = (x 1 , . . ., x p ), (y 1 , . . ., y p ) (x Proof.In order to prove the first assertion, let (x, y) ∈ O F , let i ∈ {1, . . ., p} be arbitrarily fixed and consider the set Then it readily follows from Definition 2.
and from (3.9) we deduce The second statement follows from the fact that under the stated assumption, inclusion (3.9) holds with equality, cf.[17, Proposition 1].Proposition 3.5.Let F be given by (3.7) and assume that all the mappings F i , i = 1, . . ., p, are SCD mappings.If all, but at most one, of the mappings F i , i = 1, . . ., p, have the property that gph F i is geometrically derivable at every point (x i , y i ) ∈ O F i , then F is an SCD mapping, and Moreover, for every (x, y) = (x 1 , . . ., x p ), (y 1 , . . ., y p ) ∈ gph F there holds The equality holds in the inclusions (3.10), (3.11) and (3.12) if, in addition, all but at most one of the mappings F i , i = 1, . . ., p, have the following property: For every (x i , y i ) ∈ gph F i such that T gph F i (x i , y i ) is a subspace, the dimension of this subspace is n i .
Proof.The inclusions (3.10) and (3.11) follow immediately from Lemma 3.4.Since F i , i = 1, . . ., p, are SCD mappings, O F i is dense in gph F i and from (3.10) we conclude that O F is dense in gph F. This proves that F is an SCD mapping.Now consider subspaces L i ∈ S F i (x i , y i ), i = 1, . . ., p.By taking into account the relation the inclusion (3.12) follows from (3.11) and (3.4).The statement about equality in (3.10) is again a consequence of Lemma 3.4 implying equality in (3.11) and (3.12).
The following large class of graphically Lipschitzian mappings covers many mappings important in applications; cf.[31], and is also important in the context of SCD mappings.Every multifunction F : R n ⇒ R n that is graphically Lischitzian of dimension n at some point (x, y) ∈ gph F, has the SCD property around (x, y) by [14,Proposition 3.17].We now state another property related to Proposition 3.5.Lemma 3.7.Let F : R n ⇒ R n be graphically Lipschitzian of dimension n at ( x, ȳ) ∈ gph F. Then there is an open neighborhood W of ( x, ȳ) such that for all (x, y) ∈ gph F ∩W the following properties hold: y) is a subspace, then the dimension of this subspace is n.
Proof.Regarding property (i), we refer to [14,Remark 3.18] and [32,Proposition 8.41].To show the second statement, let W , Φ, U and f be as in Definition 3.6 and consider (x, y) ∈ gph F ∩W such that T gph F (x, y) is a subspace.Denoting (u, f (u)) = Φ(x, y), we obtain that is a subspace with the same dimension as T gph F (x, y).Therefore, the graphical derivative D f (u, f (u)) is a linear mapping and f is Fréchet differentiable at u by [32,Exercise 9.25 Next, we turn to the notion of SCD regularity.Definition 3.8.
1. We denote by Z reg n the collection of all subspaces L ∈ Z n such that 2. A mapping F : R n ⇒ R n is called SCD regular around (x, y) ∈ gph F, if F has the SCD property around (x, y) and i.e., L ∈ Z reg n for all L ∈ S * F(x, y).Further, we will denote by the modulus of SCD regularity of F around (x, y).
Since the elements of S * F(x, y) are contained in gph D * F(x, y), it follows from the Mordukhovich criterion (2.2) that SCD regularity is weaker than metric regularity, and consequently SCD regularity is also weaker than strong metric regularity.
In the following proposition, we state some basic properties of subspaces L ∈ Z reg n .
Proposition 3.9 (cf.[14,Proposition 4.2]).We have L ∈ Z reg n if and only if for every (A, B) ∈ M (L) the matrix B is not singular.Thus, for every L ∈ Z Note that for every L ∈ Z reg n there is C L = AB −1 for all (A, B) ∈ M (L).Combining [14, Equation (34), Lemma 4.7 and Proposition 4.8] we obtain the following lemma.Lemma 3.10.Assume that F : R n ⇒ R n is SCD regular around ( x, ȳ) ∈ gph F. Then Moreover, F is SCD regular around every (x, y) ∈ gph F sufficiently close to ( x, ȳ) and lim sup

On semismooth* Newton methods for SCD mappings
In this section we recall the general framework for the semismooth * Newton method introduced in [13] and adapted to SCD mappings in [14].Consider the inclusion where F : R n ⇒ R n is a mapping having the SCD property around some point ( x, 0) ∈ gph F.
The following notion relaxes the semismooth * property from Definition 2.4.
Definition 4.1.We say that F : R n ⇒ R n is SCD semismooth * at ( x, ȳ) ∈ gph F if F has the SCD property around ( x, ȳ) and for every ε > 0 there is some δ > 0 such that the inequality holds for all (x, y) ∈ gph F ∩ B δ ( x, ȳ) and all (y * , x * ) belonging to any L ∈ S * F(x, y).
Clearly, every mapping with the SCD property around ( x, ȳ) ∈ gph F which is semismooth * at ( x, ȳ) is automatically SCD semismooth * at ( x, ȳ).Therefore, the class of SCD semismooth * mappings is even richer than the class of semismooth * maps.In particular, it follows from [22, Theorem 2] that every mapping whose graph is a closed subanalytic set is SCD semismooth * , cf. [14].
The following proposition provides the key estimate for the semismooth * Newton method for SCD mappings.
We now describe the SCD variant of the semismooth * Newton method.Given a solution x ∈ F −1 (0) of (4.14) and some positive scalar, we define the mappings A η, x : R n ⇒ R n × R n and N η, x : ).Assume that F is SCD semismooth * in ( x, 0) ∈ gph F and SCD regular around ( x, 0) and let η > 0. Then there is some δ > 0 such that for every x ∈ B δ ( x) the mapping F is SCD regular around every point ( x, ŷ) ∈ A η, x(x).Furthermore, for every ε > 0 there is some δ Assuming that we are given some iterate x (k) , the next iterate is given formally by x (k+1) ∈ N η, x(x (k) ).Let us take a closer look at this rule.Since we are dealing with set-valued mappings F, we cannot expect, in general, that F(x (k) ) = / 0 or that 0 is close to F(x (k) ), even if x (k) is close to a solution x.Therefore, we first perform some step that produces ( x(k) , ŷ(k) ) ∈ gph F as an approximate projection of (x (k) , 0) onto gph F. We require that for some constant η > 0, i.e., ( is true with some β ≥ 1, then and thus (4.15) holds with η = β + 1 and we can fulfill the inequality (4.15) without knowing the solution x.Furthermore, we require that S * F( x(k) , ŷ(k) ) ∩ Z reg n = / 0 and compute the new iterate as In fact, in our numerical implementation we will not compute the matrix C L , but two n × n matrices A, B such that L = rge (B T , A T ).The next iterate x (k+1) is then obtained by x (k+1) = x(k) + ∆x (k) , where ∆x (k) is a solution of the system A∆x = −B ŷ(k) .Alternatively, in view of Proposition 3.9, we can also choose a subspace L ∈ S F( x(k) , ŷ(k) ) ∩ Z reg n and compute the Newton direction as ∆x This leads to the following conceptual algorithm.
4. Newton step: Compute the Newton direction ∆x (k) by one of the following two alternatives: and calculate the Newton direction ∆x (k) as a solution of the linear system compute a solution p of the linear system and obtain the Newton direction as ∆x (k) = A (k) p.
For the choice between the two approaches to calculate the Newton direction, it is important to consider whether an element from For this algorithm, locally superlinear convergence follows from Proposition 4.3, see also [14,Corollary 5.6].
Theorem 4.4.Assume that F is SCD semismooth * at ( x, 0) ∈ gph F and SCD regular around ( x, 0).Then for every η > 0 there is a neighborhood U of x such that for every starting point x (0) ∈ U Algorithm 1 is well defined and stops after finitely many iterations at a solution of (4.14) or produces a sequence x (k) that superlinearly converges to x for any choice of ( x(k) , ŷ(k) ) satisfying (4.15) and any In particular, if F is strongly metrically regular around ( x, 0), then F has the SCD property around ( x, 0) by [14,Corollary 3.19] and it is also SCD regular around ( x, 0) as pointed out in the previous section.Therefore, if F also happens to be SCD semismooth * around ( x, 0), then the assumptions of the above statement are fulfilled.
Note that for an implementation of the Newton step, we need not know the whole derivative 5 On the implementation of the SCD semismooth* Newton method When trying to implement the SCD semismooth * Newton method directly for (1.1), it turns out that it can be rather difficult to perform the approximation step.Hence, we consider another equivalent approach which is more flexible.Consider an equivalent reformulation of (1.1) by the (decoupled) GE The new variable d acts only as an auxiliary variable.The approximation step now reads as follows: Given x (k) close to a solution x (and arbitrary d (k) , for example, d (k) = x (k) ), set x(k) := x (k) and find a point d(k) close to x (k) such that dist(0, f (x (k) ) + Q( d(k) )) is small.An approach to solving this problem could be to rewrite GE (1.1) in fixed point form x ∈ T (x) and select d(k) ∈ T (x (k) ).For example, for any λ > 0 there is (5.17) follows.In order to show that this approach is feasible as an approximation step, we have to verify that a bound of the form (4.15) holds, at least for x (k) close to x.
Proposition 5.1.Let x be a solution of (1.1) and assume that there is some λ > 0 such that the resolvent (I + λ Q) −1 has a single-valued Lipschitz continuous localization S around x − λ f ( x) for x, i.e., there are neighborhoods V of x − λ f ( x) and U of x such that S : V → U is Lipschitz continuous and Then there is some δ > 0 and some η > 0 such that for every x ∈ B δ ( x) there holds x − λ f (x) ∈ V and the vectors Proof.Choose δ > 0 such that B δ ( x) ⊆ (I − λ f ) −1 (V ) and let L f , L S > 0 denote the moduli of Lipschitz continuity of f on B δ ( x) and of S on V , respectively.Consider x ∈ B δ ( x).Then, by construction, x − λ f (x) ∈ V and hence d := S(x − λ f (x)) is well defined.Further, by (5.17), we have

In addition we have
and the assertion follows.
In particular, if Q is a maximal hypomonotone mapping, i.e., there exists some γ ≥ 0 such that γI + Q is maximal monotone, then for every 0 < λ < 1/γ the mapping (I + λ Q) is strongly monotone and hence (I + λ Q) −1 is a single-valued Lipschitz continuous function on R n , cf. [32, Proposition 12.54].However, hypomonotonicity is only a sufficient condition ensuring that (I + λ Q) −1 has this property.In Section 6 we will encounter a non-hypomonotone mapping Q, such that (I + λ Q) −1 is single-valued and Lipschitz continuous for every λ > 0.
Note that the choice d ∈ (I + λ Q) −1 x − λ f (x) corresponds to one step of the so-called Forward-Backward method for solving (1.1).
In the next proposition, we summarize some properties of G.
(i) For every x ∈ R n and (d, z) ∈ gph Q we have (iii) The mapping H is SCD regular around ( x, 0) if and only if G is SCD regular around ( x, x), (0, 0)).where P L corresponds to the orthogonal projection onto the subspace L := gph DQ(d, z).Hence, for every (d, z) ∈ gph Q and every x ∈ R n we obtain

Proof. G has the representation
and from Proposition 3.3 we conclude Similarly, we have yielding, together with Proposition 3.3, By virtue of (i) and the definition of SCD regularity, G is SCD regular around (x, d), ( f (x) + z, x − d) if and only if for every pair X,Y with rge (X,Y ) ∈ S Q(d, z) the matrix is nonsingular and the equivalence between (a) and (c) follows.
To establish (iii), just note that by Proposition 3.3 we have and the assertion follows from (ii) and the definition of SCD regularity.
Let us now consider the Newton step.Assume that, emanating from the iterate x (k) , we have computed 2 ) ∈ gph G as the result of the approximation step.Case (i): We compute the Newton direction according to step 4.a) of Algorithm 1.By Proposition 5.2, we have to compute two n × n matrices Y * (k) , X * (k) with rge )) and to solve the system Using the second equation we can eliminate ∆d (k) = ∆x (k) + ŷ(k) 2 and arrive at the linear system (5.19)

Case (ii):
The Newton direction is computed by step 4.b) of Algorithm 1.In this case we determine two n×n matrices By eliminating p 1 = X (k) p 2 − ŷ(k) 2 we obtain the linear system whose solution yields In both cases, the new iterate is given by x (k+1) := x(k) + ∆x (k) .Further, we have ∆x

Algebraic form of the discrete contact problem with Coulomb friction
We consider an elastic body represented by a bounded domain Ω ⊂ R 3 with a sufficiently smooth boundary ∂ Ω.The body is made of elastic, homogeneous, and isotropic material.The boundary consists of three non-empty disjoint parts: . Zero displacements are prescribed on Γ u , surface tractions act on Γ p , and the body is subject to volume forces.We seek a displacement field and a corresponding stress field satisfying the Lamé system of PDEs in Ω, the homogeneous Dirichlet boundary conditions on Γ u , and the Neumann boundary conditions on Γ p .The body is unilaterally supported along Γ c by some flat rigid foundation given by the halfspace R 2 × R − and the initial gap between the body and the rigid foundation is denoted by d(x), x ∈ Γ c .In the contact zone, we consider a static Coulomb Friction condition.This problem can be described by partial differential equations and boundary conditions for the displacements, which we are looking for.We refer the reader to, e.g., [11], where also a weak formulation can be found.We consider here only the discrete algebraic problem, which arises after some suitable finite element approximation.
Let n denote the number of degrees of freedom of the nodal displacement vector and let p denote the number of contact nodes x i ∈ Γ c \ Γ u .After some suitable reordering of the variables, such that the first 3p positions are occupied by the displacements of the nodes lying in the contact part of the boundary, we arrive at the following nodal block structure for an arbitrary vector y ∈ R n : In what follows, A ∈ R n×n , l ∈ R n are the stiffness matrix and the load vector, respectively.Further we are given two matrices N ∈ R p×n , T ∈ R 2p×n , where, for a given displacement vector v, Nv yields the normal components at the p contact points, and T v = (T 1 v, . . ., T p v), where T i v ∈ R 2 is the tangential nodal displacement vector at the i-th contact node.The symbol |T v| ∈ R p denotes a vector defined by We denote with α ∈ R p the vector of nodal distances with α i := d(x i ) and the friction coefficient is denoted by F .
Since the stiffness matrix A is positive definite and λ ≥ 0, condition (6.20) is equivalent to the requirement that ũ is a minimizer of the convex minimization problem Given a vector z = (z 1 , z 2 , z 3 ) T ∈ R 3 , we denote by z 12 := (z 1 , z 2 ) T ∈ R 2 the vector formed by the first two components.Using this notation, we have due to the ordering of the nodal displacements.
Next consider the transformation of variables u = ũ + d, where where l := l −Ad.Since the objective in this minimization problem is convex, this is in turn equivalent to the first-order optimality condition 0 ∈ (Au − l) i − λ i (0, 0, 1) Further, (6.21) is the same as ), i = 1, . . ., p.After eliminating λ from explicit variables, we have thus arrived at the GE where ) is dependent solely on transformed displacements u.Multipliers λ i appear only implicitly as −ϑ in the description of Q.This is a big difference with respect to other approaches, where the semismooth Newton method for equations is applied to mixed primal-dual systems or purely dual systems using some NCP-functions, see, e.g., [33], [5], [28].
Remark 6.2.We have derived the GE (6.22) for the Signorini problem with static Coulomb friction.We claim also that, for other contact problems with friction involving two elastic bodies, one can derive a GE of the same type.The interested reader is referred to [33, Section 5.2] for an algebraic transformation of a two-body problem to a one-body problem.
Note that which enables us to prove the following statement.
Proposition 6.3.H is semismooth * at each point in its graph.
Proof.Consider a point ( ū, w) ∈ gph H.By [13, Proposition 3.6] it suffices to show that Q is semismooth * at ( ū, w − A ū − l), which definitely holds true provided Q is semismooth * at all points of its graph.Thus, invoking [22, Theorem 3] and using the connection between the semismooth * property of sets and the respective distance functions, it suffices to show that gph Q is a subanalytic set.Let us pick a reference point ( v, ḡ, θ ) ∈ gph Q and consider the set Clearly, P is semianalytic (as the intersection of finitely many polynomial equalities and inequalities, it is even semialgebraic) and compact.Moreover, by construction, gph Q ∩ B 1 ( v, ḡ, θ ) is the canonical projection of P onto the space of variables (v, g, ϑ ) and hence subanalytic, cp.[4].The proof is complete.Proposition 6.4.For every γ > 0, the mapping (γI 3 + Q) −1 : R 3 ⇒ R 3 is single-valued and Lipschitz continuous on R 3 .In particular, Q is graphically Lipschitzian of dimension 3 at every point of its graph.
Proof.We have v ∈ (γI 3 + Q) −1 (w) if and only if for some ϑ ∈ N R + (v 3 ) and some v * 12 ∈ ∂ v 12 .Since γI 1 + N R + is both maximal monotone and strongly monotone, v 3 and ϑ are uniquely given by For given ϑ ≤ 0 the mapping γI 2 + F (−ϑ )∂ • is again maximal monotone and strongly monotone and thus v 12 is uniquely given by Table 1: Definitions and mechanic interpretations of the sets from (6.28).
Note that in Table 1 the impossible combinations of variables are crossed out.
Proposition 6.6.Q is an SCD mapping and for ( v, ḡ, θ and for ( v, ḡ, θ ) ∈ M 1 one has Proof.Since Q is graphically Lipschitzian of dimension 3 at every point of its graph, it is an SCD mapping.
Note that the sets L, M 1 and M + 3 exhibit a stable behavior in the sense that, for a sufficiently small ρ > 0, In particular, we have It follows from Definition 2.1 that In all three cases, we have therefore to do with linear subspaces of dimension three, which yield and the expressions for S Q( v, ḡ, θ ) in (6.29), (6.30), and (6.31).Concerning the expressions for S * Q( v, ḡ, θ ), they can be derived by first computing the respective orthogonal complements and then using the relation (3.3).The equalities S Q( v, ḡ, θ ) = S Q( v, ḡ, θ ) and S * Q( v, ḡ, θ ) = S * Q( v, ḡ, θ ) in (6.29)-(6.32)follow from the observation that the matrices that describe the subspaces continuously depend on the argument ( v, ḡ, θ ).
From Table 1 one can further infer that (i) every point from M 2 is accessible by sequences belonging solely to L or to M 1 ; (ii) every point from M − 3 is accessible by sequences belonging to M 1 or to M + 3 , and (iii) the singleton M 4 = {(0, 0, 0)} is accessible by sequences belonging to L or to M 1 or to M − 3 or to M + 3 .This implies in particular that (6.35) Formulas (6.35) are used in the implementation of the Newton step of the SCD semismooth * Newton method to the numerical solution of (6.22) in the next section.
Obviously, the mapping Q R given by (6.24) is Lipschitzian and, therefore, graphically Lipschitzian of dimension n − 3p as well.Further Combining Lemma 3.4, Proposition 3.5 and Lemma 3.7 with Proposition 6.4 and Proposition 6.6 we arrive at the following corollary.
and consequently also the mapping • γI n has these properties.Thus, for a given iterate u (k) , the choice is feasible for the approximation step by Proposition 5.1.

Numerical experiments
We treat various geometries arising from the cuboid (0, 2) × (0, 1) × (0, 1) ⊂ R 3 by modifying its bottom surface.Given a function d : (0, 2) × (0, 1) → R, we consider the elastic body represented by the domain At the left surface x 1 = 0 of the body, we impose Dirichlet boundary conditions and on the top surface x 3 = 1 and the right surface x 1 = 2 act surface tractions with densities P T and P R .The rigid obstacle is given by the half-space R 2 × R − so that the contact boundary is the bottom surface of the body.Depending on the discretization parameter lev, the elastic body is uniformly cut into n x 1 • n x 2 • n x 3 hexahedrons, where This results in a hexahedral mesh with (n x 1 + 1) • (n x 2 + 1) • (n x 3 + 1) vertices, where (n x 2 + 1) • (n x 3 + 1) vertices are in the Dirichlet part of the boundary and p := n x 1 • (n x 2 + 1) vertices belong to the contact part of the boundary.The total number of degrees of freedom of the nodal displacements is The setting is shown in Figure 1.
The resulting GE (6.22) is solved with the SCD semismooth * Newton method described in Sections 5 and 6 and the implementation was carried out in MATLAB on a PC with an i7-7700 CPU.The part of the code that describes the model is built on the original code of [3] with the accelerated assembly of the elastic stiffness matrix as described in [7].This part of the code was also applied to the Tresca friction solver of [16].The main difference between the implementations is that in the new code, no reduction is done to the nodes on the contact boundary, and the complete problem is treated with all domain nodes.
For the scalar γ used in the approximation step (6.36) we used an approximation of the largest eigenvalue of A obtained by five iterations of the power method.The system (5.19) that defines the Newton direction was solved iteratively using the GMRES method with ILU factorization as a preconditioner.We stopped the GMRES method when the relative residual (non-preconditioned) is less than the prescribed tolerance tol.Clearly, in this case, we will lose the superlinear convergence, but we can expect linear convergence with the rate tol.Of course, we can use more advanced methods to solve the linear system that determines the Newton direction, but the main task of this paper is to demonstrate the efficiency of the semismooth * Newton method and not the solution of linear systems.
As material parameters, we chose Young's modulus E := 70 GPa and Poisson's ratio ν = 0.334 (aluminum).The friction coefficient was always chosen as F = 0.23.In Table 2 we report for different discretization levels lev the number p of nodes in the contact part of the boundary and the number n of degrees of freedom.Furthermore, using the relative accuracy tol = 0.1 to calculate Newton's direction, for each of the six possible combinations of geometries d 1 , d 2 , d 3 and load cases L 1 , L 2 we list the number of Newton iterations it and the total number gmres of GMRES iterations needed to reduce the initial residual ŷ(0) by a factor of 10 −12 .The starting point u (0) for the semismooth * Newton method was always chosen as the origin.The value gmres characterizes the computational complexity.
We can see that for every geometry and every load case the iteration numbers are nearly equal and increase only slightly with finer discretizations.
In Figure 2 we depict for the different cases the undeformed bottom surface and the contact states for the deformed contact boundary.The number of iterations will decrease when a better starting point is available.In Table 3 we display the iteration numbers when we use as a starting point the solution of the previous discretization level interpolated to the finer mesh.
We observe that the number of GMRES iterations is reduced by 35 − 45% at the highest discretization level.A closer analysis shows that, as expected, we essentially avoid iterations to localize the solution.In fact, after 1-3 iterations, we have identified the correct state of all nodes in the contact part of the boundary, and the remaining iterations are only needed to reach the desired accuracy.Note that we use tol = 0.1 and therefore expect linear convergence with the rate 0.1.Since we want to reduce the residual by a factor of 10 −12 , we expect about 12 semismooth * Newton steps to achieve this goal.In most cases, we need fewer iterations: The reason is that the relative residual of the calculated direction is sometimes significantly less than tol.
Finally, we investigate the impact of the parameter tol on the performance of the semismooth * Newton method.Here, we consider only the load case L 2 and that the bottom surface is given by d 3 with the discretization level lev = 10.In Table 4 we report the iteration numbers it of the semismooth * Newton method and the total number gmres of GMRES iterations for the starting point u (0) = 0. We can see that for tol = 0.1 we need the most semismooth * Newton steps; however, the total number of GMRES iterations, which measures computational complexity, is the lowest.
We show the convergence of the semismooth * Newton method for the four values of tol in Figure 3.We see that during the first 5 or 6 iterations, when the semismooth * Newton method tries to localize the solution, the accuracy tol does not play any role in decreasing the residual, and we only need a lot of GMRES iterations to compute the search directions with higher accuracy.As soon as Table 3: Iteration numbers for tol = 0.1 and a starting point u (0) set as the interpolated solution of the previous discretization level.we are sufficiently close to the solution, the increased accuracy for computing the search direction also yields better convergence rates and consequently fewer iterations for the semismooth * Newton method.However, we also need more GMRES iterations to calculate the search direction, which defeats the advantage of a better convergence rate.

Conclusion
The paper shows the abilities of the SCD semismooth * Newton method to compute, apart from the variational inequalities of the first and second kind, also a more complicated class of equilibria which can be modeled as GEs with an SCD and semismooth * multi-valued part.This is documented by a large-scale highly complicated contact problem, where the efficiency of the method enables us, in contrast to most existing approaches, to solve the respective GE on the whole domain without the timeconsuming reduction to nodes lying on the contact boundary.We do hope that the SCD semismooth * Newton method will exhibit a similar performance also in some other mechanical problems having a similar structure as the considered contact problem with Coulomb friction.
The paper is dedicated to our friend A.L. Dontchev, for whom nonsmooth Newton methods definitely belonged to favorite research topics and who contributed to the development of this area in a remarkable way, cf., e.g., [10], [8].

Definition 3 . 6
(cf.[32, Definition 9.66]).A mapping F : R n ⇒ R m is graphically Lipschitzian of dimension d at ( x, ȳ) ∈ gph F if there is an open neighborhood W of ( x, ȳ) and a one-to-one mapping Φ from W onto an open subset of R n+m with Φ and Φ −1 continuously differentiable, such that Φ(gph F ∩W ) is the graph of a Lipschitz continuous mapping f : U → R n+m−d , where U is an open set in R d .

(
ii) Let x ∈ R n and assume that Q has the SCD property around (d, z) ∈ gph Q.Then the following statements are equivalent: (a) G is SCD regular around (x, d), ( f (x) + z, x − d) .(b) For every L ∈ S * Q(d, z) and every (X,Y ) ∈ M (L) the matrix ∇ f (x)X +Y is nonsingular.(c) For every L ∈ S * Q(d, z) and every (Y * , X * ) ∈ M (L) the matrix ∇ f (x) T Y * + X * is nonsingular.

Definition 6 . 1 ([ 3 ,
Definition 3.6]).As a solution of a discrete contact problem with Coulomb friction we declare any couple

Figure 1 :
Figure 1: An undeformed and deformed cuboid domain are shown in the left and middle pictures.The zero Dirichlet boundary condition for displacements is assumed on the blue part of the boundary Γ u , surface tractions are applied to the green part of the boundary Γ p and the contact boundary Γ c is pressed against the (red) rigid plane foundation.Front faces are not visualized.The right picture shows the deformed cuboid domain decomposed in hexahedrons.

Figure 2 :
Figure 2: Undeformed bottom surface and states in the deformed contact boundary: no contact (blue), sliding (yellow), sticking (red).

Table 2 :
Iteration numbers for tol = 0.1 and a starting point u (0) = 0. point (x(k), d(k) ).In the context of our contact problem with Coulomb friction, given an iterate u (k) of nodal displacements and d(k) by (6.36), we consider ŷ