Keywords: Absorbing boundary conditions, waveguide, finite element, periodic medium

This paper proposes a new method, in the frequency domain, to define absorbing boundary conditions for general two-dimensional problems. The main feature of the method is that it can obtain boundary conditions from the discretized equations without much knowledge of the analytical behavior of the solutions and is thus very general. It is based on the computation of waves in periodic structures and needs the dynamic stiffness matrix of only one period in the medium which can be obtained by standard finite element software. Boundary conditions at various orders of accuracy can be obtained in a simple way. This is then applied to study some examples for which analytical or numerical results are available. Good agreements between the present results and analytical solutions allow to check the efficiency and the accuracy of the proposed method.

Finite element computation of absorbing boundary conditions for time-harmonic wave problems

Denis Duhamel a, * , Tien-Minh Nguyen a

1 Introduction

Wave problems in unbounded media can occur in many applications in mechanics and engineering such as in acoustics, solid mechanics, electromagnetics, etc. It is well known that analytical solutions for such problems are available only for some special cases. On the contrary, numerical methods can be applied to many complex problems. Physically, for problems in infinite domains, the energy is produced by sources in the region to be analyzed and must escape to infinity. For methods solving the problem on a bounded domain like the finite element method, it introduces the difficulty of an artificial boundary to get a bounded domain. This boundary must be such that the energy crosses it without reflection and special conditions must be specified at the artificial boundary to reproduce this phenomena. Generally, these can be classified into local or global boundary conditions. With a global condition all degrees of freedom (dofs) on the boundary are coupled while a local condition connects only neighboring dofs.

The first global method which has been used for solving such problems was the boundary element method. This method is well adapted for infinite domains and is described in numerous classical textbooks like [START_REF] Brebbia | Boundary elements techniques in engineering[END_REF][START_REF] Crouch | Boundary element methods in solid mechanics[END_REF][START_REF] Ciskowski | Boundary element methods in acoustics[END_REF][START_REF] Chen | Boundary element methods, Computational mathematics and applications[END_REF][START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF]. It consists in solving an equation on the boundary of the domain only and the radiation conditions are taken into account analytically. It also reduces the dimension of the problem to a surface in 3D and to a curve in 2D decreasing thus the size of the linear problem to solve. However, the final problem involves full matrices which are also generally non symmetrical. It is also mainly limited to linear problems and to homogeneous domains or otherwise one has to introduce special and complex techniques to deal with non linear or non homogeneous situations. There are also singularities in the integrals which need special attention for the numerical integrations. So this method is interesting and has been extensively used but it can lead to heavy computations when the number of degrees of freedom increases. More information on such techniques can be found in the historical and review papers [START_REF] Cheng | Heritage and early history of the boundary element method[END_REF][START_REF] Beskos | Boundary element methods in dynamic analysis: Part II[END_REF].

In the other approaches, the computational domain is truncated at some distance and boundary conditions are imposed at this artificial boundary. These conditions at finite distance must simulate as closely as possible the exact radiation condition at infinity. An approach leading to a global boundary condition is the Dirichlet to Neumann (DtN) mapping proposed by [START_REF] Keller | A finite element method for large domains[END_REF][START_REF] Keller | Exact non reflecting boundary conditions[END_REF] and in an earlier version by [START_REF] Hunt | Finite element approach to acoustic radiation from elastic structures[END_REF][START_REF] Hunt | Finite element approach to acoustic scattering from elastic structures[END_REF]. It consists in dividing the domain into a finite part containing the sources and an infinite domain of simple shape. The solution in the infinite domain is solved analytically, for example by series expansions, and an exact impedance relation is obtained on the boundary between the finite and infinite domains. This relation links the variable and its normal derivative on the whole boundary. The DtN mapping is thus non local and every node on the boundary is connected to all other nodes. This gives a full matrix for the nodes of the boundary which partially destroys the sparse matrix of the FEM and increases substantially the computing resources needed to get the solution. The solution has to be found in the exterior domain by analytical or numerical methods. When the analytical solution can be found, it is generally under the form of a series expansion. The number of terms in the expansion must be sufficient for an accurate solution which can lead to heavy computations. Developments of the method can be found in [START_REF] Harari | Finite element methods for the Helmholtz equation in an exterior domain: model problems[END_REF][START_REF] Thompson | Finite element formulation of exact non-reflecting boundary conditions for the time-dependent wave equation[END_REF]. An application to the case of wave scattering in plates is also found in [START_REF] Galan | Elastodynamic guided wave scattering in infinite plates[END_REF].

The other methods are local and the condition at a node involves only neighboring nodes which make them less demanding in computing resources and much easier to implement in a finite element code but also less accurate. A first possibility of such approaches is the use of infinite elements as proposed by [START_REF] Bettess | Infinite elements[END_REF][START_REF] Bettess | More on infinite elements[END_REF][START_REF] Bettess | Numerical Methods in Coupled Systems[END_REF][START_REF] Bettess | Infinite elements[END_REF][START_REF] Astley | Infinite elements for wave problems: a review of current formulations and an assessment of accuracy[END_REF]. It consists in developing special elements with a behavior at infinity reflecting that of analytical solutions obtained for the same type of problems. For wave problems, it involves complex-valued basis functions with outwarding propagation wave-like behavior in the radial direction. The elements were further developed by [START_REF] Burnett | A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion[END_REF] to considered other coordinate systems such as expansions in prolate coordinates. This method is interesting but the inclusion of infinite elements requires the development of special elements and these elements can depend on decay parameters which have to be accurately chosen. A review of these methods has been proposed by [START_REF] Gerdes | A review of infinite element methods for exterior Helmholtz problems[END_REF].

In the perfectly matched layer proposed by [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF], originally for electromagnetic waves, an exterior layer of finite thickness is introduced around the bounded domain. The absorption in this domain is increasing as we move towards the exterior such that outgoing waves are absorbed before reaching the exterior truncation boundary. The number of elements in the layer, its thickness, the variation of the absorption properties have to be carefully chosen to optimize the efficiency of the method. This efficiency is better for a layer with a large thickness but this can lead to a significant increase in the number of elements in the finite element model. Various developments of the method can be found in [START_REF] Colino | Optimizing the perfectly matched layer[END_REF] and its optimization in [START_REF] Asvadurov | On optimal finitedifference approximation of PLM[END_REF]. Otherwise, various classes of absorbing boundary conditions were also developed by [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF]. They consist in the numerical approximation of differential operators on the boundary. For instance, examples of the application of Bayliss-Turkel conditions are presented by [START_REF] Strouboulis | The generalized finite element method for Helmholtz equation. Part II: Effect of choice of handbook functions, error due to absorbing boundary conditions and its assessment[END_REF]. However, the more accurate boundary conditions involve high order derivatives on the boundary which are difficult to implement in the finite element method [START_REF] Pinsky | Local high-order radiation boundary conditions for the two-dimensional time-dependent structural acoustics problem[END_REF]. Finally, it was proved by [START_REF] Asvadurov | On optimal finitedifference approximation of PLM[END_REF][START_REF] Guddati | Continued fraction absorbing boundary conditions for convex polygonal domains[END_REF] that, in fact, the perfectly matched layer and the absorbing boundary conditions were closely connected. The Helmholtz equation was also solved with these two boundary conditions by [START_REF] Heikkola | Fast direct solution of the Helmholtz equation with a perfectly matched layer or an absorbing boundary condition[END_REF] and these conditions were compared and optimized to minimize the reflection. Other boundary conditions involving only second order derivatives have also been proposed. They introduce auxiliary variables and systems of equations on the boundary which lead to high order boundary conditions, see [START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF] for a review of such non-reflecting boundary (NRBC) methods. They were mainly developed for acoustic problems but in [START_REF] Krenk | Local tensor radiation conditions for elastic waves[END_REF] a local boundary condition for elastic waves has been proposed. In [START_REF] Guasch | An algebraic subgrid scale finite element method for the convected Helmholtz equation in two dimensions with applications in aeroacoustics[END_REF] an impedance boundary condition in new coordinates was developed for the convected Helmholtz equation.

For fluid dynamic problems, [START_REF] Storti | Dynamic boundary conditions in computational fluid dynamics[END_REF] developed Lagrange multipliers for imposing various absorbing boundary conditions for cases where the type and the number of boundary conditions can change, for instance as the flow changes from subsonic to supersonic regimes and its direction varies with time. A general review of the methods described in the precedent paragraphs for various dynamic, acoustic and wave propagation problems can be found in [START_REF] Givoli | Non-reflecting boundary conditions: a review[END_REF][START_REF] Mesquita | Numerical methods for the dynamics of unbounded domains[END_REF][START_REF] Harari | A survey of finite element methods for time-harmonic acoustics[END_REF][START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF].

Comparisons are also made between the different methods.

In the present study, another local method is proposed. This method works on discrete systems directly, in contrast with many existing absorbing boundary conditions which are written on the continuous differential equations and discretized after. The principle of the method is to compute wave propagations in groups of elements near the boundary from the dynamic stiffness matrix of these elements. Then, a boundary condition is obtained for cancelling the reflected waves. This condition is finally written as an impedance boundary condition relating the force and displacement degrees of freedom on the boundary. The approach is based on the waveguide theory proposed by [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF][START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF][START_REF] Duhamel | Finite element computation of Green's functions[END_REF] and is used to determine absorbing boundary conditions at the truncation boundary of 2D periodic media. Only information related to one period, obtained from any standard FE software (the discrete stiffness and mass matrices and the nodal coordinates) are required to formulate the method. The advantage is that it can be applied to media with various complex behaviors.

This paper is outlined as follows. In section 2, the methodology for determining absorbing boundary conditions for periodic media is described. Then, a discussion for the application of the method to general media is proposed. In section 3, a simple application is proposed to show the results of the method in a case where detailed computations can be done. In section 4, two examples of finite element computations in acoustics and elastodynamics are presented. They allow to check the efficiency and accuracy of the proposed method for more complex cases. Finally, the paper is closed with some conclusions.

Absorbing boundary conditions

We suppose that we want to solve a mechanical problem on an infinite domain exterior to the bounded domain Ω int (see figure 1). The infinite domain is approximated by the finite domain Ω which is exterior to Ω int and is limited by the exterior boundary Γ ext . We are looking for a solution with radiating condition at infinity which means that the solution should be outgoing near the boundary Γ ext . Near this exterior boundary the solution can be seen as composed of incident waves denoted A + and reflected waves A -. For a perfectly absorbing boundary, one should have A -= 0. In fact this condition is very difficult to implement in the numerical solutions of such problems. Indeed, only the global solution is easily computed but the decomposition into incident and reflected waves is difficult to obtain. The problem is thus to find an appropriate boundary condition to impose on the exterior boundary to finally get A -≈ 0 on the solution. To be easily included in a finite element model the searched boundary condition should be local and the condition at a node of the boundary should involve only neighboring nodes.

The approach proposed in this paper consists in studying this problem by first considering the case of periodic media. For this case, positive and negative waves and their amplitudes A + and A -can be computed by the method presented below. Then an exact boundary condition can be formulated for a half-plane boundary. It is further shown how this condition can be approximated by a local condition on the boundary. As homogeneous media are special cases of periodic media, the method presented here applies also to homogeneous media. Before considering the general case, a simple example for the Klein Gordon equation will be presented.

A simple example

Consider first the stationary Klein Gordon equation given by

d 2 u dx 2 + (k 2 -m 2 )u = 0 ( 1
)
where u is the solution and k, m are real parameters. This equation is discretized with linear two nodes elements such that

u(ξ) = u 1 N 1 (ξ) + u 2 N 2 (ξ) (2)
where

N 1 (ξ) = (1 -ξ)/2, N 2 (ξ) = (1 + ξ)/2,
k = - 1 l    1 -1 -1 1    m = l 6    2 1 1 2    (3)
and the dynamic stiffness matrix of one element is given by

d = - 1 l    1 -1 -1 1    + l 6 (k 2 -m 2 )    2 1 1 2    (4)
Waves of propagating constant e µ are such that

u 2 = e µ u 1 (5) f 2 + e µ f 1 = 0 (6)
leading in an element to

e
with

d 11 d 12 = - 1 -l 2 (k 2 -m 2 )/3 1 + l 2 (k 2 -m 2 )/6 (10)
For l 2 (k 2m 2 ) << 1, one gets, at first order,

e µ ± ≈ 1 ± il √ k 2 -m 2 (11)
meaning

µ ± ≈ ±il √ k 2 -m 2 (12)
There are thus two waves in an element, such that

   u + f +    =    1 d 11 + e µ + d 12    and    u - f -    =    1 d 11 + e µ -d 12    (13)
and the general solution is given by

   u f    = a +    u + f +    + a -    u - f -    (14)
The condition for only outgoing waves is thus a -= 0 on the right boundary and a + = 0 on the left boundary leading respectively to the conditions

f /u = f + /u + on the right f /u = f -/u -on the lef t (15)
The condition in the first case is

f /u = d 11 + e µ + d 12 = d 11 +   - d 11 d 12 + ( d 11 d 12 ) 2 -1   d 12 = d 2 11 -d 2 12 = -(k 2 -m 2 ) + 1 12 ((k 2 -m 2 )l) 2 = i √ k 2 -m 2 1 - 1 12 (k 2 -m 2 )l 2 (16)
≈ i √ k 2 -m 2 (17)
In the second case, one gets

f /u = d 11 + e µ -d 12 ≈ -i √ k 2 -m 2 (18)
We recognize approximations of the classical absorbing boundary conditions which have been obtained here directly from the discretized equations. Compared to the classical boundary condition on the right (and exact in this case)

f /u = i √ k 2 -m 2 , the relative error is 1 -1 12 (k 2 -m 2 )l 2
which depends mainly on the size of the element relatively to the wavelength. The present boundary condition has been obtained entirely from the discrete matrices without any knowledge of the analytical solution of the problem.

To estimate the reflection coefficient created by such a boundary condition consider an incident wave Ae

i √ k 2 -m 2 x on the boundary. A reflected wave RAe -i √ k 2 -m 2 x
is created. The total solution and its associated force are given by

u(x) = Ae i √ k 2 -m 2 x + RAe -i √ k 2 -m 2 x f (x) = Ai √ k 2 -m 2 e i √ k 2 -m 2 x -RAi √ k 2 -m 2 e -i √ k 2 -m 2 x (19)
Writing the boundary condition ( 16), for instance by taking the boundary at x = 0, yields

f (0) u(0) = 1 -R 1 + R = 1 - 1 12 (k 2 -m 2 )l 2 (20)
So the reflection coefficient is finally given by

R = 1 -1 -1 12 (k 2 -m 2 )l 2 1 + 1 -1 12 (k 2 -m 2 )l 2 ≈ 1 48 (k 2 -m 2 )l 2 (21)
This coefficient is low and of second order when √ k 2m 2 l << 1.

General impedance boundary condition

In this section we present the general outline of the method before starting a more rigorous developement in the following section. So, to extend the precedent example to more general cases, consider a vector function u(x, y) and a force vector f(x, y) acting on a line parallel to the y axis as in figure 2. They can be decomposed by a Fourier transform as

u(x, y) = 1 (2π) 2 +∞ -∞ +∞ -∞ u(k x , k y )e i(kxx+kyy) dk x dk y f(x, y) = 1 (2π) 2 +∞ -∞ +∞ -∞ f(k x , k y )e i(kxx+kyy) dk x dk y (22)
Let us suppose that u is solution of a linear operator

L(u) = n=N n=0 α 1 +α 2 =n a α 1 α 2 ∂ n u ∂x α 1 ∂y α 2 = 0 (23)
In the Fourier domain, this relation yields

n=N n=0 α 1 +α 2 =n a α 1 α 2 (ik x ) α 1 (ik y ) α 2 u(k x , k y ) = 0 (24)
For a given value of k y , the precedent relation has non zero solutions for k x such that the determinant

n=N n=0 α 1 +α 2 =n a α 1 α 2 (ik x ) α 1 (ik y ) α 2 = 0 (25)
Let us denote by k + j the n + positive solutions such that Re(k + j ) < 0 or Re(k + j ) = 0 and the energy flux is directed towards positive values of x. We denote by k - j the other solutions. We have the decomposition

u(x, k y ) = j=n + j=1 a + j u + j e ik + j x + j=n - j=1 a - j u - j e ik - j x (26)
In the same way, for the force components

f(x, k y ) = j=n + j=1 a + j f + j e ik + j x + j=n - j=1 a - j f - j e ik - j x (27)
where f + j and f - j are the force components respectively associated to u + j and u - j . If the boundary is such that only positive waves exists at proximity, one has

u(0, k y ) = j=n + j=1 a + j u + j = U + a + f(0, k y ) = j=n + j=1 a + j f + j = F + a + (28)
where U + and F + are the matrices whose columns are respectively u + j and f + j . Eliminating the a + coefficients, one gets

f(0, k y ) = F + (U + ) -1 u(0, k y ) f(0, y) = (G * u)(0, y) (29)
with

G(k y ) = F + (U + ) -1 (k y ) ( 30
)
and * means the convolution.

In the following this boundary condition will be computed directly from the discrete equations for general linear media.

Solution in a periodic medium

Consider an infinite two dimensional periodic medium, as shown in figure 3.

The elementary period is limited by the domain (

x 1 , x 2 ) ∈ [0, b 1 ] × [0, b 2 ]. A function U(x 1 , x 2
) defined on the two-dimensional periodic medium can be decomposed as an integral of pseudo periodic functions

U(x 1 , x 2 ) = π b 2 -π b 2 e ikx 2 Û (x 1 , k, x 2 )dk (31)
where

Û(x 1 , k, x 2 ) is a periodic function in x 2 with period b 2 . From the Fourier transform Û(x 1 , k) of U(x 1 , x 2 ), one has Û (x 1 , k, x 2 ) = 1 2π +∞ m 2 =-∞ Û (x 1 , k + 2π m 2 b 2 )e i2π m 2 b 2 x 2 = 1 2π +∞ m 2 =-∞ e i2π m 2 b 2 x 2 +∞ -∞ e -i(k+2π m 2 b 2 )x U(x 1 , x)dx = 1 2π +∞ -∞ e -ikx U(x 1 , x) +∞ m 2 =-∞ e i2π m 2 b 2 (x 2 -x) dx = b 2 2π +∞ -∞ +∞ m 2 =-∞ δ(x 2 -x -m 2 b 2 )e -ikx U(x 1 , x)dx = b 2 2π +∞ m 2 =-∞ e -ik(x 2 +m 2 b 2 ) U(x 1 , x 2 + m 2 b 2 ) ( 32
)
This gives the relation inverse of [START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF]. From relation [START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF], one sees that the behavior in x 2 of the general solution can be obtained from functions as e ikx 2 Û (x 1 , k, x 2 ) with Û(x 1 , k, x 2 ) periodic in x 2 . Along direction 1, we use a decomposition in Bloch waves as it is usual in periodic media. Finally, the general solution can be obtained from functions

u(x 1 , k, x 2 ) = e ikx 2 Û (x 1 , k, x 2 ) such that: u(x 1 , k, x 2 + m 2 b 2 ) = e ikm 2 b 2 u(x 1 , k, x 2 ) (33) u(x 1 + m 1 b 1 , k, x 2 ) = e im 1 µ u(x 1 , k, x 2 ) ( 34
)
where m 1 and m 2 are integers,

k ∈ R ∩ - π b 2 , π b 2 and µ ∈ C.
The discrete dynamic equation of a cell (an elementary period) obtained from a FE model at a frequency ω and for the time dependence e -iωt is given by:

(K -iωC -ω 2 M)q = f (35)
where K, M and C are the stiffness, mass and damping matrices, respectively, f is the loading vector and q the vector of the degrees of freedom (dofs).

Introducing the dynamic stiffness matrix D = K -iωCω 2 M, decomposing the dofs into boundary (B) and interior (I) dofs, and assuming that there are no external forces on the interior nodes, result in the following equation:

   D BB D BI D IB D II       q B q I    =    f B 0    (36)
The interior dofs can be eliminated using the second row of equation [START_REF] Mesquita | Numerical methods for the dynamics of unbounded domains[END_REF], which results in

q I = -D -1 II D IB q B ( 37
)
The first row of equation [START_REF] Mesquita | Numerical methods for the dynamics of unbounded domains[END_REF] becomes

f B = D BB -D BI D -1 II D IB q B ( 38
)
which can be written as

f = Dq (39)
It should be noted that only boundary dofs are considered in the following.

The periodic cell is assumed to be meshed with an equal number of nodes on their opposite sides. The boundary dofs are decomposed into left (L), right (R), bottom (B), top (T ) dofs and associated corners (LB), (RB), (LT ) and (RT ) as shown in figure 4. The longitudinal dofs vector is defined as q l = t [ t q L t q R t q LB t q RB t q RT t q LT ]. Thus, equation ( 39) is rewritten as

       D ll D lB D lT D Bl D BB D BT D T l D T B D T T               q l q B q T        =        f l f B f T        (40)
Using the pseudo periodic condition [START_REF] Storti | Dynamic boundary conditions in computational fluid dynamics[END_REF] and the effort equilibrium at the bottom side of the cell, relations between the transverse dofs are given by

q T = e ikb 2 q B f B + e -ikb 2 f T = 0 (41)
Multiplying the third row of equation ( 40) with e -ikb 2 , taking the sum of the second and third rows of equation( 40), using conditions [START_REF] Duhamel | Finite element computation of Green's functions[END_REF], lead to

D Bl + e -ikb 2 D T l q l + D BB + D T T + e -ikb 2 D T B + e ikb 2 D BT q B = 0(42) so q B = -D BB + D T T + e -ikb 2 D T B + e ikb 2 D BT -1 D Bl + e -ikb 2 D T l q l (43)
Using ( 41) and ( 43), the first row of equation ( 40) becomes

f l = D ll -D lB + e ikb 2 D lT D BB + D T T + e -ikb 2 D T B + e ikb 2 D BT -1 × D Bl + e -ikb 2 D T l q l (44)
which can be written as

f l = D l q l (45)
Using the pseudo periodic conditions (34) also lead to the following relations between longitudinal dofs q R = e iµ q L q RB = e iµ q LB q RT = e i(µ+kb 2 ) q LB q LT = e ikb 2 q LB (46)

From the pseudo periodic conditions (46), it can be seen that all components of the vector q l depend on the set of dofs defined by q r = t [ t q L t q LB ]. This can be expressed as

q l = W 0 + e iµ W 1 q r ( 47
)
where the matrices W 0 and W 1 depend on the wavenumber k and are given by

W 0 =                  I O O O O I O O O O O e ikb 2 I                  W 1 =                  O O I O O O O I O e ikb 2 I O O                  (48)
The equilibrium conditions between adjacent cells are given by

e iµ f L + f R = 0 e iµ f LB + f RB + e i(µ-kb 2 ) f LT + e -ikb 2 f RT = 0 (49)
that can be written as

e iµ W * 0 + W * 1 f l = 0 (50)
where (.) * denotes the operator of complex conjugate and transpose.

Combining (45), ( 47) and (50), lead to

e iµ W * 0 + W * 1 D l W 0 + e iµ W 1 q r = 0 (51)
that can be written as

A 0 + e iµ (A 1 + A 2 ) + e 2iµ A 3 q r = 0 ( 52
)
where

A 0 = W * 1 D l W 0 A 1 = W * 0 D l W 0 A 2 = W * 1 D l W 1 A 3 = W * 0 D l W 1 (53)
The eigenvalue e iµ and the eigenvector q r are thus solutions of a quadratic eigenvalue problem. It is convenient to transform the problem (52) into another linear eigenvalue problem as

e iµ    A 3 O O A 3       q r q r    =    O A 3 -A 0 -(A 1 + A 2 )       q r q r    (54)
with q r = e iµ q r .

From equations ( 44) and (45), one can notice that

t D l (k) = D l (-k) (55)
Moreover, from (48), we have

W * j (k) = t W j (k) = t W j (-k) for j = 0, 1 (56)
and from (53)

t A 0 (k) = A 3 (-k) t A 1 (k) = A 1 (-k) t A 2 (k) = A 2 (-k) t A 3 (k) = A 0 (-k) (57)
It can be easily shown by taking the determinant of the matrix in relation (52) that if e iµ j is an eigenvalue for the wavenumber k, e -iµ j is also an eigenvalue for the wavenumber -k. These represent a pair of positive and negative-going waves, respectively. The 2n eigensolutions of equation ( 54) can be split into two sets of n + and n -eigensolutions with 2n = n + +n -, which are denoted by e iµ + j , q + j and e iµj , q - j respectively, with the first set such that e iµ + j ≤ 1.

In the case e iµ + j = 1, the first set of positive-going waves must contain waves propagating in the positive direction such that Re iωq H j f r j > 0 where f r j is the reduced set of boundary force dofs of left cells on right cells and is given by

f r j =    f L f LB + e -ikb 2 f LT    = W * 0 f l j = W * 0 D l W 0 + e iµ j W 1 q j (58)
In the second set of negative-going waves, the eigenvalues e iµ - j are associated with waves such that Re iωq H j f r j < 0.

With the eigenvector q j and the force component of relation (58), we introduce the state vector

x j (k) =    q j (k) f r j (k)    =    q j (k) (A 1 (k) + e iµ j (k) A 3 (k))q j (k)    (59)
In this relation q j (k) is the eigenvector associated to e iµ j (k) . One can also introduce

y j (-k) = t p j (-k)(A 2 (k) + e iµ j (k) A 3 (k)) t p j (-k) (60)
In this relation p j (-k) is the eigenvector associated to e -iµ j (k) since we have seen that e -iµ j (k) is also an eigenvalue of (52) for the wavenumber -k. From relation (52) written for the eigenvector q j (k), multiplying this relation by e -iµ j (k) and then on the left by t p i (-k), one gets

t p i (-k) e -iµ j (k) A 0 (k) + (A 1 (k) + A 2 (k)) + e iµ j (k) A 3 (k) q j (k) = 0 (61)
In the same way, writing relation (52) for the eigenvector p i (-k), taking the transpose of the relation, using relations (57) and multiplying on the right by q j (k), leads, after a global multiplication by e iµ i (k) , to the following relation

t p i (-k) e iµ i (k) A 3 (k) + (A 1 (k) + A 2 (k)) + e -iµ i (k) A 0 (k) q j (k) = 0 (62)
The difference between the two precedent relations yields

(e iµ i (k) -e iµ j (k) ) t p i (-k) A 3 (k) -e -iµ i (k) e -iµ j (k) A 0 (k) q j (k) = 0 (63)
In the case e iµ i (k) = e iµ j (k) , we get

t p i (-k) e iµ i (k) A 3 (k) -e -iµ j (k) A 0 (k) q j (k) = 0 ( 64
)
Now it is possible to compute the product y i (-k).x j (k) by

y i (-k).x j (k) = t p i (-k)(A 2 (k) + e iµ i (k) A 3 (k))q j (k) + t p i (-k)(A 1 (k) + e iµ j (k) A 3 (k))q j (k) = t p i (-k)(A 2 (k) + e iµ i (k) A 3 (k))q j (k) -t p i (-k)(A 2 (k) + e -iµ j (k) A 0 (k))q j (k) = d i δ ij (65)
The result of relation (64) has been used in the case e iµ i (k) = e iµ j (k) and d i is a factor depending on the eigenvector i. This gives orthogonality relations on the statevectors associated to the eigenvalues.

Absorbing boundary conditions

Figure 2 presents the periodic medium near the exterior boundary. In this domain the solution is described by relation [START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF], yielding, respectively for the displacement and force components,

q(x 1 , x 2 ) = π b 2 -π b 2 q(x 1 , k, x 2 )e ikx 2 dk f(x 1 , x 2 ) = π b 2 -π b 2 f(x 1 , k, x 2 )e ikx 2 dk (66)
with the force components given by relation (58). Introducing the state vector x = t ( t q, t f) and decomposing this solution into the different waves, we get

x(x 1 , x 2 ) = π b 2 -π b 2 x(x 1 , k, x 2 )e ikx 2 dk ≈ π b 2 -π b 2 j=2n j=1 a j (x 1 , k)x j (k)e ikx 2 dk (67)
The last relation is the approximation obtained by the finite element computation of wave solutions presented before. The condition of outgoing waves means that there is no incoming wave, so the amplitudes a j (x 1 , k) associated with incoming waves must equal zero. This condition is obtained by

y - l (-k). j=2n j=1 a j (x 1 , k)x j (k) = 0 f or 1 ≤ l ≤ n - (68)
In this relation y - l (-k) are the vectors associated to the negative going waves, given by relation (60). Using relation (65), one gets a - j (x 1 , k) = 0 for 1 ≤ j ≤ n -for the amplitudes of the negative going waves. Introducing the matrix Y with lines given by y - l leads to

Y(-k).x(x 1 , k, x 2 ) = 0 (69)
Decomposing now x into its displacement and force components, doing the same thing for

Y(-k) = [Q(-k) F(-k)] leads to Q(-k).q(x 1 , k, x 2 ) + F(-k). f (x 1 , k, x 2 ) = 0 (70)
The relation on the boundary is

f (x 1 , k, x 2 ) = -F -1 (-k)Q(-k)q(x 1 , k, x 2 ) ( 71
)
and then from relation (66)

f(x 1 , x 2 ) = - π b 2 -π b 2 F -1 (-k)Q(-k)q(x 1 , k, x 2 )e ikx 2 dk (72)
From the inverse relation [START_REF] Krenk | Local tensor radiation conditions for elastic waves[END_REF], one also has

q(x 1 , k, x 2 ) = b 2 2π +∞ m 2 =-∞ e -ik(x 2 +m 2 b 2 ) q(x 1 , x 2 + m 2 b 2 ) ( 73
)
which leads to

f(x 1 , x 2 ) = - b 2 2π π b 2 -π b 2 F -1 (-k)Q(-k) +∞ m 2 =-∞ e -ik(x 2 +m 2 b 2 ) q(x 1 , x 2 + m 2 b 2 )dk(74)
Introducing the function

G(x 2 ) = - b 2 2π π b 2 -π b 2 F -1 (-k)Q(-k)e -ikx 2 dk (75)
The final relation is

f(x 1 , x 2 ) = +∞ m 2 =-∞ G(x 2 + m 2 b 2 )q(x 1 , x 2 + m 2 b 2 ) (76)
This is the impedance relation on the boundary obtained with the assumption that there is no negative going wave. This relation is the absorbing boundary condition we were looking for. It can be computed from the wave vectors and the force components associated with them. Relation (76) involves an infinite number of terms on the boundary. This relation can also be written as

f(x 1 , x 2 ) = +∞ m 2 =-∞ G(x 2 + m 2 b 2 ) q(x 1 , x 2 ) + +∞ m 2 =-∞ G(x 2 + m 2 b 2 ) (q(x 1 , x 2 + m 2 b 2 ) -q(x 1 , x 2 )) = +∞ m 2 =-∞ G(x 2 + m 2 b 2 ) q(x 1 , x 2 ) + 1 2b 2 +∞ m 2 =-∞ m 2 b 2 G(x 2 + m 2 b 2 ) (q(x 1 , x 2 + b 2 ) -q(x 1 , x 2 -b 2 )) + +∞ m 2 =-∞ G(x 2 + m 2 b 2 )[q(x 1 , x 2 + m 2 b 2 ) -q(x 1 , x 2 ) - 1 2 m 2 (q(x 1 , x 2 + b 2 ) -q(x 1 , x 2 -b 2 ))] (77)
If q is slowly varying the last term should be small and for practical purposes we will use the approximate relations at various orders given by

f(x 1 , x 2 ) ≈ G 0 q(x 1 , x 2 ) + G 1 2b 2 (q(x 1 , x 2 + b 2 ) -q(x 1 , x 2 -b 2 )) + G 2 2b 2 2 (q(x 1 , x 2 + b 2 ) + q(x 1 , x 2 -b 2 ) -2q(x 1 , x 2 )) + . . . ( 78
)
with

G 0 = +∞ m 2 =-∞ G(x 2 + m 2 b 2 ) = -(F -1 Q)(0) G 1 = +∞ m 2 =-∞ m 2 b 2 G(x 2 + m 2 b 2 ) = i(F -1 Q) ′ (0) G 2 = +∞ m 2 =-∞ (m 2 b 2 ) 2 G(x 2 + m 2 b 2 ) = (F -1 Q) ′′ (0) (79)
Relation (78) involves a finite number of nodes around the point where the relation is written. It depends on the number of nodes chosen to approximate the boundary condition. This number can be 1 for a crude approximation involving only one node or can be larger. For a very large number of nodes, the condition tends towards the true absorbing condition for a half-plane in the periodic media given by (76). Up to now everything has been written for periodic media but it is clear that homogeneous media are also periodic media and so all that has been said applies also to homogeneous media. The condition (78) can be seen as a generalization of the Taylor approximation boundary condition proposed by [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF]. But, while the boundary conditions in [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] were obtained by approximation of the exact continuous relations for specific problems, they are obtained here directly and with general applicability from the discretized equations.

3 Simple examples

Estimation of the accuracy

In this section we try to estimate the quality of the proposed boundary condition compared with known relations for the simple case of the two-dimensional acoustics. Consider first a plane wave incident on the plane y = 0 at an angle θ with the normal to the plane. Let us define points at a horizontal distance D from the origin and with a vertical spacing h, see figure 5. The sound pressure at point (D, lh) is given by p l a = e iK(cos θD+sin θlh) (80

)
where K = ω/c is the wavenumber and c is the sound velocity. The analytical force at the same points is the normal derivative in direction 1 given by f l a = iK cos θe iK(cos θD+sin θlh) (81)

For a point source at origin, the pressure is solution of

∆p + K 2 p = -δ(r) ( 82
)
where r is the distance from the origin and δ(.) is the Dirac function. The solution of this equation for the time dependence e -iωt is given by

G(r) = i 4 H 0 (Kr) (83) 20
where H 0 is the Hankel function of zero order and first type. The analytical solution at each point l is

p l a = i 4 H 0 (K D 2 + (lh) 2 ) (84)
and the analytical force at the same points is the normal derivative (in direction 1) given by

f l a = - iKD 4 D 2 + (lh) 2 H 1 (K D 2 + (lh) 2 ) ( 85
)
The absorbing boundary condition described in the precedent section will allow to compute numerical forces f l n at a node from the knowledge of p l a . If the boundary condition was perfect one would have f l n = f l a but the proposed condition is approximate and one only has f l n ≈ f l a . The error can be estimated by

e = |f l n -f l a | |f l a | (86)
The next step is to compute f l n from the method proposed in this paper and the error by relation (86) to estimate the quality of the absorbing condition.

Acoustic element

Consider the rectangular four nodes acoustic element of size b 1 × b 2 . The elementary stiffness and mass matrices are given by

K = 1 6b 1 b 2           (2b 2 2 + 2b 2 1 ) (-2b 2 2 + b 2 1 ) (-b 2 2 -b 2 1 ) (b 2 2 -2b 2 1 ) (-2b 2 2 + b 2 1 ) (2b 2 2 + 2b 2 1 ) (b 2 2 -2b 2 1 ) (-b 2 2 -b 2 1 ) (-b 2 2 -b 2 1 ) (b 2 2 -2b 2 1 ) (2b 2 2 + 2b 2 1 ) (-2b 2 2 + b 2 1 ) (b 2 2 -2b 2 1 ) (-b 2 2 -b 2 1 ) (-2b 2 2 + b 2 1 ) (2b 2 2 + 2b 2 1 )           (87) M = b 1 b 2 36c 2           4 2 1 2 2 4 2 1 1 2 4 2 2 1 2 4           (88)
and the dynamic stiffness matrix can then be determined by D = Kω 2 M.

It can be noted that the reduced set of displacement dofs q r contains only q LB . Then, the matrices W 0 and W 1 have the following forms:

W 0 = t 1 0 0 e ikb 2 ; W 1 = t 0 1 e ikb 2 0 ( 89
)
The terms A j in equation ( 53) are given by

A 0 (k) = W * 1 (k)D l (k)W 0 (k) = - 1 18 12 b 2 b 1 -6 b 1 b 2 + 2K 2 b 1 b 2 + 6 b 2 b 1 + 6 b 1 b 2 + K 2 b 1 b 2 cos(kb 2 ) A 1 (k) = W * 0 (k)D l (k)W 0 (k) = - 1 9 -6 b 2 b 1 -6 b 1 b 2 + 2K 2 b 1 b 2 + -3 b 2 b 1 + 6 b 1 b 2 + K 2 b 1 b 2 cos(kb 2 ) A 2 (k) = W * 1 (k)D l (k)W 1 (k) = - 1 9 -6 b 2 b 1 -6 b 1 b 2 + 2K 2 b 1 b 2 + -3 b 2 b 1 + 6 b 1 b 2 + K 2 b 1 b 2 cos(kb 2 ) = A 1 (k) A 3 (k) = W * 0 (k)D l (k)W 1 (k) = - 1 18 12 b 2 b 1 -6 b 1 b 2 + 2K 2 b 1 b 2 + 6 b 2 b 1 + 6 b 1 b 2 + K 2 b 1 b 2 cos(kb 2 ) = A 0 (k) ( 90
)
The eigensolutions of the spectral problem (52) are then determined and are given by

e iµ = 1 2A 3 -(A 1 + A 2 ) ± (A 1 + A 2 ) 2 -4A 2 0 (91) e -iµ = 1 2A 3 -(A 1 + A 2 ) ∓ (A 1 + A 2 ) 2 -4A 2 0 ( 92
)
The signs are selected as in section 2.3. As there is only one dof in this case, one has n + = 1 and after normalization one can choose q 1 (k) = 1. From relation (60), taking also p j (-k) = 1, yields

y(-k) = A 2 + e iµ A 3 1 = A 3 2A 3 -(A 1 + A 2 ) ± (A 1 + A 2 ) 2 -4A 2 0 + A 2 1 22 = ± 1 2 (A 1 + A 2 ) 2 -4A 2 0 1
this gives, with the notation of relation (70),

Q(-k) = ± 1 2 (A 1 + A 2 ) 2 -4A 2 0 F (-k) = 1 (93)
Near k = 0 (this means in fact near the normal incidence) one has the development

                         A 0 (k) = - b 2 b 1 - 1 6 K 2 b 1 b 2 + 1 36 6 b 2 b 1 + 6 b 1 b 2 + K 2 b 1 b 2 (kb 2 ) 2 + O((kb 2 ) 4 ) A 1 (k) = b 2 b 1 - 1 3 K 2 b 1 b 2 + 1 18 -3 b 2 b 1 + 6 b 1 b 2 + K 2 b 1 b 2 (kb 2 ) 2 + O((kb 2 ) 4 ) A 2 (k) = A 1 (k) A 3 (k) = A 0 (k) (94)
and this leads to

Q(-k) = ±iKb 2 1 - 1 12 (Kb 1 ) 2 ×   1 - 1 2 1 + (Kb 2 ) 2 3 -(Kb 1 ) 2 6 -(K 2 b 1 b 2 ) 2 36 1 -(Kb 1 ) 2 12 k 2 K 2   + O((kb 2 ) 4 ) (95)
From relation (58) one also has

f r (k) = W * 0 D l W 0 + e iµ W 1 = A 1 + e iµ A 3 = Q(-k) (96)
Following the rule that the positive waves are such that Re iωq H j f r j > 0, one has to choose the minus sign in relation (93). For the case Kb 1 ≪ 1 and k ≪ K, one has the approximation:

f (0) ≃ -iKb 2 (97)
The power across the boundary is thus

P = 1 2 Re(f lef t→right v * ) = 1 2 Re(f r (0)(-iω) * ) = 1 2 Kb 2 ω > 0 (98)
For the second order approximation one has

G 0 = -(F -1 Q)(0) = iKb 2 1 - 1 12 (Kb 1 ) 2 ≈ iKb 2 G 1 = i(F -1 Q) ′ (0) = 0 G 2 = (F -1 Q) ′′ (0) = ib 2 K 1 + (Kb 2 ) 2 3 -(Kb 1 ) 2 6 -(K 2 b 1 b 2 ) 2 36 1 -1 12 (Kb 1 ) 2 ≈ ib 2 K (99)
The relation between forces and displacements dofs on the boundary of the element is thus given by using ( 78)

f (x 1 , x 2 ) ≈ iKb 2 q(x 1 , x 2 ) + i 2Kb 2 (q(x 1 , x 2 + b 2 ) + q(x 1 , x 2 -b 2 ) -2q(x 1 , x 2 )) (100)
At order 0 one finds the classical approximation of the radiating boundary condition. The factor b 2 is present because the force is calculated over an edge of an element of length b 2 .

We compare four solutions in the following [START_REF] Brebbia | Boundary elements techniques in engineering[END_REF] The zero order solution with the numerical computation of Q(0) leading to the relation

f l n = G 0 p l a (2)
The zero order solution with the simplified computation of Q(0) leading to the relation

f l n = iKb 2 p l a (3)
The second order solution with the numerical computation of

G 2 = (F -1 Q) ′′ (0) ≈ ((F -1 Q)(δ) + (F -1 Q)(-δ) -2(F -1 Q)(0))/δ 2 leading to the relation f l n = G 0 p l a + G 2 2b 2 2 (p l+1 a + p l-1 a -2p l a ) (4)
The second order solution with the simplify computation of Q(0) and (F -1 Q) ′′ (0) leading to the relation

f l n = iKb 2 p l a + i 2Kb 2 (p l+1 a + p l-1 a -2p l a )

Example

Here we compute the error of relation (86) for different cases as shown in figure 5. Case a) is for a sound pressure created by a plane wave while in case b) the sound pressure is created by a point source. The acoustic element used for the computation of the boundary condition can be of size b 1 × b 2 = 0.01m × 0.01m or b 1 ×b 2 = 0.05m×0.05m. The sound velocity is c = 340m/s and the distance between the points is h = b 2 .

The first example is for a pressure created by a plane wave at the incidence angle θ = 10 o . The error for the four cases listed in the precedent section are plotted in figure 6. The error is the same for any point on the vertical axis. It can be observed that the second order relations are much better than the first ones as expected. The comparison of the two sizes for the acoustic element shows that the size 0.05m × 0.05m can reduce the accuracy of the solution for high frequencies and the second order condition. In these cases it is better to use elements with small sizes.

In figure 7 the error is plotted versus the angle of incidence of the plane wave.

An acoustic element of size b 1 × b 2 = 0.01m × 0.01m has been used. The solution is accurate (error less than 1%) for angles up to 10 o for a zero order condition and up to 30 o for a second order condition. This clearly shows that second order conditions are much better for waves at oblique incidence.

Finally the error is plotted versus the distance along the y axis in figure 8 for a pressure created by a point source at distance D = 1m on the x axis and at the frequency 1000Hz. The reduction in accuracy can be observed as we move along the y axis leading to greater incidence angles in agreement with the precedent observation on the plane waves. All these points confirm the accuracy of the method proposed here.

Finite element examples

Acoustics

In this section we use the precedent boundary condition to solve some finite element problems for different frequencies and mesh densities. We consider first a finite element acoustic problem on a square domain with a point source excitation in its center. A two dimensional domain of size 1m×1m is generated by Ansys. The domain and an example of mesh are presented in figure 9. The size of the acoustic element is b = 0.025m leading to 40 × 40 elements for the whole domain. In each case, only square elements are used for the mesh. The sound velocity is c = 340m/s. Only mesh information, stiffness and mass matrices are picked out and are then introduced into Matlab to get the results presented below. The procedure is done over the frequency band [0, 2000Hz]. Numerical Green's functions are calculated for zero and second order boundary conditions. The excitation is at point (0, 0) and the analytical solution in infinite space is given by formula (84). The Green's functions are presented in figure 10 for a point at (0.3, 0) on the horizontal axis and in figure 11 for a point at (0.3, 0.3) along the diagonal. Good agreements between the two types of absorbing boundary conditions and the analytical solution can be observed. Both boundary conditions fail at low frequencies because the size of the domain is too small compared to the wavelength. Similarly the error for high frequencies are of same orders for both boundary conditions. For intermediate frequencies the error is lower for the second order boundary condition. This is more clearly seen in figure 12 where the relative error for the point (0.3, 0.3) is plotted versus the frequency.

The same results are presented in figure 13 for the point (0.3, 0.3) and a mesh density of 80 × 80 elements. The solution is clearly much better at high frequencies meaning that the errors seen in figures 10 and 11 can be explained by elements too large for these frequencies and not by the quality of the boundary conditions. In figure 14 the domain is now 2m × 2m with 80 × 80 elements, so with elements of the same size as for figures 10 and 11. Results are plotted for the point (0.3, 0.3). Now the improvement is clearly seen for low frequencies while there is no difference for high frequencies.

Two dimensional elastodynamics

The analytical solution in direction e n of the two dimensional elastodynamics case when submitted to a unit force at origin in direction e m , is given by

G nm (r) = i 4µ Aδ nm + B x n x m r 2 (101)
with:

A = H 0 (K T r) - 1 K T r [H 1 (K T r) -βH 1 (K L r)] B = -2A + H 0 (K T r) + β 2 H 0 (K L r)
where H 0 and H 1 are the Hankel functions of first type, of orders zero and one respectively. The wavenumbers are K L = ω/c L and K T = ω/c T for the longitudinal and transverse waves respectively. The velocities c L , c T and the ratio between them β are given by,

                 β = c T c L c 2 L = λ + 2µ ρ c 2 T = µ ρ with λ = Eν (1 + ν)(1 -2ν) and µ = E 2(1 + ν) (102)
In this example, the same global meshes as for the acoustic case are used. The boundary condition is computed from the square four nodes elements. The material is steel with E = 2.10 11 P a, ν = 0.3 and ρ = 7800kg/m 3 and plane strain conditions are used in the computation. The sizes of the periodic cell can be b = 0.025m, 0.0125m or 0.00625m.

In figure 15, numerical solutions are compared with analytical solutions for the point (0.5, 0) and (0, 0.5) for different sizes of the cell. The curves represent the real and imaginary parts of the first component of the displacement for an excitation at origin in direction 1. The same remarks as the previous examples can be made: 1) a denser mesh leads to lower errors over the high frequency band [8.10 3 Hz, 20.10 3 Hz]; 2) for the low frequency band [0, 8.10 3 Hz], numerical results are different from the analytical solutions due to the finite size of the domain (L = 1m).

In figure 16, the results for these two points are presented when the size of the domain is increased successively with L = 1m, 2m and 4m. In this case, when the size of the cell is fixed to b = 0.025m, a larger domain leads to lower errors over the low frequency band [0, 8.10 3 Hz]. The results are the same as previously in the high frequency band because in this domain the precision depends on the size of the elements and not on the size of the global domain. Some improvements are however seen for intermediate frequencies.

In figure 17 the error for boundary conditions of zero and second orders are plotted versus the frequency. The second order condition is much more accurate for intermediate frequencies as for the acoustic case.

Conclusion

In this paper, a method to determine absorbing boundary conditions for two dimensional periodic media has been presented. It works directly on the discretized equations. The boundary condition is first obtained as a global impedance relation and is then localized into boundary conditions of various orders. In two examples, good agreements were observed when compared with analytical solutions.

In any case, the proposed method is efficient because it requires only the discrete dynamic matrices which can be obtained by any standard finite element software. This method could be used for media with more complex behaviors than those presented in the precedent examples.

Fig. 2 .Fig. 5 .

 25 Fig. 1. Computational domain

Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 6789 Fig. 6. Error versus the frequency for a plane wave at 10 o for (a) an element size 0.01m × 0.01m and (b) an element size 0.05m × 0.05m

Fig. 10 .

 10 Fig. 10. Comparison of analytical and numerical Green's functions, with the present method at order 0 -and at order 2 . . for the 2D acoustics at point (0.3, 0) with L = 1m and b = 0.025m: a) real part, b) imaginary part.

  Green's function at (0.3,0.3), Imaginary part

Fig. 11 .

 11 Fig. 11. Comparison of analytical and numerical Green's functions, with the present method at order 0 -and at order 2 . . for the 2D acoustics at point (0.3, 0.3) with L = 1m and b = 0.025m: a) real part, b) imaginary part.

Fig. 14 .

 14 Fig. 14. Comparison of analytical and numerical Green's functions, with the present method at order 0 -and at order 2 . . for the 2D acoustics at point (0.3, 0.3) with L = 2m and b = 0.025m: a) real part, b) imaginary part.

Fig. 15 .Fig. 16 .Fig. 17 .

 151617 Fig. 15. Comparison of analytical and numerical Green's functions for the 2D elasticity with different sizes of elements: a) Real part at (0, 0.5), b) Imaginary part at (0, 0.5), c) Real part at (0.5, 0), d) Imaginary part at (0.5, 0).

  2µ d 12 + (d 11 + d 22 )e µ + d

						21 = 0	(7)
	where d ij are the components of the matrix d. Taking into account the sym-
	metries in the matrix d, this yields
	e 2µ + 2	d 11 d 12	e µ + 1 = 0	(8)
	whose solutions are
	e µ ± = -	d 11 d 12	± (	d 11 d 12	) 2 -1


Yoann Dieudonné

Franck Petit
email: franck.petit@lip6.fr

Vincent Villain

Leader Election Problem Versus Pattern Formation Problem

Keywords: Mobile Robot Networks, Pattern Formation Problem, Leader Election Problem

Leader election and arbitrary pattern formation are fundammental tasks for a set of autonomous mobile robots. The former consists in distinguishing a unique robot, called the leader. The latter aims in arranging the robots in the plane to form any given pattern. The solvability of both these tasks turns out to be necessary in order to achieve more complex tasks.

In this paper, we study the relationship between these two tasks in a model, called CORDA, wherein the robots are weak in several aspects. In particular, they are fully asynchronous and they have no direct means of communication. They cannot remember any previous observation nor computation performed in any previous step. Such robots are said to be oblivious. The robots are also uniform and anonymous, i.e, they all have the same program using no global parameter (such that an identity) allowing to differentiate any of them. Moreover, none of them share any kind of common coordinate mechanism or common sense of direction, except that they agree on a common handedness (chirality).

In such a system, Flochini et al. proved in [9] that it is possible to solve the leader election problem for n ≥ 3 robots if the arbitrary pattern formation is solvable for n ≥ 3. In this paper, we show that the converse is true for n ≥ 4 and thus, we deduce that both problems are equivalent for n ≥ 4 in CORDA provided the robots share the same chirality. The possible equivalence for n = 3 remains an open problem in CORDA.

Introduction

Mobile robots working together to perform cooperative tasks in a given environment is an important, open area of research. Teams (or, swarms) of mobile robots provide the ability to measure properties, collect information and act in a given physical environment. Numerous potential applications exist for such multi-robot systems, to name only a very few: environmental monitoring, large-scale construction, risky area surrounding or surveillance, and exploration of awkward environments.

In a given environment, the ability for the swarm of robots to succeed in the accomplishment of the assigned task greatly depends on [START_REF] Canepa | Stabilizing flocking via leader election in robot networks[END_REF] global properties assigned to the swarm, and (2) individual capabilities each robot has. Examples of such global properties are the ability to distinguish among themselves at least one (or, more) robots (leader), to agree on a common global direction (sense of direction), or to agree on a common handedness (chirality). The individal capacities of a robot are its moving capacities and its sensory organs.

To deal with cost, flexibility, resilience to dysfunction, and autonomy, many problems arise for handling the distributed coordination of swarms of robots in a deterministic manner. This issue was first studied in [START_REF] Suzuki | Agreement on a common x-y coordinate system by a group of mobile robots[END_REF][START_REF] Suzuki | Distributed anonymous mobile robots -formation of geometric patterns[END_REF], mainly motivated by the minimal level of ability the robots are required to have in the accomplishment of basic cooperative tasks. In other words, the faisibility of some given tasks is addressed assuming swarm of autonomous robots either devoid or not of capabilities like (observable) identifiers, direct means of communication, means of storing previous observations, sense of direction, chirality, etc. So far, except the "classical" leader election problem [START_REF] Canepa | Stabilizing flocking via leader election in robot networks[END_REF][START_REF] Dieudonné | Deterministic leader election in anonymous sensor networks without common coordinated system[END_REF][START_REF] Flocchini | Arbitrary pattern formation by asynchronous, anonymous, oblivious robots[END_REF][START_REF] Prencipe | Distributed coordination of a set of autonomous mobile robots[END_REF], most of the studied tasks are geometric problems, so that pattern formation, line formation, gathering, and circle formation-refer to [START_REF] Cohen | Local spreading algorithms for autonomous robot systems[END_REF][START_REF] Defago | Circle formation for oblivious anonymous robots with no common sense of orientation[END_REF][START_REF] Dieudonné | Circle formation of weak mobile robots[END_REF][START_REF] Flocchini | Distributed coordination of a set of autonomous mobile robots[END_REF][START_REF] Flocchini | Arbitrary pattern formation by asynchronous, anonymous, oblivious robots[END_REF][START_REF] Katreniak | Biangular circle formation by asynchronous mobile robots[END_REF][START_REF] Suzuki | Distributed anonymous mobile robots -formation of geometric patterns[END_REF] for these problems.

In this paper, we concentrate on two of the aforementioned problems: leader election and pattern formation. The former consists in moving the system from an initial configuration were all entities are in the same state into a final configuration were all entities are in the same state, except one, the leader. The latter consists in the design of protocols allowing autonomous mobile robots to form any (arbitrary) geometric pattern.

The issue of whether the pattern formation problem can be solved or not according to some capabilities of the robots is addressed in [START_REF] Flocchini | Arbitrary pattern formation by asynchronous, anonymous, oblivious robots[END_REF]. Assuming that every robot is able to observe all its pairs, the authors consider sense of direction and chirality. They show by providing an algorithm that, if the robots have sense of direction and chirality, then they can form any arbitrary pattern. They rafine their result by showing that with the lack of chirality (i.e., assuming that they have sense of direction only), the problem can be solved in general with an odd number of robots only. They also show that, assuming robots having no sense of direction, then, in general, the robots cannot form an arbitrary pattern, even with chirality. As a matter of fact, the idea of proof relies on the fact that if it is possible to solve the pattern formation problem, then the robots can form an asymmetric configuration in order to distinguish a unique robot. That means that, the ability to (deterministically) form a particular type of patterns implies the ability to (deterministically) elect a robot in the system as the leader. In other words, if it is not possible to solve the leader election problem, then it is not possible to solve the pattern formation problem. However, in [START_REF] Dieudonné | Deterministic leader election in anonymous sensor networks without common coordinated system[END_REF], assuming anonymous robots (possibly motionless) with chirality only (without sense of direction), the authors provide a complete characterization (necessary and sufficient conditions) on the robots positions to deterministically elect a leader. An interesting question arises from the above facts: "With robots devoid of sense of direction, does the (arbitrary) pattern formation problem becomes solvable if the robots have the possibility to distinguish a unique leader?" In [START_REF] Yamashita | Characterizing geometric patterns formable by oblivious anonymous mobile robots[END_REF], the authors provide a positive answer to this question assuming that robots have the chirality property. There result holds in the semi-synchronous model (SSM) introduced in [START_REF] Suzuki | Agreement on a common x-y coordinate system by a group of mobile robots[END_REF], a.k.a. Model SYm in the literature. In this paper, we show that this results also holds for n ≥ 4 robots in a fully asynchronous model, called CORDA.

Combined with the result in [START_REF] Flocchini | Arbitrary pattern formation by asynchronous, anonymous, oblivious robots[END_REF] and provided the robots have the chirality property, we deduce that Leader Election and Pattern Formation are two equivalent problems in CORDA for n ≥ 4 robots, in the precise sense that, the former problem is solvable if and only if the latter problem is solvable. The possible equivalence for the case n = 3 remains an open problem in CORDA.

The rest of the paper is organized as follows: In Section 2, we describe the distributed systems and state the problems considered in this paper. The proof of equivalence is given in Section 3 for any n ≥ 4 by providing an algorithm working in CORDA. (Due to lack of space, technical proofs have been moved in the appendix.) Finally, we make concluding remarks in Section 4.

Preliminaries

In this section, we define the distributed system and the problems considered in this paper.

Distributed Model.

We adopt the model CORDA introduced in [START_REF] Prencipe | Instantaneous actions vs. full asynchronicity : Controlling and coordinating a set of autonomous mobile robots[END_REF]. The distributed system considered in this paper consists of n robots r 1 , r 2 , • • • , r n -the subscripts 1, . . . , n are used for notational purpose only. Each robot r i is viewed as a point in a two-dimensional space unbounded and devoid of any landmark. When no ambiguity arises, r i also denotes the position in the plane occupied by that robot. Each robot has its own local coordinate system and unit measure. The robots do not agree on the orientation of the axes of their local coordinate system, nor on the unit measure.

Definition 1 (Sense of Direction) A set of n robots has sense of direction if the n robots agree on a common direction of one axis (x or y) and its orientation. The sense of direction is said to be partial if the agreement relates to the direction only -ie. they are not required to agree on the orientation.

In Figure 1, the robots have sense of direction in the cases (a) and (b), whereas they have no sense of direction in the cases (c) and (d).

Given an x-y Cartesian coordinate system, the handedness is the way in which the orientation of the y axis (respectively, the x axis) is inferred according to the orientation of the x axis (resp., the y axis).

Definition 2 (Chirality) A set of n robots has chirality if the n robots share the same handedness. In Figure 1, the robots have chirality in the cases (a) and (c), whereas they have no chirality in the cases (b) and (d). In the sequel, we assume that all the robots have the chirality property.

The robot's life is viewed as an infinite sequence of cycles. Each cycle is a sequence of four states Wait-Observe-Compute-Move characterized as follows.

Life cycle. Initially, a robot is in the waiting state (Wait). Asynchronously and independently from other robots, it observes its surroundings (Observe) by using its sensors. The latter ones return a set of all the positions occupied by at least one robot, with respect to its own coordinate system. Then, from its new observations the robot computes its next location (Compute) according to a given protocol which is the same one for all the robots. Once the computation is done, the robot moves towards its new location (Move). However, it is assumed the distance travelled by a robot in a cycle is bounded by a constant σ. So, if the destination point is too far, the robot may stop its motion before reaching the computed location . Finally, the robot returns to the waiting state.It is assumed that the amount of time spent in each phase of a cycle is finite but unpredictable and may be different for each cycle and for each robot. That is why the robots are considered to be fully asynchronous.

Finally we assume that the robots are uniform and anonymous, i.e, they all have the same program using no local parameter (such that an identity) allowing to differentiate any of them. Moreover, they have no direct means of communication and they are oblivious, i.e., none of them can remember any previous observation nor computation performed in any previous cycles.

Leader Election Problem

The leader election problem considered in this paper is stated as follows: Given the positions of n robots in the plane, the n robots are able to deterministically agree on the same robot L called the leader. Initially, the robots are in arbitrary positions, with the only requirement that no two robots are in the same position.

Arbitrary Pattern Formation Problem

In the Arbitrary Pattern Formation Problem, the robots have in input the same pattern, called the target pattern P, described as a set of positions in the plane given in lexicographic order (each robot sees the same pattern according to the direction and orientation of its local coordinate system). They are required to form the pattern: at the end of the computation, the positions of the robots coincide, in everybody's local view, with the positions of P, where P may be translated, rotated, and scaled in each local coordinate system. Initially, the robots are in arbitrary positions, with the only requirement that no two robots are in the same position, and that, of course, the number of positions prescribed in the pattern and the number of robots are the same.

Equivalence for n ≥ 4

In this section we prove the following theorem: Theorem 3 In CORDA, assuming a cohort of n ≥ 4 robots having chirality and devoid of any kind of sense of direction, if the leader election problem is solvable, then the pattern formation problem is solvable.

In a recent paper [START_REF] Flocchini | Arbitrary pattern formation by asynchronous, anonymous, oblivious robots[END_REF], the authors prove the following result:

Theorem 4 [START_REF] Flocchini | Arbitrary pattern formation by asynchronous, anonymous, oblivious robots[END_REF] In CORDA, assuming a cohort of n ≥ 3 robots devoid of any kind of sense of direction, if it is possible to solve the pattern formation problem, then the leader election problem is solvable too.

So, to prove Theorem 3, from Theorem 4, it remains to show the following lemma: Lemma 5 In CORDA, assuming a cohort of n ≥ 4 robots having chirality and devoid of any kind of sense of direction, if the leader election problem is solvable, then the pattern formation problem is solvable.

The proof of Lemma 5 is mainly based on the existence of a protocol allowing to form an arbitrary target pattern if initially the robots are in a configuration allowing the robots to deterministically elect a leader. Such a configuration is called a leader configuration.

Definition 6 (Leader configuration) A configuration allowing the robots to deterministically elect a leader is called a leader configuration.

Definitions and Basic Properties

In the rest of this paper, we assume the set of all the positions Q occupied by the robots in the plane is the set of all the coordinates expressed in a cartesian coordinate system S which is unknown for all the robots. However, all the coordinates Q expressed in S coincide with all the cordinates Q expressed in everybody's local system where Q may be translated, rotated or scaled.

Definition 7 (Smallest enclosing circle) [4] Given a set Q of n ≥ 2 positions p 1 , p 2 , • • • , p n on the plane, the smallest enclosing circle of Q , called SEC(Q), is the smallest circle enclosing all the positions in Q.
When no ambiguity arises, SEC(Q) is shortly denoted by SEC and SEC(Q) ∩ Q indicates the set of all the positions both on SEC(Q) and Q. Besides, we say that a robot r is inside SEC if, and only if, r is not located on the circumference of SEC. In any configuration Q, SEC is unique and can be computed in linear time [START_REF] Megiddo | Linear-time algorithms for linear programming in r3 and related problems[END_REF]. Note that since the robots have the ability of chirality, they are able to agree on a common orientation of SEC, denoted , in the sequel referred to as the clockwize direction.

The following property contains some details about the smallest enclosing circle:

Property 8 [START_REF] Welzl | Smallest enclosing disks (balls and ellipsoids)[END_REF] SEC passes either through two of the positions that are on the same diameter (opposite positions), or through at least three positions. SEC does not change by eliminating or adding positions that are inside it. SEC does not change by adding positions on its boundary. However, it may be possible that SEC changes by either eliminating or moving positions on its circumference.

Examples showing the latter assertion of Property 8 are proposed in Figure 2.

Definition 9 (Critical position) [START_REF] Katreniak | Biangular circle formation by asynchronous mobile robots[END_REF] Given a set Q of distinct positions. We say that a position An example of such a critical robot is given by Figure 2, Case (a). According to Property 8, a critical position cannot be inside SEC. So, we have the following corollary:

p is critical iff SEC(Q) = SEC(Q \ {p}).
Corollary 10 Let Q be a configuration. If there exists a critical position p in Q, then p is on the circumference of SEC(Q).

Before giving other properties about critical positions, we need to define extra notions.

Definition 11 (adjacent(r, C, )) Given a circle C and a group of robots located on it, we say that r ′ = adjacent(r, C, ) if r ′ is the next robot on C just after r in the clockwise direction.

In the same way, we can define adjacent(r, ) in the counterclockwise direction. When no ambiguity arises, adjacent(r, C, ) is shortly denoted by adjacent(r, ). Sometimes, if r ′ = adjacent(r, ), we simply say that r ′ and r are adjacent.

Definition 12 (angle(p, c, p ′ , )) Given a circle C centered at c and two points p and p ′ located on it, angle(p, c, p ′ , ) is the angle centered at c from p to p ′ in the clockwise direction.

In the same way, we can define angle(p, c, p ′ , ) in the counterclockwise direction. The following properties are fundammental results about smallest enclosing circles:

Lemma 13 [START_REF] Cieliebak | Gathering autonomous mobile robots[END_REF] Let r i , r j and r k be three consecutive robots on SEC centered at c such that r j = adjacent(r i , ) and r k = adjacent(r j , ). If angle(r i , c, r k , ) ≤ 180 o , then r j is non-critical and SEC does not change by eliminating r j .

Corollary 14 Let SEC(Q) be the smallest circle enclosing all the positions in Q. For all couple of positions r i and r j in SEC(Q) ∩ Q such that r j = adjacent(r i , ), we have angle(r i , c, r j , ) ≤ 180 o .

Lemma 15 [START_REF] Cieliebak | Gathering autonomous mobile robots[END_REF] Given a smallest enclosing circle with at least four robots on it. There exists at least one robot which is not critical. Definition 16 (Concentric Enclosing Circle) Given a set P of distinct positions. We say that C P is a concentric enclosing circle if and only if it is centered at the center c of SEC, has a radius strictly greater than zero and it passes through at least one position in P .

In the following, SC P and |SC P | respectively denote the set of all the concentric enclosing circle in P and its cardinality. For some k such that 1 ≤ k ≤ |SC P |, C P k indicates the k th greatest concentric enclosing circle in P and k i=1 C P i is the set of the k first greatest enclosing circles in P . Moreover, we assume that a position (or robot) located inside a concentric enclosing circle C P k is not on the circumference of C P k . C P i ∩ P indicate the set of all the positions both on C P i and P .

Remark 17 From Definition 16, SEC is the greatest concentric enclosing circle of SC (i.e., SEC = C 1 ) and the center of SEC cannot be a concentric enclosing circle.

From Definition 16, we can introduce the notion of agreement configuration:

Definition 18 (Agreement Configuration) A configuration Q is an agreement configuration if, and only if both following conditions hold:

1. There exists a robot r l in Q such that r l is the unique robot located on the smallest concentric enclosing circle

C Q |SC Q | , 2.
There is no robot at the center of SEC(Q).

In an agreement configuration, r l is called the leader Definition 19 (Equivalent agreement configuration) Two agreement configuration Q 1 and Q 2 is said to be equivalent if, and only if, both following conditions hold: 1. SEC(Q 1 ) and SEC(Q 2 ) are superimposed. 2. Let c 1 and c 2 be respectively the center of SEC(Q 1 ) and the center of SEC(Q 2 ). Let r l1 and r l2 be respectively the leader in Q 1 and the leader in Q 2 . [c 1 , r l1 ) and [c 2 , r l2 ) are superimposed.

We end with the three following definitions:

Definition 20 (M ap(Q, P)) Let Q and P be respectively an agreement configuration formed by the robots in the plane and a target pattern. M ap(Q, P) is the set of all the final positions P expressed in the plane where the robots currently lies and computed as follows: 1. First, the center of SEC(P)) is translated to the center of SEC(Q). 2. Then, let o, c, r l and s be respectively the center of SEC(Q), the center of SEC(P), the leader in Q and the first non-critical position (in the lexicographic order) located on the smallest concentric enclosing circle of P. P is rotated so that the half-line [o, r l ) is viewed as the half-line [c, s). 3. Finally, P is scaled with respect to the radius of SEC(Q) in order that all the distances are expressed according to the radius of SEC(Q). In particular SEC(Q) = SEC(P).

An example showing the construction of Definition 20 is given in Figure 3.

Definition 21 ((k, P)-partial pattern) Let Q and P be respectively an agreement configuration formed by the robots in the plane and a target pattern. We say that: In the sequel, we say that Q is a maximal (k, P)-partial pattern if Q is a (k, P)-partial pattern and not a (k + 1, P)-partial pattern.

1. Q is a (0, P)-partial pattern if the leader in Q is inside the smallest concentric enclosing circle of M ap(Q, P). 2. Q is a (k, P)-partial pattern with 1 ≤ k ≤ M in(|SC Q |, |SC P |) if the three following properties holds: a. Q is a (0, P)-partial pattern. b. C M ap(Q,P) k ∩ M ap(Q, P) C Q k ∩ Q. c. k-1 i=1 C Q i ∩ Q = k-1 i=1 C M ap(Q,P) i ∩ M ap(Q, P).
Definition 22 (Extra robots) Let P and Q be respectively a target pattern and a configuration formed by the robots in the plane such that Q is a maximal (k, P)-partial pattern. We say that a robot r is an extra robot if one of the two following properties holds:

1. k = 0, r is inside SEC(Q), and r is not the leader in Q;

The protocol

Starting from a leader configuration, the protocol, shown in Algorithm 1, allows to form any target pattern P. It is a compound of two procedures presented in the two following subsections: 1. Protocol < Leader; Agreement> transforms an arbitrary leader configuration into an agreement configuration. 2. Protocol < Agreement; Pattern> transforms an agreement configuration into a pattern P.

Algorithm 1 Form an arbitrary pattern starting from a leader configuration (n ≥ 4). In a leader configuration, we have the following corollary:

Corollary 23 [START_REF] Dieudonné | Deterministic leader election in anonymous sensor networks without common coordinated system[END_REF] If the robots are in a leader configuration, then they can distinguish a unique leader which is one of the closest robot to the center of the smallest enclosing circle of the configuration, provided that they share the property of chirality.

So, from Corollary 23, we know that we can distinguish a unique robot r l , called the leader, which is one of the robots closest to the center c of SEC(Q). However, according to Definition 18, if r l is at the center of SEC(Q) or if r l is not the unique robot closest to the center of SEC(Q), Q is not an agreement configuration. In that case, Procedure < Leader ; Agreement > allows to transform the leader configuration into an agreement configuration. Algorithm 2 describes Procedure < Leader; Agreement>.

In Algorithm 2, we use two subsoutines: Leader(Q) and M oveT o(p, →). The former returns the unique leader from a leader configuration Q. The latter allows a robot r to move towards the point p, using a straight movement. Remark 26 From Definition 24, we know that p i is not located on P -arc(p i , p i+1 , C, ).

In the remainder, we say that a P -arc is free if there is no robot located on it. In Figure 4, the circles denote the positions to achive. The crosses depict the robots. The P -arc starting after f (f excluded) and finishing at a is free. Definition 27 (Deadlock Chain) A Deadlock Chain is a consecutive sequence of P -arc starting from a free P -arc P 0 and followed in the counterclockwise direction by a P -arc P 1 such that: 1. P 1 is a P -arc(p, p ′ , C, ) such that angle(p, c, p ′ , ) = 180 o and there is only one robot r on it and r is located at p ′ , 2. and P 1 is followed in counterclockwise by a consecutive sequence (possibly empty) of P -arc(p, p ′ , C, ) such that there is only one robot r on each of them and r is located at p ′ , and that consecutive sequence (possibly empty) is followed by a P -arc(p, p ′ , C, ) such that there is at least two robots on it and one of them is located at p ′ . This P -arc is called the last P -arc of the deadlock chain.

In Figure 4, the segment staring from Position a (a included) to Position b (b excluded) forms a deadlock chain.

Definition 28 (Deadlock Breaker) Let P -arc(p, p ′ , C, ) be the last P -arc of a deadlock chain. The deadlock breaker is the robot located at p ′ .

In Figure 4, the robot located at Position c is the deadlock breaker.

Sketch of Correctness Proof of Algorithm 1

We first show that by executing Algorithm 1, the smallest enclosing circle SEC(Q) remains invariant-Lemma 30. Next, we prove that if the robots form a leader configuration which is not a final pattern P and not an agreement configuration, they eventually form an agreement configuration-Lemma 31.

In a future work, we would like to investigate the same equivalence (1) for the case n = 3 and (2) in the case without chirality.

Missing Proofs.

Lemma 30 According to Algorithm 1, the smallest enclosing circle SEC(Q) remains invariant.

Proof. Assume by contradiction SEC(Q) does not remain invariant. From Corollaries 10 and 14 and Property 8, we deduce that can occurs if, and only if:

• Either a robot r moves outside SEC(Q). However, according to Algorithm 1, no robot moves outside SEC(Q). That is a contradiction.

• Or an angle strictly greater than 180 o appears between two adjacent robots r i-1 and r i , i.e., angle(r i-1 , c, r i , ) > 180 o with r i-1 = adjacent(r i , SEC, ). This subcase can occur if, and only if

-Either a critical robot leaves SEC(Q). However, according to Algorithm 2 no critical robot leaves SEC(Q) ( only the first non-critical robot on SEC in clockwise is sometimes allowed to move). Furthermore, according to Algorithm 3 some robots are allowed to leave SEC(Q) only if these latter ones are extra robots and Q is a (1, P)-partial pattern. That implies some robots are allowed to leave SEC(Q) = C Q 1 only if these latter ones do not occupy a position ∈ M ap(Q, P) ∩ SEC(Q) and all the positions in M ap(Q, P) ∩ SEC(Q) are occupied by some robots. However, from Corollary 14 we know that for all couple of positions r i and r j on M ap(Q, P) ∩ SEC(Q) such that r j = adjacent(r i , ), we have angle(r i , c, r j , ) ≤ 180 o . Consequently, when extra robots leaves SEC(Q), SEC(Q) is not changed. So, no critical robot leaves SEC(Q). That is a contradiction.

-Or two adjacent robots r i-1 and r i , such that r i-1 = adjacent(r i , SEC, ), move along SEC(Q) so that angle(r i-1 , c, r i , ) > 180 o . That might occur only by applying Algorithm 4. However, if r i is allowed to move, it can only move in clockwise towards a position p such that angle(r i-1 , c, p, ) ≤ 180 o . Furthermore, r i-1 is never allowed to move in counterclockwise. So, angle(r i-1 , c, r i , ) is always less than or equal to 180 0 . That is a contradiction.

2

Lemma 31 If the robots form a leader configuration which is not a final pattern P and not an agreement configuration, they form an agreement configuration in a finite number of cycles.

Proof.

If the robots form a leader configuration which is not a final pattern P and not an agreement configuration, then from Corollary 23 we have two cases to consider: either (1) the leader r l is at the center c of SEC or (2) r l is not the unique robot closest to c.

• Case 1. r l is at the center of SEC. According to Procedure < Leader ; Agreement >, r l moves away from c towards a position which is closer to the center than the second robot closer to the center. Furthermore, from Lemma 30, the center c of SEC remains invariant even if r l moves. So, r l remains the unique leader and, by fairness, we deduce that an agreement configuration is formed in a finite number of cycles.

• Case 2. r l is not the unique robot closest to c. In that case, we have two subcases to consider:

-Case 2.1. r l is not a critical robot. In this subcase, r l moves towards a position which is located between c and itself (except c and itself). From Lemma 30, the center c of SEC remains invariant even if r l moves. So, by fairness we know that an agreement configuration is formed in a finite number of cycles.

-Case 2.2. r l is a critical robot. From Corollary 10, r l is on the circumference of SEC. However, by assumption r l is also one of the robots closest to the center of SEC. So, we deduce that all the robots are on SEC. Hence, by Lemma 15, we deduce there is at least one non-critical robot on SEC because there are at least four robots on it (recall that we assume the number of robot is greater than or equal to 4). According to Procedure < Leader ; Agreement > the first non-critical robot r k starting from r l on SEC in clockwise is allowed to move toward a position located between itself and c (except c and itself). From Lemma 30, the center c of SEC remains invariant even if r k moves. So, by fairness r k becomes the unique robot closest to c and it is not located at c.. So, the robots form an agreement configuration in a finite number of cycles.

2

Lemma 32 Starting from an agreement configuration Q, the robots remains in an equivalent agreement configuration or the target pattern P is formed in a finite number of cycles.

Proof.

According to Lemma 30, SEC(Q) and its center c remain invariant. Moreover, according to Algorithm 1 and more precisely Algorithm 3 no robot is allowed to pass r l .

So, if r l is not allowed to move then, according to Definition 19 all the robots remain in an equivalent agreement configuration.

If r l is allowed to move then, according to Algorithm 3 that can occur only in three cases:

• Case 1. The robots do not form any (k, P)-partial pattern. In that case, r l moves in straight line towards the middle p of the segment [c, s] in order to get closer to the center c. However, from Definition 20, we know that s is on the half line [c, r l ). So, during the motion of r l , all the robots clearly remain in an equivalent agreement configuration.

• Case 2. The center c of SEC(Q) is in M ap(Q, P) and all the positions in M ap(Q, P) are occupied except c. In that case, r l chooses to move towards c in straight line (i.e., along [c, r l )) in order to occupy the last free position in M ap(Q, P). Until r l has not reached c, the robots remain in an equivalent agreement configuration because r l is still on the same half line [c, r l ) and it remains the unique robot closest to c. So by fairness, it reaches c in a finite number of cycle and the pattern P is formed.

• Case 3. The center c of SEC(Q) is not in M ap(Q, P) and all the positions in M ap(Q, P) are occupied except the first non critical position s located on the smallest concentric enclosing circle. In that case, r l chooses to move towards s in order to occupy the last free position in M ap(Q, P). From Definition 20, we know that s is on the half line [c, r l ) and thus, until r l has not reached s the robots remain in an equivalent agreement configuration because r l is still on the same half line [c, r l ) and it remains the unique robot closest to c. So by fairness, it reaches s in a finite number of cycle and the pattern P is formed.

2

Figure 1 :

 1 Figure 1: Four examples showing the relationship between Sense of Direction and Chirality

  Critical (white) robot cannot be deleted without changing SEC. SEC SEC (b) An example showing how SEC may change by moving one robot.

Figure 2 :

 2 Figure 2: Examples illustrating Property 8 .

  Position b is the first non-critical position s in P. P mapped on the actual agreement configuration Q formed by the robots (Robots are depicted as x's in the figure).

Figure 3 :

 3 Figure 3: An example showing a pattern P mapped on an agreement configuration Q-Definition 20.

2. k ≥ 1

 1 and (a) either r is inside the enclosing circle C M ap(Q,P) k and r is not the leader in Q; (b) or r is on the circumference of C M ap(Q,P) k and r does not occupy a position in C M ap(Q,P) k ∩ M ap(Q, P).

PAlgorithm 2

 2 := the target pattern ; if the robots do not form the target pattern then if the robots do not form an agreement configuration then Execute < Leader ; Agreement >; else Execute < Agreement ; Pattern >; 3.2.1 Procedure < Leader ; Agreement>. Procedure < Leader; Agreement> for any robot r i in an arbitrary leader configuration Q := the configuration where the robots currently lies; r l := Leader(Q); c := center of SEC(Q) if r l is located at c then r k := the closest robot to c ∈ Q \ {r l }; p := the middle of the segment [r l ; r k ]; if I am r l then M oveT o(p, →); endif else if r l is not critical then p := the middle of the segment [r l ; c]; if I am r l then M oveT o(p, →); endif else /* r l is critical and r l is on SEC*/ r k := the first non-critical robot starting from r l on SEC in clockwise. if I am r k then p := the middle of the segment [r k ; c]; M oveT o(p, →);

3. 2 . 2 and 1 .

 221 Procedure < Agreement; Pattern>. Procedure < Agreement; Pattern> is shown in Algorithm 3. The routine N earest extra robot(C M ap(Q,P) k+1 , Q, M ap(Q, P)) returns an extra robot r such that r is the closest extra robot to C M ap(Q,P) k+1 which is not located on C M ap(Q,P) k+1 . If several candidats exists, then the extra robots inside C M ap(Q,P) k+1 have priority. Finally, if there is again several candidats then these latter ones are located on the same concentric circle C centered at the center c of C M ap(Q,P) k+1 C is one of the concentric enclosing circle of M ap(P, Q) 2. p i and p i+1 belong to F inal P ositions 3. p i+1 = adjacent(p i , C, )

Figure 4 :

 4 Figure 4: An example showing a Deadlock Chain and a Deadlock Breaker.

Algorithm 3 Procedure < Agreement; Pattern> for any robot r i in an agreement configuration Q := the configuration where the robots currently lies; P := the target pattern; /* P is the same for all the robots */ r l := Leader(Q); s := the first non-critical position located on the smallest concentric enclosing circle of M ap(Q, P); if the robots do not form any (k, P)-partial pattern then /*r l is not inside the smallest concentric enclosing circle of M ap(Q, P) */ p := the middle of the segment [c; s]; if I am r l then M oveT o(p, →); endif else /* the robots form a (k, P)-partial pattern */ if the center of SEC(Q) ∈ M ap(Q, P) then x := the center of SEC(Q); else x := s; endif F inal P ositions := M ap(Q, P) \ {x}; if all the positions in F inal P ositions are occupied then if I am r l then M oveT o(x, →); endif else k := the maximal k for which Q is a (k, P)-partial pattern;

if there is at least one extra robot not located on C Starting from such a configuration, M ap(Q, P) remains invariant or the target pattern P is formed-Corollary 33. It follows that from an agreement configuration which is not a (k, P)-partial pattern, the robots eventually form a (0, P)-partial pattern-Lemma 34.

From this point on, note that according to Algorithm 3, F inal P ositions is equal to all the positions in M ap(P, Q) except: 1. either the center c of SEC(Q) if c ∈ M ap(P, Q), 2. or the first non critical position located on the smallest concentric enclosing circle of M ap(P, Q) if c / ∈ M ap(P, Q) Next, we show by induction that, from a configuration being a maximal (k, P)-partial pattern, the robots eventually form a (k +1, P)-partial pattern or the target pattern P is formed-Lemmas 35 to 37. From Lemma 37 and by induction we deduce the following theorem:

Theorem 29 Starting from a leader configuration, Algorithm 1 allows to solve the pattern formation problem in CORDA among a cohort of n ≥ 4 robots having chirality and devoid of any kind of sense direction.

Conclusion

We studied the relationship between the arbitrary formation problem and the leader election problem among robots having the chirality in CORDA. We gave an algorithm allowing to form an arbitrary pattern starting from any geometric configuration wherein the leader election is possible. Combined with the result in [START_REF] Flocchini | Arbitrary pattern formation by asynchronous, anonymous, oblivious robots[END_REF], we deduce that arbitrary pattern formation problem and Leader election are equivalent, i.e., it is possible to solve the pattern formation problem for n ≥ 4 if and only if the leader election is solvable too.

Corollary 33 From an agreement configuration, M ap(Q, P) remains invariant or the target pattern P is formed.

Lemma 34 From an agreement configuration which is not a (k, P)-partial pattern, the robots form a (0, P)-partial pattern in a finite number of cycles.

Proof.

From Definition 21, we know that if an agreement configuration is not a (0, P)-partial pattern then, the leader r l is not inside the smallest concentric enclosing circle of M ap(Q, P). From Corollary 33 and according to Algorithm 3, r l is inside the smallest concentric enclosing circle of M ap(Q, P) in a finite number of cycles.

2 then, all the positions in

are occupied in a finite number of cycles.

Proof.

If all the extra robots are on C M ap(Q,P) k+1

and there exists at least one position in

which is not occupied then the robots apply the routine

, F inal P ositions) (refer to Algorithm 4). Remark that by applying this routine, no robot can collide another robot since any robot can move only in clockwise and any move of a robot on C M ap(Q,P) k+1 is only allowed in arc of circle containing no robot. Moreover, since k ≥ 1,

= SEC(Q) and thus, it is no need to prevent from creating an angle strictly greater than 180 o between two adjacent robots. In the remainder of this proof, we denote by α the number of extra robots located on C M ap(Q,P) k+1

, β the number of P -arc on C M ap(Q,P) k+1 and γ the number of free

. According to Algorithm 2, the acute reader noticed that the number α of extra robots is greater than or equal to the number β of P -arc on C M ap(Q,P) k+1

. We consider two cases.

• All the P -arcs are not free. According to Algorithm 4, each last robot on each Parc(p i , p i+1 , C, ) is allowed to move to p i+1 if it is not yet at this position. At the end of these motions, all the positions in F inal P ositions ∩ C M ap(Q,P) k+1

are occupied and remains occupied.

• At least one P -arc is free. In that case we have 1 ≤ γ < β. According to Algorithm 4, if a robot moves from a P -arc to another one then γ does not decrease because if robot r chooses to move from a P 1 -arc to a P 2 -arc, that implies that P 2 -arc is free. However, if P 1 becomes free when r reaches P 2 then the number of free P -arc remains unchanged.

We now assume by contradiction that γ never reaches the value β. So γ eventually remains inchanged. From this point on, no robot of any P -arc containing more than one robot will move towards a free P -arc. Following the algorithm, that implies that every P -arc P with more than one robot is followed in clockwise by a non free P -arc infinitely often (at least each time the last robot of P is awaked). Since the robots cannot move in counterclockwise, that also implies that every Parc with more than one robot is always followed by a non free P -arc P ′ . If this second P -arc P ′ also contains more than one robot then it is also followed by a non free P -arc. However, if P ′ contains only one robot then it is also followed by a non free P -arc P ′′ , since on the contrary, the robot of P ′ will eventually move to P ′′ and the last robot of P will eventually move to P ′ . A contradiction. So, step by step, it is clear that no P -arc can be free and γ = β which contradicts our assumption. So γ will eventually reach the value β.

When γ = β we retrieve the case where all the P -arcs are not free and the lemma holds.

2

Lemma 36 Let P be a target pattern and let Q be a configuration which is a maximal (0, P)-partial pattern. If all the extra robots are on SEC(M ap(Q, P)) then, all the positions in F inal P ositions ∩ C M ap(Q,P) k+1

are occupied in a finite number of cycles.

Proof.

If all the extra robots are on SEC(M ap(Q, P)) and there exists at least one position in F inal P ositions ∩ SEC(M ap(Q, P)) which is not occupied then the robots apply the routine Arrange(C M ap(Q,P) k+1

, F inal P ositions) (refer to Algorithm 4) for k = 1. Despite a more complicated code, the case k = 1 can be seen as the case k > 1 with an additionnal constraint on the angles and a particular statement for a deadlock configuration removal. We show (refer to last item of this proof) that the deadlock removal generates a behavior that can finally be generated by Algorithm 4 for a concentric enclosing circle which is not SEC(Q). So the aim of the proof is to show that Algorithm 4 has no deadlock. In the rest of this proof we say that a point p is a P -point if p ∈ F inal P ositions ∩ SEC(M ap(Q, P)). So assume by contradiction that there exists a deadlock and we consider the two following cases:

1. No Parc is free but there exists at least one P -point which is not occupied by a robot.

Again, we distinguish two cases:

(a) At least one P -point is occupied by a robot. Let p i be one these P -point such that its successor in clockwise p (i+1) is free. Clearly, angle(p i , c, p (i+1) , ) ≤ 180 (even if the first non critical s does not belong to F inal P ositions because, due to the fact s is not critical, from Lemma13 its absence cannot create an angle > 180). So the last robot of the p (i+1) Parc can move to p (i+1) . A contradiction.

(b) No P -point is occupied by a robot. Since there are at least three robots on SEQ, at least one of them has a predecessor with an angle less than 180 o . So it can move. A contradiction.

2. There exists at least one free Parc. Let i (0 ≤ i ≤ α -1) be an integer such that the ith Parc is free and its predecessor (the ((i -1)modα) th Parc) is not. Let us call them A and A ′ , respectively. We distinguish two cases:

(a) There exists (A ′ , A) such that A ′ contains at least two robots. We call r the last robot of A ′ and r ′ the predecessor of r on A ′ . In this case r can move to A (since angle(r ′ , c, r,

(b) Every couple (A ′ , A) is such that A ′ contains one robot only. In that case there exists at least a couple (A ′ , A) such that the predecessor A ′′ of A ′ contains at least one robot since there are at least as many robots as Ppoints on SEQ. Because the deadlock assumption, the robot r on A ′ cannot move to A so angle(r ′ , c, r, ) = 180 o where r ′ is the last robot on A ′′ . Again, we distinguish two cases:

i. r is not on p (i-1)modα . In this case r ′ also is not on p (i-2)modα since angle(p (i-2)modα , c, r, ) < 180. Since there is at least a third robot on SEQ, this robot r ′′ is such that angle(r ′′ , c, r ′ , ) < 180 so r ′ can move toward p (i-2)modα . A contradiction. ii. r is on p (i-1)modα . In this case r ′ is also on p (i-2)modα and angle(p (i-2)modα , c, p (i-1)modα , ) = 180. Since no robot can move, we can see that the configuration on SEQ is as follows: A is the first Parc of a chain starting from A in the conterclockwise such that any Parc of this chain but A contains a robot at its P -point, we call this chain P C. The last Parc of P C is followed by a free Parc (A if there exists no other free Parc). Since no robot can move we can deduce that between this free Parc and A (in the conterclockwise) all theParcs are free. So all the robots are on P C and there exists at least one Parc of P C which contains at least two robots. Let B be the first Parc of the chain (starting from A in conterclockwise) such that B contains at least two robots. Then the chain starting from A and ending to B is a deadlock chain. By definition, the robot on the P -point of B is a deadlock breaker and can move. A contradiction. Now we just focus on the behavior of the successive deadlock breakers. The aim of their behavior is to allow r to move toward the next P -point. It is easy to see that this part of the algorithm just reverses the order of the deadlock breakers and r, but once ony of these robots has started to move their behavior is the same as in the internal circle part (still with angle constraint).

2

Lemma 37 Let P be a target pattern and let Q be a configuration which is a maximal (k, P)-partial pattern. The robots form a (k + 1, P)-partial pattern or the target pattern is formed, in a finite number of cycles.

Proof. We have to consider three cases.

• k = |SC M ap(Q,P) |. In that case, C M ap(Q,P) k+1

does not exist and

∩ M ap(Q, P). That implies that it remains only one position p to occupy and p is inevitably at the center of SEC(Q) (otherwise C M ap(Q,P) k+1 would exist). According to Algorithm 2, leader r l moves toward c. From Corollary 33 and by fairness, we deduce that the target pattern is formed in a finite number of cycles.

• k = |SC M ap(Q,P) | -1. In that case, we distinguish two subcases:

1. The center c of SEC(Q) is in M ap(Q, P). In that subcase, all the positions in C which is not occupied by any robot. Of course, if we are lucky, a (k + 1, P)-partial pattern is formed during this step. Otherwise, the robots apply Algorithm 4 and, from Lemma 35 and 36, the (k + 1, P)-partial pattern is formed in a finite number of cycles. , Q, r). If we are lucky, a (k + 1, P)-partial pattern is formed during this step. Otherwise, the robots apply Algorithm 4 and, from Lemma 35 and 36, the (k + 1, P)-partial pattern is formed in a finite number of cycles.