S H Üeber

G Stadler

B I Wohlmuth
email: wohlmuth@ians.uni-stuttgart.de

A PRIMAL-DUAL ACTIVE SET ALGORITHM FOR THREE-DIMENSIONAL CONTACT PROBLEMS WITH COULOMB FRICTION *

Keywords: contact problems, dual Lagrange multipliers, inexact primaldual active set strategy, semismooth Newton methods, nonlinear multigrid

In this paper, efficient algorithms for contact problems with Tresca and Coulomb friction in three dimensions are presented and analyzed. The numerical approximation is based on mortar methods for nonconforming meshes with dual Lagrange multipliers. Using a nonsmooth complementarity function for the three-dimensional friction conditions, a primal-dual active set algorithm is derived. The method determines active contact and friction nodes and, at the same time, resolves the additional nonlinearity originating from sliding nodes. No regularization and no penalization are applied, and superlinear convergence can be observed locally. In combination with a multigrid method, it defines a robust and fast strategy for contact problems with Tresca or Coulomb friction. The efficiency and flexibility of the method is illustrated by several numerical examples.

Key words. 3D Coulomb friction,

as intersections of rotated squares in order to make the application of optimization algorithms possible. Improvements are proposed in [START_REF] Kučera | A new FETI-based algorithm for solving 3D contact problems with Coulomb friction[END_REF][START_REF] Kučera | Minimizing quadratic functions with separable quadratic constraints[END_REF]. A different idea is followed in [START_REF] Fackeldey | Solving frictional contact problems with multigrid efficiency[END_REF][START_REF] Kornhuber | Adaptive multigrid methods for Signorini's problem in linear elasticity[END_REF], where monotone multigrid methods are used to construct a globally convergent solver. The implementation of these methods relies on a multilevel hierarchy of spaces and requires the use of modified coarse grid basis functions and suitable coarse grid constraints.

Already in an early paper [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF], Newton-type methods for contact problems with friction are used. Similar methods are also studied in the more recent contributions [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][START_REF] Christensen | Reformulation: Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods[END_REF], where the performance of generalized Newton-type methods for frictional contact problems is shown to be superior to interior point methods. The methods presented in [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][START_REF] Christensen | Reformulation: Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods[END_REF] rely on the reformulation of the contact and friction conditions using nonsmooth equations and on generalized differentiability concepts. Here we apply these strategies but use a different complementarity function. Additionally, we study and exploit the structures arising in Newton-type steps and relate them to primal-dual active set strategies. We also apply our algorithm to two-body contact problems which are discretized in terms of mortar techniques and introduce different stabilizations.

Recently, the application of mortar methods as a discretization for the continuous problem has gained a considerable amount of interest (see, e.g., [START_REF] Belgacem | Numerical simulation of some variational inequalities arisen from unilateral contact problems by the finite element methods[END_REF][START_REF] Fischer | Frictionless 2D Contact formulations for finite deformations based on the mortar method[END_REF][START_REF] Fischer | Mortar based frictional contact formulation for higher order interpolations using the moving friction cone[END_REF][START_REF] Hild | Quadratic finite element methods for unilateral contact problems[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF][START_REF] Laursen | Computational Contact and Impact Mechanics[END_REF][START_REF] Puso | A mortar segment-to-segment frictional contact method for large deformations[END_REF]). By using a dual basis function for the Lagrange multiplier, the weak form of the constraints can be written as independent constraints for each node on the slave side, resulting in a quasi-variational inequality. Using a suitable local basis transformation, only the coefficients on the slave side are affected by the nonpenetration condition and the friction constraints. Thus, the algebraic structure of a two-body contact problem is exactly the same as for a one-body contact problem. For details we refer to [START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF][START_REF] Wohlmuth | Monotone multigrid methods on nonmatching grids for nonlinear multibody contact problems[END_REF]. In the mixed formulation which is motivated by the mortar approach, the Lagrange multipliers are treated as independent variables. Thus primal-dual active set strategies which use both types of variables are of special interest. We refer to [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF], where the interpretation of the active set strategy as a semismooth Newton method is given. These methods have been successfully extended and applied to contact problems without friction [START_REF] Hintermüller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF] and in two dimensions to friction [START_REF] Hüeber | Efficient algorithms for problems with friction[END_REF][START_REF] Kunisch | Generalized Newton methods for the 2D-Signorini contact problem with friction in function space[END_REF][START_REF] Stadler | Semismooth Newton and augmented Lagrangian methods for a simplified friction problem[END_REF]. However, the application to friction problems in three dimensions is not at all straightforward. In this paper, we derive a primal-dual active set method by applying a semismooth Newton method to a nonlinear complementarity function that expresses the three-dimensional (3D) Tresca friction conditions. In each step of the algorithm, we solve a linear problem with Robin boundary conditions on a suitable subset of the friction boundary and with Dirichlet conditions on the complement. This linear system can be solved using any efficient iterative or fast direct solver. Here we apply a multigrid method as an iterative solver. By solving these systems only approximately, we obtain an inexact primal-dual active set strategy. The resulting algorithm can be regarded as a nonlinear multigrid method and yields a reliable and fast convergence. Finally, we present and study a full Newton approach for the Coulomb friction problem. Our implementation is based on the finite element toolbox UG (see [START_REF] Bastian | UG -A flexible software toolbox for solving partial differential equations[END_REF]), and the direct solver PARDISO (see [START_REF] Schenk | Efficient sparse LU factorization with left-right looking strategy on shared memory multiprocessors[END_REF]) is used for the full Newton approach.

The outline of this paper is as follows: The two-body Coulomb frictional contact problem is stated in section 2, and the discretization based on mortar methods is briefly discussed. In section 3, we develop and analyze our algorithm for 3D friction problems. For simplicity of presentation, we restrict ourselves to the case of a onebody contact problem with Tresca friction. In section 4, we study the performance of the algorithm for a test problem with Tresca friction. Using fixed point ideas, this approach is extended to contact problems with Coulomb friction in section 5. In section 6, we apply our algorithms to two-body contact problems with Coulomb friction in three dimensions. Finally in section 7, we present a full Newton approach for Coulomb friction in three dimensions and compare its performance with the discussed fixed point-based methods.

Contact problem with Coulomb friction in three dimensions.

In this section, we formulate the contact problem with Coulomb friction in linear elasticity. In addition, we briefly describe our nonconforming discretization based on mortar methods and the applied basis transformation. It guarantees that the algebraic structure of the two-body system is exactly the same as in the case of a one-body contact problem.

Problem statement.

Let Ω l ⊂ R 3 , l ∈ {m, s}, denote two elastic bodies. The superscript "l" will refer to s for the slave or to m for the master body, as is common in the mortar setting. We assume that the boundaries ∂Ω l are divided into three disjoint measurable parts Γ l D , Γ l N , and Γ l C with meas(Γ l D ) = 0. We impose homogeneous Dirichlet conditions on Γ l D and Neumann data, i.e., a surface traction

p l ∈ (L 2 (Γ l N )) 3 on Γ l N .
Moreover, we denote by f l ∈ (L 2 (Ω l )) 3 the volume forces acting on Ω l . We use the constitutive law for linear elasticity, namely,

σ l := λ l tr(ε l )Id + 2μ l ε l in Ω l ,
where σ l denotes the stress field and ε l := 1/2(∇u l + ∇u l ) the linearized strain tensor that both depend on the displacement u l . The Lamé constants λ l , μ l > 0 are given by μ l := E l /(2(1 + ν l )) and λ l := (E l ν l )/((1 + ν l )(1 -2ν l )) with Young's modulus E l > 0 and Poisson ratio ν l ∈ (0, 0.5). Moreover, tr denotes the matrix trace operator and Id the identity matrix. Then the displacement u l satisfies the following equations:

Divσ l + f l = 0 in Ω l , (2.1a) u l = 0 on Γ l D , (2.1b) σ l n l -p l = 0 on Γ l N , (2.1c)
where n l denotes the unit normal outward vector on the boundary ∂Ω l . To state the contact and friction conditions, we introduce for each point of Γ s C the vectors τ 1 , τ 2 that span the tangential plane and use n := n s . We assume that {n, τ 1 , τ 2 } is an orthonormal basis in R 3 for each point of Γ s C . In order to formulate the nonpenetration condition of the two bodies, we use a predefined relation between the points of the possible contact zones Γ l C . This relation is realized by a smooth mapping

R : Γ s C → Γ m C satisfying R(Γ s C ) ⊂ Γ m C .
We assume that the mapping R is well-defined and maps any x ∈ Γ s C to the intersection of the normal on Γ s C at x with Γ m C . Then the contact conditions on Γ s C are given by (2.2)

[u] n -d ≤ 0, σ n ≤ 0, ([u] n -d)σ n = 0 on Γ s C ,
σ n = n σ s (x)n = n σ m (R(x))n on Γ s C .
Additionally, we have to ensure

τ k σ s (x)τ k = τ k σ m (R(x))τ k , k = 1, 2.
Finally, the Coulomb friction law states that (2.4)

⎧ ⎪ ⎨ ⎪ ⎩ σ τ ≤ F|σ n |, σ τ < F|σ n | ⇒ [u] τ = 0, σ τ = F|σ n | ⇒ ∃β ≥ 0 : σ τ = -β[u] τ , on Γ s C .
Here • denotes the Euclidean norm in R 3 , σ τ := σ s nσ n n the tangential component of the boundary stress, and [u] τ the relative tangential displacement given by

[u] τ := [u] -[u] n n. Furthermore, F : Γ s C → R, F ≥ 0,
is the friction coefficient. 2.2. Mortar discretization. By using dual Lagrange multipliers, we can apply locally a suitable basis transformation for the finite element basis. The shape functions on the master side are modified by adding a linear combination of nodal shape functions on the slave side. In the new basis the shape functions on the master side satisfy a weak continuity condition; i.e., the jump of these basis functions tested with a Lagrange multiplier being defined on the slave side is zero. Thus, the coefficients in this new basis on the slave side of the contact zone describe the relative displacement between the contact interfaces. In this basis, the two-body contact problem has the same structure as in the one-body case, since all constraints at the contact zone are restricted to the degrees of freedom on the slave side. For details, see [START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF][START_REF] Wohlmuth | Monotone multigrid methods on nonmatching grids for nonlinear multibody contact problems[END_REF]. Motivated by these considerations, in what follows we will write u h for the coefficient vector with respect to the new constrained basis. Therefore, the following discussion holds for both the one-body and the two-body cases. We denote the multiplier corresponding to the discretization of -σ s n on the slave side by λ h . In the new basis, (2.1) is satisfied in the discrete form (2.5)

A h u h + B h λ h = f h ,
with B h of the form (0, D) , where due to the use of dual Lagrange multipliers the matrix D is diagonal. Let us denote by S the set of all nodes of the finite element mesh belonging to Γ s C and by N the set containing all other nodes. Then, for each p ∈ S, the entry of the matrix D is given by

(2.6) D[p, q] = Id 3 Γ C φ p ψ q ds = Id 3 Γ C φ p ds δ pq =: δ pq D p Id 3 ,
where Id 3 denotes the identity matrix in R 3×3 and φ p and ψ q the primal and dual basis function associated with the node p, respectively. We remark that in (2.6) the scaled biorthogonality of the sets {φ p } p∈S and {ψ p } p∈S is exploited. If the displacements u h are known, the Lagrange multiplier can be computed directly from (2.5) as

(2.7)

λ h = D -1 (f h -A h u h ) S ,
where the subscript S on the right side indicates that we use only the entries of the vector corresponding to the nodes p ∈ S.

We next introduce the scaled normal and tangential components of the multiplier λ h and the normal and tangential components of the relative deformation u h . For a node p ∈ S, let

λ τ,p,s := λ p D p τ 1p λ p D p τ 2p ∈ R 2 and λ n,p,s := λ p D p n p ∈ R, u τ,p := u p τ 1p u p τ 2p ∈ R 2 and u n,p := u p n p ∈ R.
The discretized (and scaled) gap at the node p ∈ S is defined as

d p := 1 D p Γ s C d ψ p ds.
Then the discrete conditions for the normal contact (2.2) are given by (2.8)

u n,p ≤ d p , λ n,p,s ≥ 0, λ n,p,s (u n,p -d p ) = 0, p ∈ S.
We note that we use different scaling factors for the primal and dual variables and that D p is proportional to the local mesh size. This is motivated by the fact that the H -1/2 -norm for the Lagrange multiplier and the H 1/2 -norm for the displacements have the same error reduction. We remark that the proposed scaling factors yield better numerical convergence rates for the inexact version of the algorithms. The discrete Coulomb friction conditions are given by (2.9)

⎧ ⎪ ⎨ ⎪ ⎩ λ τ,p,s ≤ F|λ n,p,s |, λ τ,p,s < F|λ n,p,s | ⇒ u τ,p = 0, λ τ,p,s = F|λ n,p,s | ⇒ ∃β ≥ 0 : λ τ,p,s = βu τ,p ,
for all p ∈ S.

Here, for simplicity, F is assumed to be constant and independent of the solution.

Differently from the Coulomb friction conditions (2.4), where the friction bound is F|σ n |, for Tresca friction this bound is a given function g ∈ H -1/2 (Γ s C ), with g ≥ 0. The corresponding discrete bound is

g p := Γ s C g φ p ds;
i.e., for Tresca friction, we replace the friction bound F|λ n,p,s | in (2.9) by g p .

Tresca friction in three dimensions.

In this section, we introduce our algorithm for Tresca friction in the 3D case. To simplify the presentation, we restrict ourselves to the case u n,p = 0 for all nodes p ∈ S. In our iterative algorithm, we have to solve in each iteration step a linear problem with boundary conditions of Dirichlet-, Neumann-, or Robin-type on the contact zone Γ C . The detection of the zones for the different types of boundary conditions is based on a primal-dual approach. We recall that the Lagrange multiplier playing the role of the dual variable can be locally computed in a postprocess from the primal variable u h ; see (2.7). This algorithm for 3D friction is closely related to the primal-dual active set method for 2D friction. Consequently, it perfectly fits into the abstract framework used in [START_REF] Hüeber | Efficient algorithms for problems with friction[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF] and inherits the advantages of these methods, i.e., their simple implementation and their numerical efficiency. The method can be regarded as an active set strategy or alternatively as a semismooth Newton method guaranteeing fast local convergence; see Theorem 4.1.

Iterative solver.

To start with, we review the conditions for 3D Tresca friction as given by (2.9). Since for g p = 0 these conditions simplify to homogeneous Neumann boundary conditions in the tangential direction, we first assume that g p > 0 and later comment on the case g p = 0. By a straightforward calculation, it can be verified that (2.9) with g p instead of F|λ n,p,s | is equivalent to C(u τ,p , λ τ,p,s ) = 0 for all p ∈ S, c τ > 0, where the nonlinear complementarity function C(• , •) is defined by

(3.1) C(u τ,p , λ τ,p,s ) := max(g p , λ τ,p,s + c τ u τ,p )λ τ,p,s -g p (λ τ,p,s + c τ u τ,p ) .
In the following, we take (3.1) as the starting point for our algorithm. Our main idea is to apply a Newton-type algorithm for the solution of C(u τ,p , λ τ,p,s ) = 0. Unfortunately, both the Euclidean norm and the max-function are not smooth and not differentiable in the classical sense. However, they are semismooth in the sense of [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF] which justifies the application of a semismooth Newton method. As generalized derivative for f (b

) := max(a, b), we use G f (b) = 0 if a ≥ b and G f (b) = 1 if a < b.
We remark that the arbitrary choice for the case a = b does not influence the convergence rate of the proposed algorithm. Note that, in the first term of (3.1), the Euclidean norm appears for nonzero arguments only. This is due to the fact that if λ τ,p,s + c τ u τ,p = 0, we obtain max(g p , λ τ,p,s + c τ u τ,p ) = g p , and the Euclidean norm vanishes. Therefore, the only nondifferentiability that matters in (3.1) is the max-function. In the semismooth Newton step the derivative of the Euclidean norm occurs only for points that are differentiable in the classical sense.

We now compute, for p ∈ S, the generalized derivative

G C of C(• , •). For the variation (δu τ,p , δλ τ,p,s ) ∈ R 2 × R 2 ,
I k τ := {p ∈ S : λ k-1 τ,p,s + c τ u k-1 τ,p -g p ≤ 0}, (3.4a) A k τ := {p ∈ S : λ k-1 τ,p,s + c τ u k-1 τ,p -g p > 0}. (3.4b)
Using this notation and (3.2) in (3.3), a straightforward computation shows that the new iterate (u k τ,p , λ k τ,p,s ) satisfies

u k τ,p = 0 for p ∈ I k τ , (3.5a) Id 2 -M k-1 p λ k τ,p,s -c τ M k-1 p u k τ,p = h k-1 p for p ∈ A k τ , (3.5b) where M k-1 p := e k-1 p (Id 2 -F k-1 p
), with the scalar value e k-1 p and the 2 × 2-matrix F k-1 p given by, respectively,

(3.6) e k-1 p := g p λ k-1 τ,p,s + c τ u k-1 τ,p , F k-1 p := λ k-1 τ,p,s λ k-1 τ,p,s + c τ u k-1 τ,p g p λ k-1 τ,p,s + c τ u k-1 τ,p

, and the vector h

k-1 p ∈ R 2 (3.7) h k-1 p := e k-1 p F k-1 p (λ k-1 τ,p,s + c τ u k-1 τ,p ) = λ k-1 p .
Note that, while on the inactive set I k τ Dirichlet conditions are imposed, the condition on the active set A k τ is of Robin type since it involves the displacement u k τ,p and the surface traction λ k τ,p,s . The sets A k τ and I k τ approximate the sets of slippy and of sticky nodes, respectively. In the more general setting of Signorini conditions, we have also to set the boundary condition in the normal direction. The Robin condition (3.5b) can be easily handled if we rewrite (3.5b) as

(3.8) -λ k τ,p,s + L k-1 p u k τ,p = r k-1 p .
We note that r k-1 p enters in the right-hand side of the linear system and L k-1 p gives a contribution to the system matrix. Comparing (3.8) with (3.5b) and (3.7), we get, under the assumption that

Id 2 -M k-1 p is regular, L k-1 p := c τ (Id 2 -M k-1 p ) -1 M k-1 p = c τ (Id 2 -M k-1 p ) -1 -Id 2 , (3.9a) r k-1 p := -(Id 2 -M k-1 p ) -1 h k-1 p = γ k-1 p λ k-1 p , (3.9b)
where γ k-1 p is a suitable scaling factor. We mention that this Robin condition only guarantees positive definiteness of the system matrix if L k-1 p is positive definite. Note that for L k-1 p = 0 we find a pure Neumann condition. The degeneration of Robin-type to Dirichlet-type boundary conditions is not included in the form (3.9). However, this is not required in our situation, since nodes p with a Dirichlet condition belong to the set A k τ and therefore are not handled by (3.8). One can easily see that defining the iteration process Id 2 -M k-1 p is not necessarily regular and therefore not invertible. However, in the case of convergence, Id 2 -M k-1 p tends to a positive definite and symmetric matrix. This observation motivates the introduction of three possible modifications of the Robin system (3.5b) such that a regular matrix Id 2 -M k-1 p is obtained. Two of these modifications give a positive definite and symmetric matrix L k-1 p . We remark that all modifications converge in the limit case to the original system (3.5b) and thus preserve the local convergence properties of the algorithm.

Modifications of the Robin system.

To obtain a robust and convergent scheme, we have to replace the matrix F k-1 p by a scaled matrix F k-1 p,l , l = 1, 2, 3. The index l stands for one of the three possibilities considered in this paper. Our numerical results show that the scaling is essential for the robustness of the iteration scheme. According to the definition of M k-1 p and (3.7), we replace in

(3.5b) M k-1 p by M k-1 p,l and h k-1 p by hk-1 p,l given by M k-1 p,l := e k-1 p (Id 2 -F k-1 p,l ), hk-1 p,l := e k-1 p F k-1 p,l (λ k-1 τ,p,s + c τ u k-1 τ,p ), l = 1, 2, 3.
In a second step, we replace the matrix

Id 2 -M k-1 p by Id 2 -β k-1 p,l
M k-1 p,l with a scaling factor β k-1 p,l > 0, such that the resulting matrix is regular. Then the form (3.8) of the Robin boundary conditions is written as

-λ k τ,p,s + Lk-1 p,l u k τ,p = rk-1 p,l , with (3.10) Lk-1 p,l := c τ (Id 2 -β k-1 p,l M k-1 p,l ) -1 -Id 2 , rk-1 p,l := -(Id 2 -β k-1 p,l M k-1 p,l ) -1 hk-1 p,l
for l = 1, 2, 3. We mention that similar modifications as done for the matrix F k-1 p are used in primal-dual algorithms for the minimization of functionals involving Euclidean norms; see, e.g., [START_REF] Andersen | An efficient primal-dual interior-point method for minimizing a sum of Euclidean norms[END_REF][START_REF] Chan | A nonlinear primal-dual method for total variationbased image restoration[END_REF][START_REF] Hintermüller | An infeasible primal-dual algorithm for total variationbased inf-convolution-type image restoration[END_REF].

Since for all modifications F

k-1 p,l → F k-1 p and β k-1 p,l → 1 as (u k τ,p , λ k τ,p,s
) converges to the solution, the modifications do not degrade the local superlinear convergence. Next we present the three possible modifications used in this paper.

First modification. We use a parameter β

k-1 p,1 = 1 in Id 2 -β k-1 p,1 M p in (3.10) only if (λ k-1 τ,p,s ) (λ k-1 τ,p,s + c τ u k-1 τ,p ) < 0.
This condition is equivalent to the fact that the angle between the two vectors is greater than 90 degrees. Since in the limit case both vectors are parallel, the modification applies only when the iterates are far away from the solution. We define the scaled matrix

(3.11) F k-1 p,1 := λ k-1 τ,p,s λ k-1 τ,p,s + c τ u k-1 τ,p max g p , λ k-1 τ,p,s λ k-1 τ,p,s + c τ u k-1 τ,p . Note that F k-1 p,1 only differs from F k-1 p if λ k-1 τ,p,s > g p , i.e.
, if the Lagrange multiplier is not in the feasible set given in (2.9).

α k-1 p := (λ k-1 τ,p,s ) λ k-1 τ,p,s + c τ u τ,p λ k-1 τ,p,s λ k-1 τ,p,s + c τ u τ,p , δ k-1 p := min λ k-1 τ,p,s g p , 1 ,
and it is easy to see that γ 1

Fp,1

= α k-1 p δ k-1 p is an eigenvalue of F k-1 p,1 with λ k-1 τ,p,s
as an eigenvector. The second eigenvalue is given by γ 2

Fp,1 = 0, where the corresponding eigenvector is the vector which is orthogonal to

λ k-1 τ,p,s + c τ u k-1 τ,p . Therefore, the eigenvalues of M k-1 p,1 are γ 1 Mp,1 = e k-1 p (1 -α k-1 p δ k-1 p ) and γ 2 Mp,1 = e k-1 p . Since -1 ≤ α k-1 p ≤ 1 and 0 ≤ δ k-1 p ≤ 1, we get due to 0 < e k-1 p < 1 the relation 0 ≤ γ 1 Mp,1 < 2 and 0 < γ 2 Mp,1 < 1. Using (3.10) with β k-1 p,1 := 1 1-α k-1 p δ k-1 p if α k-1 p < 0,
1 otherwise yields a unsymmetric matrix Lk-1 p,1 with positive eigenvalues. Second modification. In contrast to the first modification, we use a symmetric matrix F k-1 p,2 . Here we need a parameter β k-1 p,2 = 1 in (3.10) for all cases with

(λ k-1 τ,p,s ) (λ k-1 τ,p,s + c τ u k-1 τ,p ) = 0. We replace F k-1 p
by the symmetrization of (3.11), namely,

(3.12) F k-1 p,2 := λ k-1 τ,p,s λ k-1 τ,p,s + c τ u k-1 τ,p + λ k-1 τ,p,s + c τ u k-1 τ,p λ k-1 τ,p,s 2 max g p , λ k-1 τ,p,s λ k-1 τ,p,s + c τ u k-1 τ,p . One can prove that the eigenvalues of F k-1 p,2 are γ 1,2 Fp,2 = 1 2 (α k-1 p ±1)δ k-1 p ∈ [-1, 1],
and therefore the eigenvalues of the matrix

M k-1 p,2 are γ 1,2 Mp,2 = e k-1 p (2 -(α k-1 p ± 1)δ k-1 p )/2.
Using the same arguments as before, we get 0 ≤ γ 1,2 Mp,2 < 2. Setting

β k-1 p,2 := 2 2 -(α k-1 p -1)δ k-1 p
results in a symmetric and positive definite matrix Lk-1 p,2 .

Third modification. In the third modification, we use the matrix

(3.13) F k-1 p,3 := λ k-1 τ,p,s λ k-1 τ,p,s max g p , λ k-1 τ,p,s 2 instead of F k-1 p
. Obviously this matrix is symmetric and positive semidefinite with eigenvalues γ 1 defined by (3.10) with β k-1 p,3 = 1 is symmetric and positive definite. We remark that the matrix F k-1 p,3 converges in the limit case to the matrix F k-1 p since, in the solution for a node p ∈ A τ , we have

λ τ,p,s + c τ u τ,p λ τ,p,s + c τ u τ,p = λ τ,p,s g p .
Furthermore as for the unmodified case (3.7) and (3.9b), we observe that h k-1 p and r k-1 p are proportional to λ k-1 p for the first and the third modification.

3.3.

Extension to the case g p = 0. In the above discussion, we assumed g p > 0 for all p ∈ S. The reason for this assumption is that for g p = 0 one cannot stringently deduce λ τ,p,s = 0 (which follows directly from (2.9)) from C(u τ,p , λ τ,p,s ) = 0. However, if a Tresca friction combined with fixed point ideas is used to solve Coulomb friction problems, g p = 0 naturally occurs for all noncontact points, which makes the case g p = 0 rather important. Fortunately nodes p with g p = 0 can also be handled using (3.5) within the setting (3.11)-(3.13).

In the following, we consider p ∈ S with g p = 0. First, we assume = 0. In particular, these matrices are well-defined for this case. So only for the cases λ k-1 τ,p,s + c τ u k-1 τ,p = 0 or λ k-1 τ,p,s = 0 do we have to enforce the node p to be in A k-1 τ .

λ k-1 τ,p,s + c τ u k-1 τ,p > 0. Then p ∈ A k τ , since g p = 0. In the case λ k-1 τ,p,s + c τ u k-1 τ,p

General remarks.

Computing generalized derivatives of nonsmooth functionals is a delicate issue. While in [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF] mainly intuitive arguments are used, the related papers [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][START_REF] Christensen | Reformulation: Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods[END_REF] use the concept of Bouligand differentiability. This concept allows the use of globalization (e.g., linesearch) strategies, but calculating the search direction requires one to solve a nonlinear system in each Newton step. In these papers, this problem is circumvented by substituting this nonlinear system by a linear one. Despite this heuristic step, the authors report on good numerical results. The concept of semismoothness [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF] used in the present paper has the advantage that the search direction can be found by solving a linear system. Nevertheless, one can also prove local superlinear convergence of the iterates; see Theorem 4.1.

Furthermore, we remark that there is some freedom in choosing the nonlinear complementarity function to express the complementarity conditions for the Tresca friction law. For g p > 0, we can also work with (3.1) replaced by

(3.14) C(u τ,p , λ τ,p,s ) := λ τ,p,s -g p (λ τ,p,s + c τ u τ,p ) max(g p , λ τ,p,s + c τ u τ,p )
.

Complementarity functions closely related to (3.14) for dealing with friction conditions have been used in [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][START_REF] Christensen | Reformulation: Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods[END_REF]. A semismooth Newton iteration for the solution of C(u τ,p , λ τ,p,s ) = 0 results in an iteration rule that also uses the active and inactive sets defined in (3.4) but results in a modified iteration step on A k τ . Our numerical experience yields that algorithms based on (3.1) perform more robustly compared to those based on (3.14).

Algebraic representation.

Now we give the matrix representation of the algebraic system we have to solve in each iteration step. As mentioned above, we restrict ourselves to the case u n,p = 0. The tangential conditions are either Dirichlet conditions for the inactive nodes p ∈ I k τ or Robin conditions for the active nodes p ∈ A k τ . To rotate the stiffness matrix A arising from standard linear elasticity and introduced in (2.5), we define the matrix with the normal vectors by

N := ⎛ ⎜ ⎜ ⎝ . . . 0 0 0 0 • • • 0 n p 0 • • • 0 0 0 0 . . . ⎞ ⎟ ⎟ ⎠ ∈ R |S|×3|S|
and the matrix with tangential vectors by

T := ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ . . . 0 0 0 0 • • • 0 τ 1 p 0 • • • 0 0 • • • 0 τ 2 p 0 • • • 0 0 0 0 . . . ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ∈ R 2|S|×3|S| .
As above, for each p ∈ S, n p , τ 1p , and τ 2p denote the unit outward normal and the tangential vectors, respectively. Moreover, for p ∈ A k τ the matrices Lk-1 p,l are assembled into

L A k τ := diag Lk-1 p,l p∈A k τ ∈ R 2|A k τ |×2|A k τ | .
Similarly, we define the vector r A k τ by the entries rk-1 p,l for p ∈ A k τ . We now use the decomposition of S into A k τ and I k τ . The vector u h can then be partitioned into

(u N , u I k τ , u A k τ ) , the multiplier λ h into (λ I k τ , λ A k τ )
, and the right-hand side

f h into (f N , f I k τ , f A k τ )
. Correspondingly, we can decompose the matrices A h , N , and T . Then the arising linear system has the form

(3.15) ⎛ ⎜ ⎜ ⎝ A NN A NI k τ A NA k τ 0 Id I k τ 0 0 0 N A k τ T A k τ A A k τ N T A k τ A A k τ I k τ T A k τ A A k τ A k τ + L A k τ T A k τ ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ u k N u k I k τ u k A k τ ⎞ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎝ f N 0 0 T A k τ f A k τ + r A k τ ⎞ ⎟ ⎟ ⎠ .
We remark that the Dirichlet boundary condition (3.5a) on I k τ is reflected in the second row and the Robin-type conditions for p ∈ A k τ (see (3.5b) or (3.8)) are included in the fourth row. In the case of the more general Signorini condition also used in the next section, we have to replace the second and third rows by more complex ones; see [START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF].

Algorithm and numerical examples for Tresca friction.

In this section, we present an efficient primal-dual algorithm based on a multigrid method to solve the resulting linear problems. Furthermore, we study numerically the performance of our algorithm and compare the three modifications of section 3.2. While in the previous section, we have derived the iteration rule for Tresca friction without taking into account Signorini contact, the applied algorithm can handle Tresca friction and Signorini contact (2.2). The discrete Signorini conditions (2.8) are realized in terms of the nonlinear complementarity function

A k n := {p ∈ S : λ k-1 n,p,s + c n (u k-1 n,p -d p ) > 0}, I k n := {p ∈ S : λ k-1 n,p,s + c n (u k-1 n,p -d p ) ≤ 0}.
We next state the inexact primal-dual active set (IPDAS) strategy for contact with Tresca friction. Initialize u 0,0 h and λ 0 h as an initial solution.

(1) Define the active and inactive sets by

A k n := p ∈ S : λ k-1 n,p,s + c n u k-1,m n,p -d p > 0 , I k n := p ∈ S : λ k-1 n,p,s + c n u k-1,m n,p -d p ≤ 0 , A k τ := p ∈ S : λ k-1 τ,p,s + c τ u k-1,m τ,p -g p > 0 , I k τ := p ∈ S : λ k-1 τ,p,s + c τ u k-1,m τ,p -g p ≤ 0 .
(2) For i = 1, . . . , m, compute

u k,i h = M G(u k,i-1 h , A k n , I k n , A k τ , I k τ , u k-1,m h , λ k-1 h ). (3) If u k,m h -u k,0 h / u k,m h < ε u , stop.
(4) Compute the Lagrange multiplier as

λ k h = D -1 f S -A S u k,m h . ( 5
) Set u k+1,0 h = u k,m h , k = k + 1
, and go to step [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF].

Above, we denote by

u k,i h = M G(u k,i-1 h , A k n , I k n , A k τ , I k τ , u k-1,m h , λ k-1 h
) the iterate after one multigrid cycle with the actual computational refinement level of the mesh as the top level and refinement level zero as the bottom level for the linear system (3.15). We mention that we start with a given mesh on level 0. To obtain the next finer mesh on level l + 1, we decompose each element on level l into 8 subelements. For all examples presented in this paper, we use a W-cycle with three pre-and postsmoothing steps for the multigrid. As a smoother, a symmetric Gauß-Seidel iteration is applied. By A S , we denote the rows of the stiffness matrix A corresponding to the nodes in the set S. In the case m ∈ N, we solve the linear system approximately; i.e., we update the active and inactive sets after m multigrid steps. For m = "∞," we get the exact version of our algorithm. Furthermore we mention that in each Newton step we have to change only the lines belonging to the contact nodes as described by (3.15). Therefore there is no need of reassemble the whole stiffness matrix in each iteration step if we store the original part A S . Recalling the derivation of our iteration rule as a semismooth Newton method, we obtain the following local convergence result; see [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF].

Theorem 4.1. For m = "∞," Algorithm 1 converges locally superlinear.

Our numerical experience reveals that Algorithm 1 also converges globally, and thus, there is no need for a globalization strategy. This has also been observed for the application of primal-dual active set strategies for related problems in, e.g., [START_REF] Hintermüller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hüeber | Efficient algorithms for problems with friction[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF][START_REF] Kunisch | Generalized Newton methods for the 2D-Signorini contact problem with friction in function space[END_REF]. A rigorous analysis for a simper model can be found in [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF].

Example. As a first example for a 3D contact problem with Tresca friction, we consider a one-body contact problem.

The problem setting. We consider a linearly elastic cube Ω = [0, 1] 3 with material parameters E = 200 and ν = 0.3 and take the xy-plane as the rigid obstacle which implies that the initial gap is d = 0. The cube is subject to the Dirichlet distortion u = (0, 0.2, 0.06 -0.15x) on its upper surface [0, 1] 2 × {1}; see the left part of Figure 4.1. For Tresca friction, the friction bound g is given a priori and does not depend on the distortion as in the case of Coulomb friction. For this problem, we choose g = 800xy(1x)(1y). In the right of Figure 4.1, we show the distorted body with the effective von Mises stress σ eff given by σ 2 eff :=

3 i,j=1 |σ ij -δ ij p| 2
, where the pressure p is given by p := 1 3 tr(σ). The Tresca friction law. To get a better understanding of the different types of nodes that occur for contact problems with Tresca friction, we show in Figure 4.2 the nodes of the contact surface on level 5. Different types of nodes are marked differently; see the legend in Figure 4.2. Note that for each node the displacement is parallel to the multiplier λ τ as required. We remark that, using the Tresca friction law, nodes can stick in the tangential direction without being in contact with the obstacle.

Performance of Algorithm 1. To investigate the performance of the algorithm, we first solve the linear system arising in each iteration step exactly. We note that, in each step, our algorithm updates the contact/noncontact sets A k n and I k n as well as the slip/stick sets A k τ and I k τ and at the same time performs a Newton step to adopt the distortion's direction u τ to the direction of λ τ for sliding nodes. Table 4.1 Performance of Algorithm 1 for the exact strategy (m = "∞"). Initialization: u 0,0 h = 0, λ 0 h = 0 for all levels. Tolerance ε u = 10 -9 with the modification (3.11). For our tests, we initialize the algorithm on each level with u 0,0 h = 0 and λ 0 h = 0. This leads to A 1 n = A 1 τ = ∅. We terminate the IPDAS iteration if the relative change in the solution is less than 10 -9 . In the complementarity function, we use c n = c τ = 100; Algorithm 1 yields a fast and stable converge on all levels. Table 4.1 shows the number of iterations needed on different refinement levels and the number of nodes belonging to the active sets A k n and A k τ for the first modification (3.11). We remark that for the second (3.12) and the third modification (3.13) only minor differences occur. Note that in each iteration step one linear system has to be solved. The number of iterations increases only weakly on finer levels; it seems to depend linearly on the level. Usually, after the exact active sets for both friction and contact conditions are found, the method requires about 3-4 more iterations to converge. In these steps, the algorithm adjusts, for p ∈ S, the direction of the tangential traction λ τ,p to the tangential displacement u τ,p .

l |A k n | / |A k τ | for k = 2,
Comparison between the modifications. Now we compare the convergence and the behavior of the factors α k-1 p and β k-1 p,l for the three modifications. In the left of Figure 4.3, we show the errors λ k h -λ h on level 5 in a logarithmic plot for the initialization u 0,0 h = (1, 1, 0) and λ 0 h = (-1, -1, 0) . Here λ h denotes the Lagrange multiplier of the solution. We remark that additional iteration steps, compared to Table 4.1, are needed due to the different initialization and the smaller tolerance (middle) and β k-1 p,l (right) at the node (0.0625, 1, 0) ∈ S; u 0,0 h = (1, 1, 0) , λ 0 h = (-1, -1, 0) . Table 4.2 Comparison between the exact (m = "∞") and the inexact (m = 1) strategy with c n = c τ = 100 using (3.11). Tolerance for IPDAS strategy: ε u = 10 -10 . ε u = 10 -14 . In the middle of Figure 4.3, we present the cosinus α k-1 p of the angle between the vectors λ k-1 τ,p,s and λ k-1 τ,p,s + c τ u k-1 τ,p and on the right the behavior of the scaling factor β k-1 p,l , l = 1, 2, for the node p = (0.0625, 1, 0) ∈ A k τ for all k. We observe for all modifications a superlinear convergence. Comparing the behavior of the factor α k-1 p,l , the first and the second modifications show almost the same behavior. Although the factor α k-1 p,l tends faster towards 1.0 for the third modification, we observe a slower convergence. For the factor β k-1 p,l we observe a better behavior for the first modification than for the second one. From now on, we use modification one for all computations.

Strategy

Exact versus inexact method. We next compare the exact version, i.e., m = "∞," of Algorithm 1 with the inexact version. In the inexact version, we use m = 1; i.e., we update the active and inactive sets after each multigrid step. We denote by K l the iteration step in which the correct active and inactive sets are found for the first time and kept afterwards. For the inexact approach, we denote this step by M l . Table 4.2 shows the numbers K l and M l on each level and the necessary numbers of multigrid iterations to solve the full nonlinear problem on level l. We observe that the numbers K l and M l are almost the same. They seem to depend linearly on the level l. Therefore, there is no need to solve the linear system exactly. Furthermore, we compare the inexact approach, where we start with u 0,0 h = 0 and λ 0 h = 0 on each level, with the nested approach in which we inherit u 0,0 h and λ 0 h on level l + 1 from level l. The values M l and the necessary numbers of multigrid (MG) iteration steps are shown in the last column of Table 4.2. We note that the inexact primal-dual active set strategy can be interpreted as a nonlinear multigrid method.

Influence of the parameter c τ . In a last test, we investigate the influence of the parameter c τ on our algorithm. Recall that c τ can be seen as weight for the tangential distortion u τ,p in the sum with the tangential component λ τ,p,s of the Lagrange multiplier. Thus, it plays a similar role for the tangential component as c n Table 4.3 Comparison between different values for c τ for the inexact strategy (m = 1) using c n = 100 and ε u = 10 -10 . Initialization: u 0,0 h = 0, λ 0 h = 0 for all levels l. for the normal component. In Table 4.3, we compare the numbers K l and the numbers MG of necessary multigrid steps for different values of c τ , where we fix c n = 100 and use the inexact approach m = 1. As can be seen, the algorithm behaves quite stably and independently of c τ , if c τ is large enough. For this example, we find c τ ≥ 10. For c τ = 1 the algorithm does not converge on level 5. A very similar behavior is observed in [START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF] with respect to the parameter c n . Thus, in general it appears advantageous to choose both c τ and c n in order to balance the different scales of the distortion u h and the Lagrange multiplier λ h . A closer look reveals that the scaling parameters c τ and c n should reflect the material parameter.

c

Fixed point algorithm and numerical examples for Coulomb friction.

In this section, we extend Algorithm 1 to contact problems with Coulomb friction and give a numerical example. For Coulomb friction, the friction bound g p = F|λ n,p,s | needs to be iteratively adjusted using the normal component of the Lagrange multiplier. Therefore, we get an additional outer loop for the update of the friction bound. Initialize u 0,0 h and λ 0 h . (1) If mod k f (k -1) = 0, set k c = k -1 and update the friction bound by

g k c p = F max 0, λ k c n,p,s , p ∈ S.
(2) Define the active and inactive sets by

A k n := p ∈ S : λ k-1 n,p,s + c n u k-1,m n,p -d p > 0 , I k n := p ∈ S : λ k-1 n,p,s + c n u k-1,m n,p -d p ≤ 0 , A k τ := p ∈ S : λ k-1 τ,p,s + c τ u k-1,m τ,p -g k c p > 0 , I k τ := p ∈ S : λ k τ,p,s + c τ u k-1,m τ,p -g k c p ≤ 0 . ( 3
) For i = 1, . . . , m, compute u k,i h = M G(u k,i-1 h , A k n , I k n , A k τ , I k τ , u k-1,m h , λ k-1 h ).
(4) Compute the Lagrange multiplier as (5

λ k h = D -1 f S -A S u k,m h .
) If k > k f and u k,m h -u k c ,m h / u k,m h < ε u , stop. (6) Set u k+1,0 h = u k,m
h and k = k + 1, and go to step [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF].

In this algorithm, we denote by mod * * the modula-operator. Comparing this algorithm with Algorithm 1 for Tresca friction, we remark that here we update the friction bound after k f steps of the (inexact) active set strategy. This update is done in step (1). Since we do not solve the resulting linear problems exactly, it is not guaranteed that λ k n,p,s ≥ 0 for all p ∈ S. Therefore, we set g k p = F max 0, λ k n,p,s . For the choice m = k f = 1, the friction bound and the active and inactive sets are updated after each multigrid step. As the stopping criterion, we use the relative error between the actual solution u k,m h and the solution for the last friction bound u k c ,m h . For the choice m = k f = "∞," we get the exact version of the algorithm. In this case, we solve the resulting Tresca friction problem exactly for each friction bound. Obviously, this approach is rather costly. However, for a small friction coefficient F it can be shown that this discrete fixed point mapping is contractive and thus converges; see [START_REF] Nečas | On the solution of the variational inequality to the Signorini problem with small friction[END_REF].

Example.

In this section, we study the performance of Algorithm 2 for Coulomb friction. This friction model is physically more realistic than the Tresca model, since only points that are in contact with the obstacle are points for which friction occurs. Points on the contact boundary with a positive distance to the obstacle are traction-free; i.e., we apply homogeneous Neumann conditions.

Recall that, in the Coulomb law, the friction bound becomes g p = F|λ n,p,s | and thus depends on the actual distortion. In the following tests, we consider the same geometry and data as for the example of section 4.1 and choose the friction coefficient F ≡ 1. The distorted cube can be seen in the left plot of Figure 5.1. Note that its distortion is significantly different from the one obtained with the Tresca law (Figure 4.1). In the right of Figure 5.1, we visualize the constraint λ τ,p,s ≤ F|λ n,p,s | that holds for all p ∈ S. The few nodes where this inequality holds in a strict sense are nodes in contact with the obstacle that, at the same time, stick to this obstacle, i.e., u n,p = d p = 0 and u τ,p = 0. We remark that the solution has a singularity at the node (1, 1, 0) .

The Coulomb friction law. As for Tresca friction (Figure 4.2), we also visualize the different types of contact nodes for the Coulomb law. Note that, in contrast to the Tresca model, the Coulomb law allows only three types of nodes, since nodes that are not in contact with the obstacle are not subject to any (friction) constraints. The size of the nodes in Figure 5.2 is proportional to the normal contact force |λ n,p |. Performance of the algorithm. We apply the inexact version of Algorithm 2 with m = k f = 1 to solve the contact problem with Coulomb friction. We use c n = c τ = 100, ε u = 10 -10 and initialize the iteration with u 0,0 h = 0 and λ 0 h = 0 on each level. Figure 5.3 shows the behavior of our algorithm on various levels. The number of nodes contained in the active sets are plotted over the iteration steps k. We denote by k max the necessary number of iteration steps and by K l the iteration step in which the correct active sets are found for the first time and do not change afterwards. The number K l is marked by a dashed vertical line in Figure 5.3. We observe that both k max and K l appear to be almost independent of the level l. Looking at Figure 5.3 more closely, we note that on each level there are only minor changes of the active nodes after k = 10. Table 5.1 shows a comparison between the initialization u 0,0 h = 0 and λ 0 h = 0 on each level and the nested approach. As expected, the correct active sets are found earlier, and as a result fewer iterations are required for the nested approach.

Numerical examples for two-body contact with Coulomb friction.

Now we consider a curved contact interface subjected to Coulomb friction. A twodimensional cross section of our geometry is shown in Figure 6.1. The lower domain Ω m , assumed to be the master side, models a spherical shell that is fixed at the outer boundary. Against this shell, we press the body modeled by the domain Ω s , which is assumed to be the slave side. At the top surface of Ω s , we apply the Dirichlet data (0, 0, -0.2) . The geometry is given by r i = 0.7, r a = 1.0, r = 0.6, h = 0.5, and d = 0.3. In Ω s , we use a Young modulus E s = 300 and a Poisson ratio ν s = 0.3, while in Ω m we have E m = 400 and ν m = 0.3.

Discussion of the results and comparison for various friction coefficients. The results for the friction coefficient F = 0.5 are shown in Figures 6.1-6.2. Figure 6.1 shows the deformed body with the effective von Mises stress σ eff on level 3 and a two-dimensional cross section on level 2. In Figure 6.2, we show the nodes being in contact, their relative tangential slip [u τ,p ] (lines), and the tangential contact pressure λ τ,p (arrows). We remark that the nodes in the middle with only an arrow are the sticky nodes; the others are slippy. In the two pictures on the right, the normal and the tangential part of the Lagrange multiplier together with the friction bound are plotted. The comparison of the result for various friction coefficients F is shown in Figures 6. 3-6.4. In the first one, the visualization of the contact nodes together with the relative tangential slip and the tangential contact pressure are shown. contact zone on the surface of Ω s . The small oscillations in the plot occur due to the fact that all nodes p ∈ S over the whole contact zone are projected onto a straight radial line and that we work with an unstructured mesh. We observe that the number of sticky nodes decreases for smaller F. Performance of the algorithm. We consider the performance of Algorithm 2 for the friction coefficient F = 0.5. For the initialization, we set u 0,0 p = 0.1n p for p ∈ S, u 0,0 p = 0 for p / ∈ S, and λ 0 p = 0.0001n p on each level. Using the parameters c n = c τ = 100 , m = k f = 1, and ε u = 10 -9 , we get the performance of Algorithm 2 shown in Figure 6.5. Here the number of nodes in |A k n | and |A k τ | is shown for levels 1-3. Again, the dashed vertical line marks the step K l in which the correct active sets are found for the first time and remain unchanged. A comparison between this approach and the nested approach for various friction coefficients is shown in Table 6.1. Although for F = 0.2 the convergence rates of the multigrid method for levels 1 and 2 are worse, the results obtained show qualitatively the same behavior as for the cube; see section 5.1.

Problem with nonsymmetric boundary traction. For the next example, we use the same data and geometry as above (see Figure 6.1) but use a nonsymmetric boundary traction at the top of Ω s instead of the Dirichlet data; namely, we apply the surface traction (15, 0, -150 exp(-100r 2 )) , where r denotes the distance to the midpoint of the top surface of Ω s . The results for F = 0.7 are presented in Figure 6.6.

Full Newton approach for multibody contact with Coulomb friction.

While our algorithm in section 5 for the solution of the Coulomb friction problem is based on fixed point ideas, we now present a full Newton approach. The main advantage of this approach is its fast convergence, which is due to the fact that the friction bound is updated in the Newton iteration and not via a fixed point loop. In this section, we apply a fast direct solver [START_REF] Schenk | Efficient sparse LU factorization with left-right looking strategy on shared memory multiprocessors[END_REF] to solve the linear system in each Newton step. To derive the full Newton iteration, one replaces g p in (3.1) either by g p (u n,p , λ n,p,s ) := F max(0, λ n,p,s ) (see also [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][START_REF] Christensen | Reformulation: Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods[END_REF]) or by g p (u n,p , λ n,p,s ) := F max(0, λ n,p,s +c n (u n,p -d p )) before deriving the Newton iteration step. The equivalence of these two choices follows from (4.1). In what follows, we use the latter replacement, since then the Newton-type iteration automatically takes the form of an active set method that in each iteration estimates the three relevant sets for Coulomb friction (no contact, contact and stick, contact and slip). The resulting nonlinear complementarity function is

D(u h , λ h ) := max F(λ n,p,s + c n (u n,p -d p )), λ τ,p,s + c τ u τ,p λ τ,p,s -F max 0, λ n,p,s + c n (u n,p -d p ) (λ τ,p,s + c τ u τ,p ) . (7.1)
Similarly as in section 3, we derive a semismooth Newton step for the solution of D(u h , λ h ) = 0 and (4.1). Using the notation g k-1 p := F(λ k-1 n,p,s + c n (u k-1 n,pd p )), we obtain the following settings.

• On I k n := {p ∈ S : g k-1 p ≤ 0} (estimation for a set of nodes not in contact):

λ k n,p,s = 0 and λ k τ,p,s = 0.

Note that

I k n ⊂ A k τ , where A k τ := {p ∈ S : λ k-1 τ,p,s + c τ u k-1 τ,p -g k-1
p ≥ 0}. We remark that the setting for λ k τ,p,s is derived directly from (7.1) for that case, namely, from λ τ,p,s + c τ u τ,p λ τ,p,s = 0.

• On I k τ := {p ∈ S : λ k-1 τ,p,s + c τ u k-1 τ,p -g k-1
p < 0}: (estimation for a set of sticky nodes):

(7.2) u k n,p = d p and u k τ,p + Fu k-1 τ,p /g k-1 p λ k n,p,s = u k-1 τ,p . Note that I k τ ⊂ A k n , where A k n := {p ∈ S : g k-1 p > 0}. • On A k τn := A k τ ∩ A k n = {p ∈ S : λ k-1 τ,p,s + c τ u k-1 τ,p ≥ g k-1
p > 0} (estimation for a set of slippy nodes):

u k n,p = d p and -λ k τ,p,s + L k-1 p,l u k τ,p + Fv k-1 p,l λ k n,p,s = r k-1 p,l + g k-1 p v k-1 p,l , (7.3)
where Comparing this system with (3.15), we observe that here the normal and tangential components are pointwise coupled by lines four and five. Numerical example. As example, we consider the situation presented in the left of Figure 7.1, where a two-dimensional cross section of the problem definition is shown. The ring is fixed on its upper outer edge and the tool on its bottom. Note that the bodies penetrate in their reference configuration. We use Young modulus E m = 8.13 × 10 8 and a Poisson ratio ν m = 0.3 for the inner tool modeled by Ω m and E s = 9 × 10 7 and ν s = 0.3 for the outer ring being the slave domain Ω s . The friction coefficient is F = 0.7. The deformed mesh with effective von Mises stress σ eff is shown in the middle and the right of Comparison between the fixed point Newton and the full Newton approaches. To show the behavior of the full Newton approach, we compare the convergence rates of the Lagrange multiplier with those obtained from the fixed point Newton approach given by Algorithm 2 with m = "∞" and k f = 1. Our algorithm is initialized with u 0,0 h = 0 and λ 0 h = 0, and the parameters c n = c τ = 10 8 and the tolerance ε u = 10 -9 are used. The finite element mesh consists of 66.600 degrees of freedom. In the right of Figure 7.2, the relative error λ k h -λ * h / λ * h is shown for the fixed point and the full Newton approach in a logarithmic scale. We observe superlinear convergence of the full Newton approach, and only half of the iteration steps compared with the fixed point approach are required. For the fixed point approach the convergence rate

v k-1 p := (Id 2 -M k-1 p ) -1 (λ k-1 τ,p,s + c τ u k-1 τ,p )/ λ k-1 τ,p,s + c τ u k-1 τ,p ∈ R
λ k+1 h -λ * h / λ k h -λ * h
tends to approximately 0.3. In a last test we investigate the influence of the friction coefficients on the full Newton approach. The convergence rates for the example presented in Figure 6.1 for the friction coefficients F = 0.2, 0.4, 0.6 are shown in Figure 7.3. As initial values we use the same as in Figure 6.5. Again we observe the same behavior as before for all friction coefficients. Therefore the behavior of the full Newton method is independent of the size of the friction coefficient.

  where d ≥ 0 denotes the gap between the two elastic bodies, σ n := n σ s n the normal component of the boundary stress, [u] := (u s (x) -u m (R(x))) the jump and [u] n := [u]n the jump in the normal direction. For the boundary stress in the normal direction on the possible contact part, we have to satisfy the condition (2.3)

Fp, 3 = 0 and γ 2 Fp, 3 = (δ k- 1 p) 2 .

 32312 Therefore we get for the eigenvalues of the matrix M k-1 p,3 due to 0 < e k-1 p < 1 the relation 0 ≤ γ Mp,3 < 1, and the matrix Lk-1 p,3

  set p ∈ A k τ . In both cases, we set M k-1 p,l = 0, and (3.5b) leads to the desired homogeneous Neumann condition λ k τ,p,s = 0. We mention that, in the case λ k-1 τ,p,s + c τ u k-1 τ,p > 0 and λ k-1 τ,p,s > 0, we get due to (3.4b) and (3.6) automatically p ∈ A k τ and e k-1 p = 0, and therefore M k-1 p,l

(4. 1 )

 1 λ n,p,smax(0, λ n,p,s + c n (u n,pd p )) = 0, with c n > 0. The strategy uses the active and inactive sets

Algorithm 1 :

 1 IPDAS for contact with 3D Tresca friction (0) Set k = 1, choose c n > 0, c τ > 0, and m ∈ N.

Fig. 4 . 1 .

 41 Fig. 4.1. Problem definition with prescribed Dirichlet distortion (left) and distorted body with effective von Mises stress σ eff ; the lower surface in the plot is subject to unilateral contact and Tresca friction.

Fig. 4 . 2 .

 42 Fig. 4.2. Visualization of the solution at the nodes p ∈ S for Tresca friction on level 5 (left), legend (upper right), and cutout (lower right).

Fig. 4 . 3 .

 43 Fig. 4.3. Convergence of λ k h on level 5, c n = c τ = 100, and ε u = 10 -14 (left); behavior of α k-1 p,l

Algorithm 2 :

 2 Fixed point IPDAS for contact with 3D Coulomb friction (0) Set k = 1 and choose c n > 0, c τ > 0, m ∈ N, and k f ∈ N.

Fig. 5 . 1 .

 51 Fig. 5.1. Distorted body with effective von Mises stress σ eff (left); the lower surface in the plot is subject to contact and Coulomb friction. Visualization of the friction bound F|λ n,p | (small balls) and of λ τ,p (right).

Fig. 5 . 2 .

 52 Fig. 5.2. Visualization of the solution at the nodes p ∈ S for Coulomb friction on level 5 (left), legend (upper right), and cutout (lower right).

Fig. 5 . 3 .

 53 Fig. 5.3. Behavior of Algorithm 2: Numbers |A k n | and |A k τ | of active nodes in each iteration step k on each level l. We used the parameters c n = c τ = 100, m = k f = 1, and ε u = 10 -10 . Initialization: u 0,0 h = 0 and λ 0 h = 0.

Fig. 6 . 1 .

 61 Fig. 6.1. Problem definition (left), deformed mesh with effective von Mises stress σ eff on level 3 (middle), and two-dimensional cross section of the deformed mesh on level 2 (right) for F = 0.5.

Fig. 6 . 2 .Friction coefficient F=0. 6 Fig. 6 . 3 .Fig. 6 . 4 .

 6266364 Fig. 6.2. Visualization of the nodes being in contact together with the deformation and stress vectors on level 3 for F = 0.5 (first) and a cutout (second)-for the legend we refer to Figure 5.2and visualization of λ n (third) and λ τ (fourth) with the friction bound 0.5|λ n | (dotted line).

Fig. 6 . 5 .

 65 Fig. 6.5. Behavior of Algorithm 2: Numbers |A k n | and |A k τ | in each iteration step k on levels 1, 2, 3 for F = 0.5. We use the parameters c n = c τ = 100, m = k f = 1, and ε u = 10 -9 .

Fig. 6 . 6 .

 66 Fig. 6.6. Nonsymmetric boundary data: Deformed mesh with effective von Mises stress σ eff on level 3 (left); visualization of the nodes being in contact together with the deformation and the stress vectors (middle) and of λ τ (right) with the friction bound 0.7|λ n | (dotted line).

2 and g k- 1 p 1 p 1 pFig. 7 . 1 .

 11171 Fig. 7.1. Problem definition (left), deformed mesh with effective von Mises stress σ eff (middle), and two-dimensional cross section of the deformed mesh (right) for F = 0.7. solved in each full Newton step for the Coulomb problem as ⎛ ⎜ ⎜ ⎜ ⎜ ⎝

Figure 7 . 1 .

 71 The first two pictures of Figure7.2 show the possible contact nodes on the ring being the slave side. The nodes without a line are nodes not being in contact with the inner tool.

Fig. 7 . 2 .

 72 Fig. 7.2. Visualization of the nodes being in contact, together with the deformation and stress vectors for F = 0.7 (left) and a cutout (middle)-for the legend we refer to Figure 5.2-and convergence rates for the fixed point Newton and the full Newton approaches (right).

Fig. 7 . 3 .

 73 Fig. 7.3. Convergence rates for the example of Figure 6.1 for the fixed point Newton and the full Newton approaches for F = 0.2, 0.4, 0.6.

  we obtain G C (u τ,p , λ τ,p,s )(δu τ,p , δλ τ,p,s ) = max(g p , λ τ,p,s + c τ u τ,p )δλ τ,p,s + c τ δu τ,p )g p (δλ τ,p,s + c τ δu τ,p ) . A denotes the characteristic function of the set A τ = {p ∈ S : λ τ,p,s + c τ u τ,p > g p }, i.e., χ A = 1 if λ τ,p,s + c τ u τ,p > g p and χ A = 0 if λ τ,p,s + c τ u τ,p ≤ g p . We note that λ τ,p,s (λ τ,p,s + c τ u τ,p ) is a 2 × 2-matrix, either zero or of rank one.

	+ χ A (δλ τ,p,s (3.2) λ τ,p,s (λ τ,p,s + c τ u τ,p ) λ τ,p,s + c τ u τ,p
	Here χ Performing a semismooth Newton step at a current iterate (u k-1 τ,p , λ k-1 τ,p,s ), one derives
	the new iterates (u k τ,p , λ k τ,p,s ) from
	(3.3)	Solve G C (u k-1 τ,p , λ k-1 τ,p,s )(δu k τ,p , δλ k τ,p,s ) = -C(u k-1 τ,p , λ k-1 τ,p,s ) and update (u k τ,p , λ k τ,p,s ) = (u k-1 τ,p , λ k-1 τ,p,s ) + (δu k τ,p , δλ k τ,p,s ).
	The characteristic function χ A in (3.2) separates the nodes of S into the inactive set I k τ and the active A k τ set according to

  MG K l MG K l MG M l MG M l MG

	τ	1		10			100		1000	10000
	Level l K l 1 3	12	3	11	3	11	3	11	3	11
	2	3	12	3	12	3	13	3	13	4	13
	3	4	14	4	14	4	14	4	14	4	14
	4	6	15	6	14	6	14	6	14	6	14
	5	-	-	8	17	8	17	9	17	9	17

Table 5 . 1

 51 Behavior of Algorithm 2: K l and k max (see text for explanation) on each level l for c n = c τ = 100, m = k f = 1, and ε u = 10 -10 .

	Level l DOF	u 0,0 h = 0, λ 0 h = 0 Nested approach K l k max K l k max
	0	8	10	32	10	32
	1	27	9	38	3	34
	2	125	10	45	5	32
	3	729	21	47	3	37
	4	4913	12	46	4	38
	5	35937	14	45	5	36

  1 n p , λ 0 p = 0.0001 n p

										Nested approach
				K l			k max			K l			k max
			0.2 0.4 0.6 0.2 0.4 0.6	0.2 0.4 0.6 0.2 0.4 0.6
	0	104	4	4	4	12	16	21	4	4	4	12	16	21
	1	541	10	5	5	111 20	33	1	3	2	91	20	25
	2	3384	5	7	6	67	27	26	2	3	2	62	21	22
	3	23694	9	11	9	33	32	32	3	3	2	29	21	21

* This work wa s supported in part by the Deutsche Forschungsgemeinschaft, SFB 404, B8, and SPP 1146. http://www.siam.org/journals/sisc/30-2/67106.html