Pierre Alart
email: pierre.alart@lmgc.univ-montp2.fr

David Dureisseix
email: david.dureisseix@lmgc.univ-montp2.fr

A scalable multiscale LATIN method adapted to nonsmooth discrete media

Keywords: keyword domain decomposition, scalability, multilevel, nonsmoothness, homogenization PACS 45.10.Db, 46.15.Cc, 46.70.Lk

The simulation of discrete systems often leads to large scale problems, for instance if they result of a discretization technique, or a modeling at a small scale.

A multiscale analysis may involve an homogenized macroscopic problem, as well as a coarse space mechanism to accelerate convergence of the numerical scheme. A multilevel domain decomposition technique is used herein as both a numerical strategy to simulate the behaviour of a non smooth discrete media, and to provide a macroscopic numerical behaviour of the same system.

Several generic formulations for such systems are discussed in this article. A multilevel domain decomposition is tested and several choices of the embedded coarse space are discussed, in particular with respect of the emergence of weak interfaces, characteristics of the discrete media substructuration. The application problem is the quasi-static simulation of a large scale tensegrity grid.

Introduction: discrete multiscale systems

The study of the discrete systems ensues from various motivations. Some systems are naturally discrete on the scale of the usual modeling; it is there so in mechanics of the spatial reticulated structures (trusses of bars); the discrete approach supplies then an effective model. When the number of elements increases, it can be useful to opt for a continuous model via an homogenization procedure [START_REF] Tollenaere | Continuous modeling of lattice structures by homogenization[END_REF]; we thus distinguish the discrete scale of the elements from the global scale of the structure. But this approach remains exceptional.

Discrete-to-continuum approaches

More usually the discrete models intervene to enrich locally initially continuous models; it is about a natural outcome of the micro-macro approach become classical in computational material science. It was initially a question of taking into account the material heterogeneousness by making interact various models remaining continuous at their scales. The development of the computational means led to a descent towards scales close to the atom where the models are discrete. The jump in scales is however considerable and the ab initio computations, very used by physicists and chemists [START_REF] Cimpoesu | The ab initio analytical approach of vibronic quatities : Application to inorganic stereochemistry[END_REF][START_REF] Cancès | Computational quantum chemistry: A primer[END_REF], provide only qualitative information on the macroscopic behaviour of the material. The numerical simulation is used here as a quasi-experimental tool to investigate the possible collective behaviour of complex systems comprising a very large number of elements, molecules, atoms, and various interactions between them. The molecular dynamics draws from this issue the main part of its success.

To go up from the atomic scale to the continuum one several stages and corresponding scales are often necessary. In these so-called mesoscopic scales, systems, sometimes still discrete as the dislocations, can interact to supply the plastic behaviour of the single-crystal [START_REF] Bulatov | Dislocation modelling at atomistic and mesoscopic scales[END_REF][START_REF] Lemarchand | Homogenization method for a discrete-continuum simulation of dislocation dynamics[END_REF]. Without resorting always to well identified scales, a discrete model can help to report exotic behaviour as metastable states in phase transformations [START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF][START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF]. Some authors tried to propose a discrete-to-continuum bridging [START_REF] Bendhia | Problème de mécanique multiéchelle: la méthode arlequin[END_REF][START_REF] Fish | Discrete-to-continuum bridging based on multigrid principles[END_REF], with the aim of carrying out a discrete structural zoom on zones of large deformations or defects, as fracture or buckling in the nanotubes [START_REF] Arroyo | A finite deformation membrane based on interatomic potentials for the transverse mechanics of nanotubes[END_REF][START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF].

The granular media lend themselves hardly to a continuous modeling, especially when they are dynamically loaded while remaining dense. Their behaviour is neither fluid nor solid, but takes features of these two models. The discrete approach turns out consequently particularly useful and, associated to sometimes large calculations, supplied fine results that are experimentally unobservable as the bimodal character of the stress network in granular packings [START_REF] Radjai | Bimodal character of stress transmission in granular packings[END_REF].

In some cases, very different scales are to be taken into account. The case of TexSol, mixture of sand and wires is particularly interesting. Very long textile wires are incorporated into the sand; the wire diameter is much smaller than the average of that of the sand grains and the wire length is almost infinite with respect to the grain size and induces long-distance interactions. A straightforward approach consists in modeling the wire as a continuous one-dimensional medium, i.e. a bending beam. But the coupling in the same model of a granular model for the sand and a continuum one for the wire imposes an unwieldy finite element discretisation of the wire. In [START_REF] Laniel | Numerical modelling of reinforced geomaterials by wires using the non smooth contact dynamics[END_REF], we opt directly for a discrete model of wire as "chain of pearls": every pearl is connected to the other ones by a bilateral link or a spring and is in frictional contact with the nearby grains.

Then a micro-macro approach cannot aspire still to a predictive status which requires a rigorously mathematical scale change. Such a topic is essential for granular media for which the definition of a representative elementary volume stays an open question. Especially since the choice of the equivalent effective continuum model is not evident. Is a first gradient theory relevant? Does it need to introduce a second gradient theory? Or a Cosserat like medium? These questions are valid for standard granular media, and even more for TexSol for which the interactions at a distance via the wire have to be considered [START_REF] Frémond | Non smooth thermomechanics[END_REF].

Diffuse non smoothness

If we define discrete systems as a set of nodes and of links between them, the non smoothness may only occur in the constitutive relations of the links. The molecular dynamics involves only smooth interaction laws between particles even if these interactions may be strongly non linear using Lennard-Jones, Landau-Gizburg or well potentials. A granular medium constitutes a typical example of a fully non smooth system for which all the links involve non smooth behaviour: unilateral contact and dry friction between the grains. Such a system is the main and ultimate motivation of the present numerical study, but it is quite difficult to define rigorously onto it a multiscale approach and to proceed to validation tests. Indeed a substructuring of a granular medium has to be frequently updated to take into account the changes of the connectivity of the particles. It is a technical difficulty but so an obstacle to the definition of a permanent subdomain having an homogenized behaviour. Then it is more convenient to consider a structural system that preserves the connectivity of the elements with a static or a dynamical loading. The non smoothness may be introduced by the unilateral behaviour of cables. When such a behaviour occurs through the whole structure, the non smoothness is called diffuse. The concept of tensegrity structure [START_REF] Fuller | The dymaxion world of Buckingham Fuller[END_REF][START_REF] Motro | [END_REF] is then recovered and the modeling of a large range of reticulated structures is investigated in the next section.

Domain decomposition of discrete systems

Discrete systems, especially for large scale problems, often require specific solvers, according to the formulation (static or dynamical...). Section 2 is devoted to the presentation of a wide range of formulations, from static to dynamical behaviours, and accounting for different degrees of non smoothness.

Section 3 recalls a domain decomposition approach that serves herein as a basis, as detailed in a previous paper [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF]. One particularity for discrete systems is the arising of weak interfaces between substructures. These interfaces are used to select the coarse space of the multilevel domain decomposition method in Section 4. The bridge with the FETI-DP method is also discussed.

Finally, Section 6 proposed some numerical tests on a 3D grid modeled as a linear truss or as a nonsmooth tensegrity grid. Numerical scalability is exemplified when increasing the size of the problem together with the number of subdomains, while maintaining the nonsmoothness ratio of the problem constant.

2 From a smooth system to a fully non smooth one via tensegrity

Even if the quasi-static situation is only investigated in the next sections using a domain decomposition strategy, both quasi-static and dynamical modelings of a structure are discussed in this section. Indeed the features of the final equations are quite similar in the two cases and a comparable strategy may be performed. To illustrate the concepts introduced above we consider then a structure composed with bars and more or less cables; such a discrete structure may be described as a set of nodes and links between them, the non smoothness only occuring in the constitutive relations of the links. In the set of nodes Ω, we distinguish the subset Γ u of the nodes where the displacement is prescribed to clamp the structure to the support. The basic equations are given in the two first subsections for the quasi-static and dynamical cases. Thus some situations with more or less cables are investigated underlining the common features of the formulations, and some differences.

Quasi-static process

Three configurations are to be considered: the current one Ω 1 for which the tensions and displacements are unknown, the prestressed configuration Ω 0 before applying additional external loading and the relaxed configuration Ω -1 for which the selfstresses are virtually vanished. The three configurations are assumed to be close enough to preserve the principle of small perturbations and the prestresses are assumed to be given. In Table 1, the main notations are introduced.

b, c

Bar and cable subscripts τ 0 , τ 0 b , τ 0 c such that B t τ 0 = 0 Self balanced prestress e 0 , e 0 b , e 0 Unknown and prescribed nodal displacements B t , B Link-to-node and node-to-link mappings

K b = B t b k b B b , K c = B t c k c B c
Bar network and cable network stiffnesses

K = K b + K c
Global stiffness Table 1: General notations for a system of bars and cables and a quasi-static process.

The balance equation links the external nodal forces F d to the internal tensions via the link-to-node mapping B t .

-F + F d = 0 with F := B t τ. (1)
The strain admissibility equation connects the nodal displacements to length variations of the links via the node-to-link mapping B as defined in Table 1.

The behaviour relation for a bar indiced by α involves a local stiffness between the tension and the length variation taking into account the prestress,

τ α = k α (e α + e 0 α ). (2)
An inextensible cable may be modeled with complementary conditions expressed as follows, 0 ≤ -e α ⊥ τ α ≥ 0.

The previous notation summarizes three relations: the two inequalities and the complementary condition using the orthogonality symbol. The modeling of extensible cables is a combination of a bar and an inextensible cable; the behaviour law takes the form of a piecewise linear function.

τ α = k α (e α + e 0 α ) if e α + e 0 α > 0 0 if e α + e 0 α ≤ 0 (4)
But we can easily prove that this relation is equivalent to a complementary condition between the tension and a corrected length variation λ α , defined in Table 1,

0 ≤ λ α ⊥ τ α ≥ 0. ( 5
)

Dynamical process

When the dynamics is considered, and when the masses are concentrated in the nodes, some additional variables have to be introduced as summarized in Table 2.

η = -ė = -B U = -BV
Relative velocities (strain rate admissibility)

π = t + t -τ dν Average impulsions λ + c = η + c + k -1 c h -2 π c -h -1 (e - c + e 0 c ) Corrected relative velocities in cables M Mass matrix W = BM -1 B t = W bb W bc W cb W cc Delassus operator W = W + h -2 diag(k -1 b , k -1 c ) Corrected Delassus operator = Wbb W bc W cb Wcc
Table 2: Additional notations for a system of bars and cables and a dynamical process.

For smooth motions the dynamical equation involves the time-derivative of the velocities. Since shocks are expected, it is more convenient to write this equation as a measure differential equation ( [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], [START_REF] Jean | The non-smooth contact dynamics method[END_REF]),

M dV + B t τ dν = F d dt. ( 6
)
where dt is a Lebesgue measure, dV is a differential measure representing the acceleration, dν a non-negative real measure relative to which dV happens to possess a density function, and τ is a representative of local density of tension forces. A dual (or reduced to links) formulation of the dynamics may be preferred using the Delassus operator W ,

dη -W τ dν = -F d dt = -BM -1 F d dt. (7)
When a time discretisation is performed an elementary subinterval [t -, t + ] of length h is considered. The variables evaluated at t -, respectively t + , have the -, respectively + , superscript. Since discontinuous velocities are expected, high order integration schemes are not necessary and even troublesome; first-order schemes are enough when many shocks may occur simultaneously. We consider here the implicit Euler scheme underlining the impulsion π over the time step as the product of the time step h by an average tension τ + considered at the end of step (note that the + superscript may be then omitted),

η + -η --W π = -h F d with π = hτ + . ( 8
)
The dynamical version of the behaviour relation ( 2) for a bar involves the relative velocity η + α and the impulsion π α ,

π α = h 2 k α [-η + α + 1 h (e - α + e 0 α )]. (9)
According to the approach of [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], a dynamical discrete version of (3) for an inextensible cable is derived involving complementary conditions between relative velocity and impulsion,

if -e - α > 0 then τ + α = 0, if -e - α ≤ 0 then 0 ≤ η + α ⊥ π α ≥ 0. ( 10
)
An integration lemma given by [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] proves that the iterates verifying (10) tend to verify (3) when the time step h tends to zero. For an extensible cable the combination of the behaviour law (4) and the integration scheme provides the following relation,

τ + α = k α (e - α -hη + α + e 0 α ) if e - α -hη + α + e 0 α > 0 0 if e - α -hη + α + e 0 α ≤ 0 (11)
According to previous developments we can still recover a complementary formulation linking the impulsion to a corrected relative velocity λ + c defined in Table 2,

0 ≤ λ + α ⊥ π α ≥ 0. ( 12
)
With these ingredients we can postulate some problems with different mathematical features.

Truss

If the system is only composed by bars, the equilibrium of the so obtained truss is classically characterized by a linear system with the nodal displacements as unknowns,

KU = F d -B t ke 0 . (13)
The (smooth) dynamical behaviour is governed by a system of second-order differential equations,

M Ü + KU = F d -B t ke 0 . ( 14
)
If some inextensible cables are added -for instance to hang up the structure to the support because it is too flexible and too heavy -the equilibrium depends on the tension in these cables according to a few complementary conditions,

K b U + B t c τ c = F d -B t b k b e 0 b 0 ≤ -B c U ⊥ τ c ≥ 0. ( 15
)
It is then a structural type problem as defined in [START_REF] Alart | Parallel computational strategies for multi-contact problems: Applications to cellular and granular media[END_REF]; the global stiffness matrix of the truss of bars is not singular and the system (15) is equivalent to minimize a lower bounded quadratic bulk energy under convex constraints i.e. a well-posed problem. Such a problem does not involve a diffuse non smoothness.

Net and granulate

If a lot of bars are replaced by inextensible cables (to make the structure lighter for instance) the matrix K b may be singular and the problem is no more wellposed. The dynamical problem may be solved more easily and takes the form of a linear complementary problem (LCP),

Ŵc π c -η + c = F d c 0 ≤ η + c ⊥ π c ≥ 0. ( 16
) with Ŵc = W cc -W cb W -1 bb W bc F d c = -η - c + h F d c + W cb W -1 bb [η - b -h F d b -1 h (e - b + e 0 b )] (17)
We can imagine to replace all the bars by cables -think of a catenary arch model of Gaudi (see Figure 1) or of a fishing net. The dynamical discrete behaviour then derives from the previous system with Ŵc equal to W and this system is similar to the one issued from the modeling of a granular system with frictionless contact between grains i.e. a granular type problem as defined in [START_REF] Alart | Parallel computational strategies for multi-contact problems: Applications to cellular and granular media[END_REF]. Such a problem is clearly fully non smooth.

Tensegrity

The stiffness of the cables is often weaker than the one of the bars. It is then convenient to consider extensible cables instead of inextensible ones. Such a structure is a selfstressed tensegrity system if the set of compressed components is discontinuous and the set of tensioned components is continuous, [START_REF] Motro | [END_REF]. Such a problem is then characterized by a diffuse non smoothness. The matrix K b is then singular with a large kernel composed of the rigid modes of all the bars. Since the global stiffness of the associated truss is invertible, the equilibrium of the system may be characterized by a LCP involving the tensions in the cables and the corrected length variations defined in Table 1,

A c λ c -τ c = -k c B c K -1 F d -τ 0 c 0 ≤ λ c ⊥ τ c ≥ 0. ( 18
)
with

A c = k c -k c B c K -1 B t c k c . (19)
The one step dynamical discrete problem is still a LCP when considering the corrected relative velocities λ + c defined in Table 2,

Wc π c -λ + c = F d c 0 ≤ λ + c ⊥ π c ≥ 0. with Wc = Wcc -W cb W -1 bb W bc F d c = F d c + 1 h (e - c + e 0 c ) (20)

A multiscale domain decomposition approach

The extension to 3D discrete system of the micro-macro strategy in [START_REF] Ladevèze | A micro / macro approach for parallel computing of heterogeneous structures[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF] (initially developed for continuum media, and tested for 2D problems), can be found in [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF]. As it serves as a basis for the present experimentations, and in order for this article to be self-contained, it is briefly recalled in the following.

Substructuring

The first step consists of a decomposition of the structure into substructures and interfaces. There are two common ways to partition a discrete domain: (i) distribute the elements among substructures, or (ii) distribute the nodes among substructures. Both strategies were used in the litterature. For a topological point of view, they correspond to the splitting of two different graphs: the connectivity graph or its dual graph. The node distribution if often considered in algebraic partitioning where the left hand side matrix is splitted according to its columns (degrees of dreedom).

In our case, since the non smoothness may occur within the element behaviours, we choose to distribute the elements among the substructures (see Figure 2). Indeed, with such a choice, some of the 'boundary' nodes are duplicated in the data structures arising for the subdomains. The interface between two subdomains is defined to be the set of these nodes, joining a substructure to the other. The nonsmoothness is therefore localized within the substructures only. This modeling choice is identical to [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF] and somehow the dual of the one proposed in [START_REF] Nouy | A multiscale computational approach for contact problems[END_REF] where the non linearities (contact in cracks) are isolated in the interfaces.

Such a decomposition allows to reformulate the reference problem. Basically, only the admissibility equations have to be reformulated. A substructure E is submitted to the action of its neighboring interfaces Γ EE ′ : forces F EE ′ and displacements U EE ′ . Extended to all the interfaces connected to the substructure E, the assembling of the previous fields are denoted with F EΓ and U EΓ . The balance equation for this substructure is

-F E + F d E + C t EΓ F EΓ = 0 with F E = B t E τ E ( 21
)
where C EΓ is a boolean mapping matrix that selects the trace on the local interfaces. The strain admissibility is:

e E = B E U E , U EΓ = C EΓ U E and U E | Γu = U d E ( 22
)
Figure 2: Element oriented partitioning (left) and perfect discrete interface between the substructures (right).

Moreover, the behaviour of the interfaces has to be settled. Due to the previous splitting choice, we now deal with perfect discrete interfaces: an interface Γ EE ′ transfers the forces F EE ′ and the displacements U EE ′ on each of its sides, i.e.

F EE ′ + F E ′ E = 0 and U EE ′ = U E ′ E ( 23
)

Micro-macro description

Following the approach proposed in [START_REF] Ladevèze | A micro / macro approach for parallel computing of heterogeneous structures[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF], we define the macroscopic quantities as generalized averages of mechanical fields on each interface Γ EE ′ : the generalized macro forces are f EE ′ = R t EE ′ F EE ′ , and the generalized macro displacements are u EE ′ = R t EE ′ U EE ′ , where R t EE ′ is a projector onto the macro space. An orthogonal projector which is identical for displacements and forces can be built, see [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF], such that

U EE ′ = R EE ′ u EE ′ + U m EE ′ and F EE ′ = R EE ′ f EE ′ + F m EE ′
, the superscript m denoting the micro additional quantities, which are orthogonal to the macro ones. Similar notations are used with a subscript EΓ that denotes the assembly of all the quantities on the neighboring interfaces of the subdomain E. The macro displacements are usually selected as mean translations, rotations, mean stretchings and shearing of the interface displacement, see Figure 3; the macro forces are the corresponding dual quantities. Therefore, there are 9 parameters on each interface to describe each macro field; they are stored in u EE ′ and f EE ′ . The basis vectors of these particular fields are stored in R EE ′ . Figure 3: Some macro displacement fields in R EE ′ on a 3D discrete interface: a mean translation, a mean rotation and a mean stretching.

The LATIN method as a solver

The LArge Time INcrement method (LATIN, [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]) is a template that can be used to design solving algorithms; for steady state problems, it is closely related to an augmented Lagrangian approach [START_REF] Alart | Using nonsmooth analysis for numerical simulation of contact mechanics[END_REF]. For the problem we intend to solve here, the algorithm exhibits two steps at each iteration:

• once an approximation of the solution (e E , τ E , U EE ′ , F EE ′ ) is known, the local stage looks for (ê E , τE , ÛEE ′ , FEE ′ ) satisfying a search direction for the substructure fields (τ E -τ E )+l(ê E -e E ) = 0

the constitutive relations of the links (2), (5);

a search direction for the interface fields

( FEE ′ -F EE ′ ) -d( ÛEE ′ - U EE ′ ) = 0
the perfect behaviour of the interfaces [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF].

d and l are two positive scalar parameters. This stage leads to local linear or non-linear problems, whose solutions are explicit.

• with the previous approximation, the linear stage looks for a new approximation (e E , τ E , U EE ′ , F EE ′ ) satisfying the admissibility constraints: balance (21) and compatibility [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF];

a search direction for the substructure fields (τ E -τ E )-l(e E -ê E ) = 0 -additional admissibility constraints: continuity of macro displacements throughout interfaces (this will lead to the coarse problem); -a weakened search direction for the interface fields to be consistent with the additional admissibility constraints.

The linear stage requires to be further detailed. In [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF], a force-oriented approach was used, though both a displacement and force oriented are available for continuum media [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF]. Here, the displacement-oriented one happens to be more suited to the treatment of weak interfaces in the next Section. Therefore, the presentation will summarize the micro-macro approach within this framework.

The macro continuity on interface Γ EE ′ is:

u EE ′ = u E ′ E .
Therefore there is only one macro displacement for each interface. For the whole set of interfaces, they can be stored in a macro vector denoted with u Γ . For the interfaces connected to the subdomain E, one gets

u EΓ = c E u Γ ( 24
)
where c E is a boolean matrix selecting the entries in u Γ (f Γ = E c t E f EΓ is the macro unbalance through the interfaces).

The search direction on interfaces is weakened for the macro part only. The micro part remains:

(F m EE ′ -F m EE ′ ) + d(U m EE ′ -Û m EE ′ ) = 0.
But for the macro part, the corresponding weakened expression, for an optimal choice of the macro parameter [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF][START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF] is merely

f Γ = 0 (25)
(i.e. no macro unbalance). We will still use this particular search direction in all of the following. Using the search directions and the compatibility,

τ E = (τ E -lê E ) + lB E U E (26)
F m EΓ = ( F m EΓ + d Û m EΓ ) -d(C EΓ U E -R EΓ u EΓ ) (27)
the balance equation now reads:

K E U E -C t EΓ dR EΓ (u EΓ + d -1 f EΓ ) = F d E ( 28
)
where K E = B t E lB E + C t EΓ dC EΓ is a regular matrix thanks to the regularizing effect of the search direction stiffness d on the interfaces, and

F d E = F d E -B t E (τ E - lê E ) + C t EΓ ( F m EΓ + d Û m EΓ
) is known at this stage. The relationship between macro displacement and internal displacement is:

u EΓ = R t EΓ C EΓ U E .
This allows to condense the whole problem onto macro dofs:

H E u EΓ = f EΓ + B -1 E f d E (29)
with

B E = R t EΓ dC EΓ (K E ) -1 C t EΓ dR EΓ which can be proved to be invertible, f d E = R t EΓ dC EΓ (K E ) -1 F d E is given, and H E = B -1 E (d 1 -B E )
is the homogenized stiffness of the subdomain. One has to note that the left hand side of ( 29) is not always regular. Indeed, if there are floating subdomains (i.e. subdomains for which the prescribed dofs on Γ u do not clamp every rigid body motion), the trace of their rigid body motions on the interfaces are in the coarse space, and since no energy is associated, they belong to the kernel of H E . This homogenized behaviour is involved in the domain decomposition coarse problem; indeed, to close the problem, one has to add the additional admissibility [START_REF] Nouy | A multiscale computational approach for contact problems[END_REF], and the corresponding search direction [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]. This leads to:

( E c t E H E c E )u Γ = - E c t E B -1 E f d E (30)
This is the coarse problem that allows to get the global macro displacements u Γ at each iteration. If the reference problem is well posed, this coarse problem is regular as well.

The relocalization into each subdomain independently is ensured by ( 28) and ( 29).

Weak interfaces and coarse space selection

For discrete systems, the substructuring may lead to 'weak' interfaces, i.e. interfaces containing a too small number of nodes to define a macro part as described above. Such a situation raises the question of the meaning of the coarse space (involving the macro quantities) with respect to the homogenized behaviour.

Figure 4 illustrates the point for a 2D discrete problem: in all of the depicted cases, there are 4 'strong' interfaces (namely 12, 23, 34 and 41); in the case on the left, there are 2 'weak' interfaces (13 an 24) since the corresponding couples of subdomains share 1 node each. In the other cases, only 1 weak interface ( 13) is detected. The continuity relationships enforced by the weak interfaces in the left case, are redundant with the perfect behaviour of the strong interfaces. Therefore, these weak interfaces can be discarded without changing the reference problem. This is what has been done up to now [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF], and this happens for a regular decomposition of a tensegrity grid designed with module duplication.

On the contrary for the other cases of Figure 4, the weak interfaces cannot be discarded. This is what happens for a random tensegrity network, or a granular media for which the links between the grains are the potential contact interactions. There are mainly two alternatives: (i) no macro part is selected on the weak interfaces, and the macro homogenized interpretation is not modified; (ii) an additional macro field has to be designed. This will be the topic of the next sections.

The 3D case is more complicated since a weak interface can be identified as 'corner' or 'wedge' between subdomains, while a strong interface is more likely to be located on a 'face'. 'Wedge-type' interfaces may have a sufficient number of nodes to define averaged fields (macro fields of Figure 3), but some of their components may be ill-defined due to the spatial distribution of the nodes. In such cases, these interfaces will still be selected as weak ones.

No macro field on weak interfaces

No change has to be done on the previous algorithm if no macro field is defined on weak interfaces where

U EE ′ = U m EE ′ , F EE ′ = F m EE ′ .
If no macro field is defined on weak interface, the problem becomes illconditioned since the rigid body motions of floating subdomains are determined by these weak interfaces only (roughly, a floating subdomain is kept with springs on its corners only). As a consequence, the convergence rate decreases.

Full macro field on weak interfaces

The second possibility considers the whole behaviour of the weak interfaces as a part of the coarse space. The difference with the previous one lies in the fact that no more micro quantities on weak interfaces remain.

If we define the macro field on the weak interfaces as the whole field, we change the homogenization interpretation of the macro problem, but we allow again the rigid body motions of floating subdomains as part of the macro field. This approach will be denoted with the 'macro' version.

To be compatible with the previous approach without weak interface, we define the macro part on weak interface with

U EE ′ = R EE ′ u EE ′ , F EE ′ = R EE ′ f EE ′ , with R EE ′ = 1.
Figure 5, center, illustrates this 'macro' version of the admissibility constraints in a 2D continuous problem.

The consequence is the increase of the macro problem size. A modification of this situation is possible: it is inspired by the corner modes of the FETI-DP method, that will be briefly recalled in the next Section.

The continuum media case: example of the FETI-DP method

For a linear continuum media, discretized with finite elements, the internal forces are obviously F E = K E U E where K E is here the finite element stiffness of the subdomain E. The genuine FETI-DP method [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part I: a faster alternative to the two-level FETI method[END_REF] proposes to enforce the displacement continuity of the so-called 'corner' nodes (usually, the nodes at the intersection of more than two subdomains in 2D, but other choices may be possible [START_REF] Lesoinne | A FETI-DP corner selection algorithm for three-dimensional problems[END_REF]). These corner modes are denoted with a subscript 'C', while the other ones are denoted with a subscript 'r':

U E = U Er U EC
Therefore, one prescribes U EC = p E u C where p E is a boolean matrix and u C stores all the 'corner' dofs. The interfaces between subdomains are concerned with 'r' dofs only, and the equilibrium of all subdomains may reads:

∀U ⋆ Er , u ⋆ C E U ⋆ Er p E u ⋆ C t (-K E U E + F d E + C EΓ F Er ) = 0 ( 31
)
F Er is the action of the neighboring interfaces on the subdomain r dofs. Though the solving algorithm used in FETI-DP method is a conjugate gradient acting on interface forces, there is an equivalent of the linear stage for which F Er is given by a search direction, that reads from [START_REF] Quirant | Designing tensegrity systems: the case of a double layer grid[END_REF], once the splitting in C and r dofs is used:

K E rr U Er + K E rC p E u C = F d Er + C EΓ F Er (32) ( E p t E K E CC p E )u C + E p t E K E Cr U Er = E p t E F d EC ( 33
)
Condensing the r dofs on the corner dofs gives the coarse problem:

E p t E (K E CC -K E Cr (K E rr ) -1 K E rC )p E u C = E p t E [F d EC -K E Cr (K E rr ) -1 (F d Er + C EΓ F Er )] (34)
Such a multilevel domain decomposition method proved to be scalable for 2D and plate problems. The numerical scalability of a domain decomposition method is its ability to solve with a quasi-constant number of iterations, a set of problems of increasing size and an increasing number of subdomain (such that the subdomain sizes remain constant). This is an important property for the approach to be able to manage efficiently large scale problems.

For 3D problems, the scalability is lost. It is back recovered with an augmented algorithm featuring the addition of an averaged continuity enforcement throughout the interfaces [START_REF] Farhat | A scalable dual-primal domain decomposition method[END_REF]. Note that when the corner nodes are dropped, the scalability is somehow lost again [START_REF] Klawonn | Some computational results for dual-primal FETI methods for three dimensional elliptic problems[END_REF].

Corner modes on weak interfaces

The FETI-DP method, based on corner modes, may be enriched with interface average continuity conditions, whereas we proposed herein to enrich the interface-based macro behaviour with weak interface modes.

We take the opportunity of the weak interface appearance to build an augmentation of the previous approach. Since there will be a global continuity of the weak interfaces, this version will be denoted with the 'corner' version. Figure 5, right, illustrates this version of the admissibility constraints in 2D.

In the FETI-DP method the corner displacement continuity is expressed within the equilibrium. Nevertheless, the framework we used here is the generalization of the LATIN micro/macro approach: the weakening of the search direction on the weak interfaces. This approach is detailed in the next Section.

5 A multiscale approach embedding a corner mode constraint

Kinematics description

The displacement U E of the substructure E is split into C dofs (that stands for 'corner' dofs, here: the dofs related to the nodes on the weak interfaces of substructure E) and the remaining ones, called r dofs.

U E = U Er U EC (35)
The trace of the displacement on the strong interfaces of the subtructure E is still U EΓ = C EΓ U E . Note that U E is used in this expression, and not U Er , which is a slight difference when compared to the continuum media case of the FETI-DP approach, but it is mandatory for a discrete media (especially when there is no redundancy between strong and weak interfaces).

Sthenics description

The equilibium of the substructure E is now,

-F E + F d E + C t EΓ F EΓ + 0 F EC = 0 (36)
F EC can be viewed as a Lagrange multiplier that will be used to enforce the corner displacement continuity. Note that both for the displacements and the forces, the micro/macro description on strong interfaces is unmodified.

Modifying the admissibility conditions within the linear stage

The admissibility conditions are augmented with (i) the macro displacement continuity

u EΓ = c E u Γ (37)
on the strong interfaces, and (ii) the corner displacement continuity

U EC = p E u C (38)
The optimal search direction corresponding to the first augmentation, (25), is still

E c t E f EΓ = 0 (39)
The search direction corresponding to the second augmentation is also weakened:

E U * T EC d C (U EC -ÛEC ) + d -1 C (F EC -FEC ) = 0 for all U * EC satisfying U * EC = p E u * C .
This leads to:

a(u c -ûC ) + E p t E (F EC -FEC ) = 0 (40)
where a = E p t E d C p E . One property of the solution arising from local stage is that E p t E FEC = 0. Moreover, if one chooses d C = 0, then a = 0 and the search direction is merely:

E p t E F EC = 0 (41)
The admissibility finally reads:

K E U E -Q E u * EΓ - 0 F EC = F d E (42)
where

• K E is the same matrix as in [START_REF] Lesoinne | A FETI-DP corner selection algorithm for three-dimensional problems[END_REF],

• Q E = C t EΓ dR EΓ , • u * EΓ = u EΓ + d -1 f EΓ , • F d E = F d E -B t E (τ E -lê E ) + C t EΓ ( F m E + d Û m E
) is given at this stage. In addition, we get the macro gluing within a subdomain u EΓ = R t EΓ U EΓ , that can be expressed as:

-Q E U E + du * EΓ = f EΓ (43)
With the splitting into r and C dofs, the problem (42)-(43) reads:

  K E rr -Q Er K E rC -Q t Er d I d -Q t EC K E Cr -Q EC K E CC     U Er u * EΓ U EC   =   F d Er 0 F d EC   +   0 f EΓ F EC   ( 44
)
where I d is the identity matrix.

BE fE = KE⋆ cE u Γ u C -F ⋆ E (48) with BE = I d -KE⋆ d -1 I d 0 0 0 (49)
BE is not symmetric, but it is invertible: indeed if one prescribes the macro displacements u EΓ and U EC , there are unique corresponding forces f EΓ and F EC given by (48) (the reverse is not true if the clamped corner dofs do not prevent all rigid body motions on the floating subdomain). These macro forces are therefore given by:

fE = HE cE u Γ u C -B-1 E F ⋆ E (50)
with HE = B-1

E
KE⋆ which is the homogenized compliance of the substructure E, relating macro forces to macro displacements. Though BE is not symmetric, HE is symmetric: with algebraic manipulations, it reads:

HE = d I d N t E d -1 B -1 E d I d N E + -d I d 0 0 K E, * CC (51)
with

N E = -Q t EC + Q t Er (K E rr ) -1 K E rC and B E = (Q t Er (K E rr ) -1 Q Er )d -1
. HE is the homogenized macro stiffness of the substructure E. It is not obviously invertible. Indeed it is not if the clamped corner dofs do not prevent all rigid body motions on the floating subdomain.

The assembly of all subdomain macro contributions is completed with the optimal macro search directions (39) and (41) that reads:

E ct E fE = 0 and leads to:

E ct E HE cE u Γ u C = E ct E B-1 E FE (52)
This is the coarse homogenized problem, for which the macro unknowns are the generalized displacements on the strong interfaces, and the corner displacements on the weak interfaces.

Numerical results

In order to test the proposed approaches, a first problem deals with a 3D truss build upon the tensegrity grid network of Figure 6. This is a linear problem, since there is no more unilateral behaviours on the elements of the truss (all are bars).

As boundary conditions, the lower nodes on two opposite bottom edges of the domain are clamped, and a uniform vertical force field is prescribed on every node. The numerical scalability test consists of increasing the number of subdomains, while keeping the problem per subdomain at constant size (16 modules per subdomain). Table 3 recalls the characteristics of the problem set for a grid that will be used with two modelings: one for a regular truss, the second for a non smooth tensegrity grid.

In all of the following tests, the search direction parameter is kept fixed, in order to compare the convergence rates. The proposed approach, at least without a multiscale feature, has been identified as an augmented lagrangian type approach [START_REF] Alart | Using nonsmooth analysis for numerical simulation of contact mechanics[END_REF]. In such cases, there is an optimal value for the search direction parameter d. As a first guess, it is chosen as the stiffness of springs at interface nodes that gives the same global stiffness as a whole subdomain (without slack cables). This evaluation is performed once on a typical subdomain. With the subdomain size and the chosen tensegrity module, whose characteristics are given in [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF], one gets d = 3.9 10 6 N/m. Table 4 reports the convergence for the two previous approaches named as 'macro' and 'corners'. Numerical scalability is exemplified for both of them. Moreover, their convergence rates are identical, though the number of additional degrees of freedom for the weak interfaces increase slower for the second approach; therefore, the coarse problem is smaller for the 'corner' approach when the size of the problem increases.

Note that for the tested grid, built upon an elementary module depicted in Figure 7 (left), there are 3 internal infinitesimal mechanisms per module; their assembly cancels all but one global grid mechanism (i.e. without strain energy). This last mechanism is finally clamped with the boundary conditions [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF]. Since such a mechanism can be interpreted as a global mode that participates to the corse space, a second grid is built upon a second module that does no possess such mechanisms, Figure 7 (right). Convergence rates are very similar for this second design, as reported in Table 5. Clearly, the optimal parameter d for the second module should be stiffer. Nevertheless, the same value as for the first test case is still used. Finally, the grids built upon the first are used with a tensegrity modeling: the initial module is composed of 12 cables and 4 bars [START_REF] Quirant | Designing tensegrity systems: the case of a double layer grid[END_REF]. When increasing the size of the grid, the external loading is decreased accordingly in order to keep the non smoothness 'level': the percentage of slack cables is roughly identical (the target ratio for slack cables was 9 %). Table 6 reports the number of iterations to reach convergence, here defined with an error of 0.1 % (and for an error of 0.01 % for comparison purpose). The error is estimated with respect to a refence solution that is computed once for all with the same algorithm but a stopping error indicator at least two order of magnitude below the current required error to detect convergence.

A reference solution could be computed with other solvers should use the formulation ( 18), [START_REF] Jean | The non-smooth contact dynamics method[END_REF]. It becomes rapidly unusable due to the term K -1 B t c that is very costly to compute as the size of the problem increases. A mono-domain strategy similar to the one used here could also be used to solve the reference problem. This have been done in [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF], where the reader can find the comparison between both solvers.

For the two strategies that are proposed herein, the convergence rate is still similar. For the nonsmooth case, the numerical scalability is less obvious, though the number of iterations increase slowly with respect to the number of subdomains.

To assess the influence of the search direction parameter d, the last case with 144 subdomain is used. Several simulations are performed with different values of the parameter, for the 'corner' strategy, and the iteration numbers at convergence are gathered in Table 7. There is clearly an optimal value, which is located below the approximate guess selected herein (let us recall that the guess is computed from a subdomain without any slack cable, which is obviously stiffer than with 21 % of slack cables).

Figure 8 illustrates the results obtained with an extreme loading (for which 21 % of the cables slacks), with the deformation of the grid, and with 3D plots reporting as the vertical direction the level of internal forces in each element, i.e. compression level in bars and tension level in cables, Figure 9. The structure still possess a stiffness reserve.

Conclusion and prospects

The proposed approach is concerned with large scale non smooth problems solved with a multiscale substructuring approach. The application is the quasistatic behaviour of highly loaded tensegrity structures, when a significant number of cables slacken.

Weak interfaces (containing a too small number of nodes to define a macro part), are characteristics of a discrete system substructuring. The proposal in this article is to use them to increase the admissibility constraints in the previously designed approach. This is connected to the use of corner modes in the multiscale domain decomposition FETI-DP. This last method, based on corner nodes, may be enriched with interface average continuity conditions, whereas we enrich the interface-based macro behaviour with weak interface modes.

These weak interfaces are automatically produced by the domain decomposition. Additionally, other similar interfaces could be added in an artificial way, but this is not used herein.

As a second-hand product, the algorithm also provides a numerically homogenized behaviour of the substructures. The context of algorithm augmentation with weak interface constraints raises the question of the meaning of the coarse space involving the macro quantities with respect to this homogenized behaviour.

A future development concerns the the dynamical behaviour of such tensegrity structures, as well as large scale granular media simulations where non smoothness occurs with contact and friction between grains.

c

  Related prestrain τ , τ b , τ c Internal tensions, in bars, in cables e = BU , e b = B b U , e c = B c U Length variations (strain admissibility) k b , k c Local stiffnesses (in tension for cables) λ c = -(e c + e 0 c ) + k -1 c τ c Corrected length variations in cables F , F d Internal and external nodal forces U , U d

Figure 1 :

 1 Figure 1: Catenary arch model (Gaudi museum) and the Needle Tower (tensegrity).

Figure 4 :

 4 Figure 4: Several cases for interface detection once a decomposition has been performed on a discrete media.

Figure 5 :

 5 Figure 5: Illustration of different macro gluings on interfaces: no gluing on weak interfaces (left), macro modes (center), corner modes (right).

Figure 6 :

 6 Figure 6: Substructuring of a tensegrity grid with 16 subdomains.

Figure 7 :

 7 Figure 7: Used modules (initial module on the left and modified module on the right).

Figure 8 :

 8 Figure 8: 144 subdomain tensegrity grid, vertical displacement (amplified twice) / m.

Figure 9 :

 9 Figure 9: 144 subdomain tensegrity grid, compressions in bars (left) / N, tensions in cables (right) / N.

  The local coarse quantities are therefore the macro displacements on strong interfaces u EΓ as well as the local corner dofs U EC . If we use the concise notations:

	then the problem (44) reads
			K E rr -Qt E	-QE KE		U Er ūE	=	F d Er FE	+	0 fE	(45)
	with QE = -Q Er K E rC and KE = Solving this system for all macro loadings, i.e. condensing dofs U Er onto d I d -Q t EC . -Q EC K E CC
	macro dofs ūE reads:	KE⋆ ūE = F ⋆ E + fE	(46)
	with KE⋆ = KE -Qt E (K E rr ) -1 submatrix K E rr is invertible), and F ⋆ QE (since K E is symmetric, positive definite, its E = FE -Qt E (K E rr ) -1 F d Er .
	The additional admissibility conditions (37) and (38) reads:
			ūE = cE	u Γ u C	+	d -1 I d 0 0 0 fE	(47)
	with cE =	c E 0 0 p E	. Therefore, one gets:
		ūE =	u * EΓ U EC		fE =	f EΓ F EC	FE =	0 F d EC

Table 3 :

 3 3D regular grid -characteristics of a serie of problems with increasing size.

			16	25	36	64	100	144
	nb of modules	144	256	400	576	1024	1600	2304
	total nb of dofs	1443 2499 3843 5475 9603 14883 21315
	total nb of elements	2304 4096 6400 9216 16384 25600 36864
	nb of strong macro dofs	108	216	360	540	1008	1620	2376
	nb of additional weak macro	24	54	96	150	294	486	726
	dofs for 'macro' version
	nb of additional weak macro 120	27	48	75	147	203	363
	dofs for 'corners' version
	nb of subdomains 9 16 25 36 64 100 144
	'macro'	27 28 29 30 30 31	31
	'corners'	27 28 29 29 31 31	31

Table 4 :

 4 3D regular truss -number of iterations to get an error less than 0.01 %.

Table 6 :

 6 3D regular tensegrity grid -number of iterations to get an error less than 0.1 % (within parenthesis: less than 0.01 %).

	. A direct solve

Table 7 :

 7 3D regular tensegrity grid with 144 subdomains -number of iterations to get an error less than 0.1 % (within parenthesis: less than 0.01 %) versus search direction parameter.

	d / (10 6 N/m)	7.84	3.92	1.96	0.98	0.49	0.245
	nb of iterations 68 (108) 39 (60) 23 (36) 17 (37) 20 (63) 29 (103)