Jonas Koko

an automatic penalty parameter selection

Keywords: Linear elasticity, unilateral contact problem, alternating direction, augmented Lagrangian

Alternating direction method of multiplier for the unilateral contact problem with

Introduction

Unilateral contact (frictionless) problem is a very common problem in engineering and poses serious challenges. It differs from the classical linear elasticity only by the presence of a linear constraint (the non-penetration condition). Various numerical methods have been developed, see, e.g., [START_REF] Kikuchi | Contact problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Hüeber | A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction[END_REF][START_REF] Ito | Semi-smooth Newton methods for the Signorini problem[END_REF][START_REF] Koko | Uzawa block relaxation for the unilateral contact problem[END_REF][START_REF] Stadler | Path-following and augmented Lagrangian methods for contact problems in linear elasticity[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF] and references therein. The method proposed in [START_REF] Koko | Uzawa block relaxation for the unilateral contact problem[END_REF] is an alternating direction method of multiplier (ADMM) and it is related to augmented Lagrangian operator splitting methods [START_REF] Fortin | Augmented Lagrangian Methods: Application to the Numerical Solution of Boundary-Value Problems[END_REF][START_REF] Glowinski | Augmented Lagrangian and Operator-splitting Methods in Nonlinear Mechanics[END_REF]. The main idea is to separate the differentiable part of the problem (i.e. linear elasticity) from the nondifferentiable part (i.e. the non-penetration condition) by introducing an auxiliary unknown. Applying an alternating direction method to the corresponding augmented Lagrangian leads to a simple two-step iterative method in which the matrix of the linear elasticity is constant. This property makes ADMM algorithm outperforms, see [START_REF] Koko | Uzawa block relaxation for the unilateral contact problem[END_REF], another optimization based method, the popular active-set semi-smooth Newton method [START_REF] Hüeber | A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction[END_REF][START_REF] Ito | Semi-smooth Newton methods for the Signorini problem[END_REF][START_REF] Stadler | Path-following and augmented Lagrangian methods for contact problems in linear elasticity[END_REF].

Alternating direction method of multiplier is a powerful operator-splitting algorithm for solving structured optimization problems encoutoured in a wide range of areas such as compressed sensing [START_REF] Yang | Alternating direction algorithms for l 1 -problems in compressive sensing[END_REF], image restoration [START_REF] Koko | An augmented lagrangian method for T V g +L 1norm minimization[END_REF], machine learning [START_REF] Boyd | Distributed optimization and statistical learning via alternating direction method of multipliers[END_REF], etc. The method was initiated by Glowinski et al. [START_REF] Glowinski | Sur l'approximation par éléments finis d'ordre un, et la résolution par pénalisation-dualité, d'une classe de problèmes de Dirichlet non linéaires[END_REF][START_REF] Fortin | Augmented Lagrangian Methods: Application to the Numerical Solution of Boundary-Value Problems[END_REF][START_REF] Glowinski | Augmented Lagrangian and Operator-splitting Methods in Nonlinear Mechanics[END_REF], who systematically developed ADMM algorithms (also known as Uzawa block relaxation methods) for solving nonlinear partial differential equations. But ADMM based algorithms have a severe drawback: the choice of the penalty parameter. Indeed, the rate of convergence of the algorithm is strongly related to the penalty parameter and the choice of this parameter is a non trivial problem. We propose, in this paper, a simple procedure for the numerical approximation of the optimal penalty parameter for use with ADMM algorithm in the case of a frictionless contact problem. To this end, we carefully derived a convergence analysis on a pure dual ADMM algorithm. The pure dual version is obtained by eliminating the primal and auxiliary unknowns. Note that the approximate formula is based on extreme eigenvalues, and is therefore suitable for coarse meshes only.

The paper is organized as follows. In Section 2 the ADMM algorithm for the unilateral frictionless contact problem is presented. In Section 3 the finite dimensional versions are presented. The convergence analysis is carried out in Section 4. The penalty parameter approximation is outlined in Section 5, followed by the numerical experiments in Section 6.

The model problem and ADMM algorithm

We consider an elastic body occupying in its initial (undeformed) configuration a bounded domain Ω of R d (d = 2, 3) with a boundary Γ = Γ D ∪Γ C . We assume that the elastic body is fixed along Γ D with meas(Γ D ) > 0. Γ C denotes a portion of Γ which is a candidate contact surface between Ω and a rigid foundation. The normalized gap between Γ C and the rigid foundation is denoted by g. In this paper, we consider the small strains hypothesis so that the strain tensor is E(u) = (∇u+∇u t )/2, where u = (u 1 (x), . . . , u d (x)) is the displacement field. Hooke's law is assumed, i.e. the stress tensor is linked to the displacement through the linear relation

σ(u) = CE(u)
where C = (C ijkl ) is the (fourth order) elastic moduli tensor, assumed to be symmetric positive definite. Let n be the outward unit normal to Ω on Γ. We consider the normal component of the displacement field and the stress tensor given by

u n = u • n, σ n (u) = (σ(u)n) • n.
The unilateral contact problem consists, for a given volume force f , of finding the displacement field u satisfying (i) the equilibrium equations

-divσ(u) = in Ω, (2.1)
u = 0 on Γ D , (2.2)
(ii) and the contact (i.e. non-penetration) conditions

u n -g ≤ 0, σ n (u) ≤ 0, (u n -g)σ n (u) = 0, on Γ C . (2.3)
Let us introduce the functions space

V = v ∈ H 1 (Ω) 2 , v = 0 on Γ D and the set of admissible displacements K = {v ∈ V, v n -g ≤ 0 on Γ C }. Let a(•,
•) be the symmetric, continuous and coercive bilinear form which corresponds to the virtual work in the elastic body

a(u, v) = Ω σ(u) : E(v)dx.
We denote by f (•) the linear form of external forces

f (v) = Ω f • vdx.
Using the potential energy functional

J(v) = 1 2 a(v, v) -f (v)
we can formulate the unilateral contact problem as the constrained minimization problem

Find u ∈ K such that J(u) ≤ J(v), ∀v ∈ K. (2.4)
To achieve a solution of (2.4) by an ADMM algorithm, we need additional steps. Following Glowinski and Le Tallec [START_REF] Glowinski | Augmented Lagrangian and Operator-splitting Methods in Nonlinear Mechanics[END_REF], we introduce the set

C = p ∈ L 2 (Γ C ), p -g ≤ 0 on Γ c
and its characteristic functional

I C : L 2 (Γ c ) → R ∪ {+∞} defined by I C (p) = 0 if p ∈ C, +∞ if p ∈ C.
It is clear that (2.4) is equivalent to the following constrained minimization problem

Find (u, p) ∈ V × L 2 (Γ C ) 2 such that J(u) + I C (p) ≤ J(v) + I C (q) ∀(v, q) ∈ V × L 2 (Γ C ) 2 ,
(2.5)

u n -p = 0 on Γ C . (2.6)
With (2.5)-(2.6) we associate the augmented Lagrangian functional

L r (v, q; λ) = J(v) + I C (q) + (v n -q, λ) Γ C + r 2 v n -q 2 Γ C , (2.7)
where r > 0 is the penalty parameter. ADMM algorithm for (2.4) is based on (2.7). The method performs successive minimizations in u and p followed by the multiplier λ update as follows

u k+1 = arg min v L r (v, p k ; λ k ), (2.8
)

p k+1 = arg min p L r (u k+1 ; q; λ k ), (2.9)
λ k+1 = λ k + r(u k+1 n -p k+1 ).
(2.10) Minimization (2.8) leads to an equilibrium problem which always has a unique solution even though the original problem allows rigid body motions. Minimization (2.9) is solved explicitly. The corresponding ADMM algoririthm is presented in Algorithm 1 (see [START_REF] Koko | Uzawa block relaxation for the unilateral contact problem[END_REF] for detailed derivation). We iterate until the relative error on u k and p k becomes "sufficiently" small.

Algorithm 1 ADMM algorithm for the unilateral contact problem Initialization. φ 0 and λ 0 are given.

Iteration k ≥ 0. Compute successively u k+1 , p k+1 and λ k+1 as follows

• Find u k+1 ∈ V such that a(u k+1 , v) + r(u k+1 n , v n ) Γ C = f (v) + (rpk -λ k , v n ) Γ C , ∀v ∈ V.
• Compute the auxiliary unknown

p k+1 = u k+1 n + 1 r λ k -λ k + r(u k+1 n -g) + .
• Update the Lagrange multiplier

λ k+1 = λ k + r(u k+1 n -p k+1 ).

Finite dimensional algorithms

We assume that Ω ⊂ R d is a polyhedral domain and therefore can be exactly triangulated. We consider a (continuous piecewise linear) finite element triangulation T h of Ω consistent with the decomposition of its boundary Γ into Γ D and Γ C . If n is the number of nodes of the triangulation, the dimension of the finite element subsapce V h ⊂ V is dimV h = dn. Let x j (j = 1, . . . , m) be a contact node (i.e. x j lies on Γ C ). The displacement vector and the unit outward normal at x j are denoted by u j and n j , respectively. We introduce the linear mappings N : R dn → R m , such that Nu is the vector of the normal components of u at contact nodes, that is, (Nu) j = u j n j , j = 1, . . . , m.

The finite element discretization leads to the following matrices and vectors

• A, (dn) × (dn) stiffness matrix (symmetric positive definite), i.e. from the bilinear form a(u, v);

• M normal mass matrices (m × m symmetric positive definite);

• f ∈ R dn the (discrete) external forces;

• g ∈ R m the (discrete) normalized gap;

• p ∈ R m the (discrete) auxiliary unknown;

• λ ∈ R m the (discrete) lagrange multiplier.
With the notations above, the augmented Lagrangian (2.7) now reads in the discrete setting

L r (u, p; λ) = J(u) + λ M(Nu -p) + r 2 (Nu -p) M(Nu -p), (3.1)
where

J(u) = 1 2 u Au -f u.

Discrete ADMM algorithm

The derivation of the discrete version of Algorithm 1 from the augmented Lagrangian function (3.1) is straigthforward, Algorithm 2.

Algorithm 2 Algebraic ADMM algorithm for the unilateral contact problem Initialization k = 0 r > 0, p 0 and λ 0 are given.

Iteration k ≥ 0 . Compute successively u k+1 , p k+1 and λ k+1 as follows.

Step 1 Compute u k+1 ∈ V h such that

(A + rN MN)u k+1 = f + N M(rp k -λ k ) (3.2)
Step 2 Compute p k+1 ∈ P h

p k+1 = Nu k+1 + 1 r λ k -(λ k + r(Nu k+1 -g)) + (3.3)
Step 3 Update Lagrange multiplier

λ k+1 = λ k + r(Nu k+1 -p k+1 )

ADMM algorithm without the auxiliary unknown

Algorithm 2 can be simplified by eliminating the auxiliary unknown p. Indeed, using the update formula (3.3) we have

rp k = r(p k -Nu k ) + rNu k = -(λ k-1 + r(Nu k -p k ) + rNu k + λ k-1 = rNu k + λ k-1 -λ k , that is rp k -λ k = λ k-1 -2λ k + rNu k . (3.4) If we set λ k -= min(0, λ k-1 + r(Nu k -g))
, the right-hand side of (3.4) becomes

λ k-1 + rNu k -2λ k = λ k-1 + r(Nu k -g) + rg -2λ k = λ k --λ k + rg. It follows that rp k -λ k = λ k --λ k + rg i.
e., we have eliminated the auxiliary unknown p from the right-hand side of (3.2). It remains to remove p from the multiplier update formula. From (3.3) we deduce that

r(p k+1 -Nu k+1 ) = λ k -(λ k + r(Nu k+1 -g)) + .
(3.5) Substituting (3.5) into the Lagrange multiplier update formula, we obtain the new multiplier update formula λ k+1

+ := (λ k + r(Nu k+1 -g)) + .
Gathering the results above, the ADMM algorithm without the auxiliary unknown p is described in Algorithm 3. In Algorithm 3, λ 0 + and λ 0 -must be consistent. Indeed, if, e.g., u 0 = 0, then we can set λ 0 + = r(Nu 0 -g) + = 0. It follows that λ 0 -= r(Nu 0 -g) -= -rg.

Algorithm 3 ADMM without auxiliary unknown for the unilateral contact problem Initialization k = 0 r > 0, λ 0 + and λ 0 -are given.

Iteration k ≥ 0 . Compute successively u k+1 , λ k+1 + and λ k+1 -as follows.

Step 1 Compute u k+1 such that

(A + rN MN)u k+1 = f + rN Mg + N M(λ k --λ k + ) (3.6)
Step 2 Update the Lagrange multipliers

λ k+1 + = (λ k + + r(Nu k+1 -g)) + λ k+1 - = (λ k + + r(Nu k+1 -g)) -

Pure dual version

The pure dual version of the ADMM algorithm for the unilateral contact problem is obtained by eliminating the displacements vector u. From (3.6) we have

u k+1 = A -1 r f + rN Mg + N M(λ k --λ k + )
where we have set

A r = A + rN MN.
Substituting u k+1 into the Lagrange multiplier update formulas, we obtain Algorithm 4.

λ k+1 + = (I -rNA -1 r N M)λ k + + rNA -1 r N Mλ k -+ r(Na -g) + (3.7) λ k+1 - = (I -rNA -1 r N M)λ k + + rNA -1 r N Mλ k -+ r(Na -g) - (3.8)
Algorithm 4 is unpraticable since in A -1 r can be a full matrix of large size. However, the recursive relation (3.7)-(3.8) is important because it gives an easy way to prove the convergence of the ADMM algorithm. Indeed, Algorithm 4 is useful as a theoretical basis for the study of the influence of r.

Convergence

Since min(0, x) = -max(0, -x), the recursive relation (3.7)-(3.8) can be rewritten in the useful form

λ k+1 + = (I -rNA -1 r N T M)λ k + + rNA -1 r N T Mλ k -+ r(Na -g) + (4.1) λ k+1 - = --(I -rNA -1 r N T M)λ k + -rNA -1 r N T Mλ k --r(Na -g) + . (4.2)
Since for any x ∈ R n x + ≤ x we have

λ k+1 + 2 ≤ (I -rNA -1 r N M)λ k + + rNA -1 r N Mλ k - 2 + r(Na -g) 2 (4.3) λ k+1 - ≤ -(I -rNA -1 r N M)λ k + -rNA -1 r N T Mλ k - 2 + r(Na -g) 2 (4.4)
so that

λ k+1 + 2 + λ k+1 - 2 ≤ (I -rNA -1 r N M)λ k + + rNA -1 r N T Mλ k - 2 + -(I -rNA -1 r N M)λ k + -rNA -1 r N Mλ k - 2 +2 r(Na -g) 2 that is λ k+1 ± 2 ≤ Aλ k+1 ± 2 +2 r(Na -g) 2 ,
where we have set

λ k+1 ± = λ k+1 + λ k+1 - and A = I -rNA -1 r N T M rNA -1 r N T M -(I -rNA -1 r N T M) -rNA -1 r N T M .
The convergence of Algorithm 4 therefore depends on the spectral radius of A. We have the following proposition for the non-zeros eigenvalues of A.

Proposition 4.1 The non-zero eigenvalues of A are eigenvalues of I -2rNA -1 r N M and conversely.

Proof. Let us set

A 1 = I -rNA -1 r N M and A 2 = rNA -1 r N M so that A = A 1 A 2 -A 1 -A 2 .
We have to show that the non-zero eigenvalues of A are eigenvalues of A 1 -A 2 = I -2rNA -1 r N M and conversely. Let α be an eigenvalue of A and (X 1 , X 2 ) the corresponding eigenvector. We have

A 1 X 1 + A 2 X 2 = αX 1 , (4.5)
-A 1 X 1 -A 2 X 2 = αX 2 (4.6)
that is, α(X 1 + X 2 ) = 0. For non-zero eigenvalues we deduce that the eigenvectors are such that X 1 = -X 2 . Substituting in (4.5)-(4.6) we obtain

(A 1 -A 2 )X 1 = αX 1 , (A 1 -A 2 )X 2 = αX 2 .
We deduce that α is an eigenvalue of A 1 -A 2 . From the relations above, it is obvious that non-zero eigenvalues of A 1 -A 2 are also eigenvalues of A. 2

To study the eigenvalues of I -2rNA -1 r N M, we consider the following sequence

µ k+1 = (I -2rNA -1 r N M)µ k = µ k -2rNA -1 r N Mµ k , k ≥ 0. (4.7)
We follow the same procedure as [3, p. 46] for a quadratic programming problem with equality constraints. From (4.7), we deduce that

A -1 N Mµ k+1 = A -1 N Mµ k -2rA -1 N MNA -1 r N Mµ k . (4.8) Since A -1 r = (A + rN MN) -1 = (I + rA -1 N MN) -1 A -1 it follows that A -1 N Mµ k+1 = A -1 N Mµ k -2rA -1 N MN(I+rA -1 N MN) -1 A -1 N Mµ k (4.9)
By introducing the auxiliary unknown η k = A -1 N Mµ k , the sequence (4.9) becomes

η k+1 = (I -2rA -1 N MN(I + rA -1 N MN) -1 )η k , ∀k ≥ 0. (4.10)
The convergence of the sequence {η k } will gives us informations about the convergence of the sequence {µ k } and then {λ k + , λ k -}. From [5, th. 1.1.6 and cor. 2.4.

3.4], B = (I -2rA -1 N MN(I + rA -1 N MN) -1 ) is a polynomial function of A -1 N MN, with p(t) = (1 -2rt(1 + rt) -1
) as polynomial. Then the eigenvalues of B are polynomial functions of the eigenvalues of A -1 N MN and thus λ i (B) = p(λ i (A -1 N MN), that is

λ i (B) = 1 -2r λ i (A -1 N MN) 1 + rλ i (A -1 N MN) = 1 -rλ i (A -1 N MN) 1 + rλ i (A -1 N MN) . (4.11)
The eigenvalues (and eigenvectors) of A -1 N MN are solution to the generalized symmetric eigenvalue problem N MNw = λAw Since A is symmetric positive definite and N MN is symmetric positive semi-definite, the eigenvalues of A -1 N MN are non-negative, and the eigenvectors corresponding to two distinct eigenvalues are A-otthogonal. It follows that the eigenvalues of B are such that 0 ≤ λ i (B) ≤ 1, ∀i. (4.12)

Note that A and M are square and non singular matrices. Then if 0 is an eigenvalue of A -1 N MN , then the corresponding eigen-subspace is ker(N). Thus R(A -1 N M) is spanned by the eigenvectors of A -1 N MN associated with strictly positive eigenvalues (using the property R(N) = (ker(N)) ⊥ ). Since {η k } ⊂ R(A -1 N M), we can write (4.10) in a basis of R(A -1 N M) formed with eigenvectors w i of A -1 N MN associated with strictly positive eigenvalues λ i . We get

η k+1 i = 1 -rλ i 1 + rλ i η k i , ∀i = 1, . . . , N 1 and k ≥ 0.
We deduce the following convergence theorem. Corollary 4.3 (optimal step-size and convergence rate) The optimal penalty parameter for r is

r * = 1 √ λ m λ M (4.13)
where λ m and λ M are, respectively, the smallest nonzero eigenvalue and the largest eigenvalue of A -1 N MN. With the choice (4.13), the convergence of Algorithm 4 is linear with an asymptotic constant θ satisfying

θ ≤ 1 -(λ m /λ M ) 1/2 1 + (λ m /λ M ) 1/2 . ( 4
κ 2 (A -1 N MN) = A -1 N MN 2 (N MN) -1 A 2 = λ M /λ m .
We obtain the following corollary.

Corollary 4.4 If N MN is nonsingular, then the eigenvalues of A -1 N MN are strictly positive and Algorithm 4 with the optimal choice (4.13) converges linearly with an asymptotic constant θ satisfying

θ ≤ κ 2 (A -1 N MN) -1 κ 2 (A -1 N MN) + 1

Penalty parameter approximation

We now study how to compute λ m and λ M , the smallest nonzero eigenvalue and the largest eigenvalue of A -1 N MN. For this purpose, it is convenient to arrange the indices in such a way that the contact indices C and the interior indices I occur in a consecutive order. This leads to the block matrix representation of the stiffness matrix as

A = A II A IC A CI A CC and N MN = O II O IC O CI M CC ,
where O II , O IC and O CI are matrices of zeros; and M CC is the contact boundary (diagonal) mass matrix. Consequently, if we set

A -1 = ÃII ÃIC ÃCI ÃCC ,
then the nonzero eigenvalues of A -1 N MN are eigenvalues of ÃCC M CC . To compute ÃCC , we can apply block-Gaussian elimination to A with A II as block-pivot

ÃCC = A CC -A CI A -1 II A IC -1 .
It follows that ÃCC is the Schur complement of A II in A.

We now detail the practical steps for the approximation of the penalty parameter. Since computing A -1 II is unpracticable even for medium size problem, we compute ÃCC by solving equivalent linear systems. We obtain Algorithm 5.

Algorithm 5 Algorithm for computing the approximate penalty parameter r *

Step 1. Solve for X the system A II X = A IC (5.1)

Step 2. Compute the schur complement of

A II S = A CC -A CI X Step 3. Solve for X SX = M CC (5.2)
Step 4. Compute λ m and λ M the exterme eigenvalues of X and the optimal penalty parameter approximation (4.13).

Solving (5.1) or (5.2) is equivalent to solving several linear systems that have the same matrix but different right-hand sides. Since the different right-hand sides are available, all the systems of equations can be solved at the same time using a single Gaussian elimination (e.g. using MATLAB \ operator : X = S\M CC for (5.2)).

Numerical experiments

We have implemented the algorithms described in the previous sections in MATLAB, using piece-wise linear finite element, vectorized assembling functions and the mesh generator provided in [START_REF] Koko | A MATLAB mesh generator for the two-dimensional finite element method[END_REF][START_REF] Koko | Fast MATLAB assembly of fem matrices in 2d and 3d using cell array approach[END_REF], on a computer equipped running Linux (Ubuntu 16.04) with 3.00GHz clock frequency and 32GB RAM. The ADMM solver used in the experiments is based on Algorithm 1. The test problems used are designed in order to illustrate the behavior of the algorithms more than to model contact actual phenomena. Since our approximation procedure combines matrix inversion and computation of extreme eigenvalues, we restrict our penalty parameter estimation to small size problems.

We use the following notations

• r * approximate optimal penalty parameter obtained by (4.13).

• r standard numerical optimal penalty parameter obtained using the uniform sampling of interval (.1, αE), where α > 1 and E is the Young modulus.

• m c the size of matrix ÃCC M CC used for computing r * . To reduce the computational cost m c ≤ 100 in the numerical experiments.

• Iter. The number of iterations required for convergence in Algorithm 2.

• CPU Time. CPUT time in Seconds.

Example 1: 2D problem

We consider a rectangular elastic body Ω = (0, 2) × (0, 1) with the boundary partition We first consider the meshes of size h=1/8 (Figure 1), 1/16 and 1/32. We compare the approximate penalty parameter r * obtained with (4.13) with the numerical penalty parameter r obtained by sampling the interval (0.1, 10) * E with 200 points uniformly spaced. The evolution of the number of iterations with respect to the penalty parameter is shown in Figure 2. The normal stress distribution on Γ C is shown in Figure 3 We summarize in Table 1 the (optimal) penalty parameters and the behavior of Al-gorithm 2 in terms of CPU time and number of iterations. Note that the CPU time for r * includes the time spent in both Algorithm 2 and Algorithm 5. The CPU time for r is the computational time of the whole sampling procedure. We can notice that the sampling process gives better results in terms of the number of iterations required by the ADMM solver. But sampling is a time-consuming process. For the largest problem (h = 1/32), obtaining the optimal penalty parameter (4.13) with Algorithm 5 is about 80 times faster than the standard sampling procedure. Table 2 shows that, for r * = 3134.1289 and r = 3510 (obtained using a mesh of size 1/32), the performances of Algorithm 2 are almost comparable for large size meshes. .

Γ D = {0} × (0,

Example 2: Hertz problem

The Hertz problem is a classical test problem in the numerical simulation of unilateral contact problem. It consists of an infinitely long cylinder resting in a rigid foundation, and subjected to a uniform load along its top of intensity f = -(0, 1600). The radius of the cylinder is R = 8. The material constants for the cylinder are E = 2000 and ν = 0.3. For symmetry reason, only quarter of the cylinder is consider (Figure 4) and we set u 1 = 0 on Γ D = {0} × (0, 8). The contact surface is Γ C = {x ∈ (0, 8) 2 | x 2 1 + x 2 2 = 16}. We first consider meshes with 442 and 1692 with 35 and 69 nodes on Γ C , respectively. Note that the problem allows a rigid body motion in the vertical direction. In ADMM Algorithm 2, the mass terms provided by the normal integral prevent infinite displacements at the initial step. But for Algorithm 5, A II is only positive semidefinite and then non invertible. We simply add a small positive constant on the vertical components of Γ D . As in Example 1, we compute r by sampling (.1, 10) * E with 200 points, uniformly spaced. Figure 5 shows the variation of the number of iterations with respect to the penalty parameter. The deformed configuration and the stress distribution on the contact Γ C are shown in Figure 6 and 7, respectively. We summarize in Table 3 the (optimal) penalty parameters and the behavior of Algorithm 2 in terms of CPU time and number of iterations. For m c = 35, the sampling procedure gives better results in terms of the number of iterations for the ADMM solver but the computational time is much greater. For m c = 69 the number of iterations is almost equivalent for both r and r * but the computational time is again much greater for r. We report in Table 4 the performances of Algorithm 2 using r * = 201.7825 and r = 200. We can notice that the performances of Algorithm 2 using both penalty parameters are almost equivalent.

Example 3: 3D problem

We now study the behavior of a 3D rectangular elastic body pressed onto a solid hemisphere, Figure 8 ( [START_REF] Stadler | Path-following and augmented Lagrangian methods for contact problems in linear elasticity[END_REF][START_REF] Koko | Uzawa block relaxation for the unilateral contact problem[END_REF]). The elastic body occupies the domain Ω = (-0.5, 0.5)×(-1, 1)× (-0.5, 0.5) and the obstacle is a half-ball with radius r = 0.5 and center (-0.3, 0, -1). The contact surface is the lower surface Γ C = (-0.5, 0.5) × (-1, 1) × {-0.5}. A displacement u D = (0, 0, -0.2) is prescribed on the upper surface Γ D = (-0.5, 0.5) × (-1, 1) × {0.5}. The Young modulus is E = 10 6 and the Poisson ratio ν = 0.45. Figure 8 shows geometry and the obstacle. We first consider a non uniform mesh with 198 nodes, 773 tetrahedrons and 45 nodes on Γ C . The deformed configuration obtained using ADMM Algorithm 2 is shown in Figure 9 and the normal (contact) stress distribution in Figure 10. For r, we sample (0.1 10 * E) with 200 points, uniformly spaced. We summarize in Table 5 the computed penalty parameters r * and r and the behavior of Algorithm 2 in terms of number of iterations and CPU time. Even though r * is 20% larger than r, the number of iterations required by the ADMM solver with both values is almost the same. We notice again that computing the penalty parameter with Algorithm 5 is about 100 times faster. reported in Table 6. As in the 2D case, the number of iterations is virtually independent of the mesh size, for both penalty parameters. We can notice that, for both penalty parameter, the difference of computational time for the largest problem is less than 5%. The performances of Algorithm 2 with r * and r are therefore comparable.

Conclusion

We have investigated the automatic parameter selection for the alternating direction method of multiplier applied to the unilateral contact problem. Numerical experiments have shown that the proposed procedure for the parameter selection is a good alternative to the standard sampling procedure.

Algorithm 4

 4 Dual ADMM algorithm for the unilateral contact problem Initialization k = 0 r > 0, λ 0 + and λ 0 -are given. Set A r = A + rN MN and a = A -1 r (f + rN Mg) Iteration k ≥ 0 . Compute successively λ k+1 + and λ k+1 -as follows.

Theorem 4 . 2 (

 42 Convergence) ADMM Algorithm 4 converges for any r > 0.

  1), Γ N = (0, 2) × {1} and Γ C = (0, 2) × {0}. On Γ D the body is clampsed (u D = 0) and on Γ N non-homogeneous Neumann boundary conditions are prescribed σ(u) • n = -10. On Γ C , the normalized gap betwenn Ω and the rigid foundation is g(x) = 0.01. The material constants are E = 2000 and ν = 0.3.

Figure 1 :

 1 Figure 1: Mesh sample for Ω = (0, 2) × (0, 1) of size h = 1/8

Figure 2 :

 2 Figure 2: Number of iterations versus the penalty parameter for Example 1.

Figure 3 :

 3 Figure 3: Normal Stress distribution on Γ C for Example 1.

Figure 4 :

 4 Figure 4: Mesh sample for the Hertz problem.

Figure 5 :

 5 Figure 5: Number of iterations versus penalty parameter for the hertz problem.

Figure 6 :

 6 Figure 6: Deformed configuration and Von Mises effective stress for the Hertz problem

Figure 7 :

 7 Figure 7: Normal Stress distribution on Γ C for the Hertz probel

Figure 8 :

 8 Figure 8: Geometry of the 3D obstacle problem

Figure 9 :

 9 Figure 9: Deformed configuration and Von Mises effective stress for Example 3.

Table 1 :

 1 r obtained by sampling Versus r * obtained by (4.13).

	h		1/8	1/16	1/32
	size of m c	17	33	65
	r *		1591.4878 2220.5867 3134.1289
	Iter.	35	27	27
	CPU Time	0.057	0.177	0.269
	r		3510	2810	3510
	Iter.	19	18	19
	CPU Time	1.797	4.957	23.766
			.
	h		1/8	1/16 1/32 1/64 1/128 1/256
	r * = 3134.1289	Iterations CPU Time 0.009 0.018 0.096 0.479 1.603 9.127 24 19 27 25 19 21
	r = 3510	Iterations CPU Time 0.008 0.014 0.078 0.339 1.646 8.713 19 20 19 20 20 20

Table 2 :

 2 Performances of Algorithm 2 with r * = 3134.1289 and r = 3510 for Example 1.

Table 3 :

 3 r obtained by sampling Versus r * obtained by (4.13), for the Hertz problem.

	.

Table 4 :

 4 Performances of Algorithm 2 with r * = 201.7825 and r = 200 for the Hertz problem .

Table 5 :

 5 r obtained by sampling Versus r * obtained by (4.13) for the 3D obstacle problem To study the beahivior of Algorithm 2, the initial mesh of 198 nodes and 773 tetrahedrons is successively refined to produce meshes with 1289, 9245 and 69897 nodes and 6144, 49472 and 395776 tetrahedrons, respectively. The performances of Algorithm 2 areFigure 10: Normal stress distribution on Γ C Number of nodes on Ω/Γ C 198/45 1289/153 9245/603 69897/2341

	r * = 3.1026 × 10 6	Iter. CPU Time	28 0.053	30 0.613	31 10.267	32 363.534
	r = 2.6000 × 10 6	Iter. CPU Time	26 0.046	28 0.577	33 10.568	30 349.122

Table 6 :

 6 Performances of Algorithm 2 with r and r * for the 3D obstacle problem.