D Dureisseix

C Farhat

A numerically scalable domain decomposition method for the solution of frictionless contact problems

Keywords: decomposition, FETI, scalability

We present a domain decomposition method with Lagrange multipliers for solving iteratively frictionless contact problems. This method, which is based on the FETI method and therefore is named here the FETI-C method, incorporates a coarse contact system that guides the iterative prediction of the active zone of contact. We demonstrate numerically that this method is numerically scalable with respect to both the problem size and the number of subdomains.

INTRODUCTION

Contact problems are frequent in structural analysis, particularly when investigating the assembly of substructures. They are characterized by constraints such as non-penetration conditions, and an active area of contact -that is, an area where contact effectively occurs -that is unknown a priori. For these reasons, these problems lead to stiff non-linear systems of equations. Even in the absence of friction, several approaches exist for solving static contact problems [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Zhong | Static contact problems -a review[END_REF][START_REF] Wriggers | Finite element algorithms for contact problems[END_REF]. In most of them, the numerical methods that are employed for enforcing the contact constraints can be grouped into Lagrange multiplier and penalty methods [START_REF] Wriggers | Penalty and augmented Lagrangian formulations for contact problems[END_REF].

The penalty methods [START_REF] Kikuchi | Penalty/finite element approximations of a class of unilateral contact problems[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for dynamic contact / impact problems[END_REF] are closely related to the regularization of the contact constraints. They are usually formulated in terms of the displacement variables, and therefore are primal methods. They allow treating contact as a material behavior, as examplified by the method of joint finite elements [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]. Penalty methods can experience various numerical difficulties, especially illconditioning, when a too large or too small penalty parameter is introduced.

Lagrange multiplier methods are dual methods where the multipliers, which represent the contact reaction forces, are introduced in order to enforce exactly the non-penetration conditions. These methods can be further classified into (a) trial and error approaches for determining the active zone of contact, (b) mathematical programming techniques [START_REF] Klarbring | A mathematical programming approach to three-dimensional contact problems with friction[END_REF][START_REF] Hung | Frictionless contact of elastic bodies by finite element method and mathematical programming technique[END_REF], (c) search gradient methods [START_REF] May | The conjugate gradient method for unilateral problems[END_REF][START_REF] Dilintas | A conjugate projected gradient method with preconditioning for unilateral contact problems[END_REF][START_REF] Raous | Conjugate gradient for frictional contact[END_REF][START_REF]Proceedings of the Second Contact Mechanics International Symposium[END_REF], and (d) a combination of all of these approaches. Augmented Lagrange multiplier methods [START_REF] Arora | Multiplier methods for engineering optimization[END_REF][START_REF] Klarbring | Mathematical programming and augmented lagrangian methods for frictional contact problems[END_REF][START_REF] Radi | Treatment of the frictional contact via a Lagrangian formulation[END_REF][START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF] result in mixed formulations involving both displacement and force unknowns. The numerical solution schemes underlying both the Lagrange multiplier and augmented Lagrange multiplier methods are often related to the Uzawa algorithm [START_REF] Arrow | Studies in nonlinear programming[END_REF][START_REF] Simo | An augmented Lagrangian treatment of contact problem involving friction[END_REF][START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF].

Recently, various forms of domain decomposition have been introduced in the formulation of contact problems in order to speedup the solution of these problems, particularly when they are large-scale, on both serial and parallel computing platforms [START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF][START_REF] Dostal | Solution of coercive and semicoercive contact problems by FETI domain decomposition[END_REF]5]. Of particular interest is the approach proposed in [START_REF] Dostal | Solution of coercive and semicoercive contact problems by FETI domain decomposition[END_REF], not only because of its numerous innovative concepts, but also because it is based on the FETI domain decomposition method [START_REF] Farhat | A Lagrange multiplier based divide and conquer finite element algorithm[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF], which has well established numerical scalability properties for both second-order elasticity [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF][START_REF] Mandel | Convergence of a substructuring method with Lagrange multipliers[END_REF] and fourth-order plate and shell problems [START_REF] Farhat | The two-level FETI method for static and dynamic plate problems -Part I: an optimal iterative solver for biharmonic systems[END_REF][START_REF] Farhat | The two-level FETI method -Part II: extension to shell problems, parallel implementation and performance results[END_REF][START_REF] Mandel | A scalable substructuring method by Lagrange multipliers for plate bending problems[END_REF]. We remind the reader that an iterative method is said to be numerically scalable if its convergence properties -for example, the number of iterations for convergence -are asymptotically independent of the size of the problem to be solved.

The objective of this paper is to present an alternative domain decomposition method for the solution of frictionless contact problems that is also based on the FETI method, and which for this reason is named the FETI-C method (as in FETI-Contact).

FETI-C is in many aspects different from the domain decomposition method recently proposed in [START_REF] Dostal | Solution of coercive and semicoercive contact problems by FETI domain decomposition[END_REF] for solving frictionless contact problems. Indeed, the solution strategy described in [START_REF] Dostal | Solution of coercive and semicoercive contact problems by FETI domain decomposition[END_REF] is organized around two levels of iterations, the first one aimed at satisfying the contact conditions, and the second one at satisfying, among other things, the equilibrium conditions. In the outeriterations, an active zone of contact is updated by a mathematical programming technique. In the inner-iterations, a minimal subregion of the previously predicted area of contact is frozen and a state of equilibrium is sought after. In the FETI-C method described in this paper, both contact and equilibrium conditions are updated by a single iterative procedure. Most importantly, the FETI-C method incorporates an auxiliary "coarse contact problem" which not only guides the prediction of the active zone of contact, but also appears to ensure the numerical scalability of the proposed solution method with respect to both the number of subdomains, and the size of the problem. These numerical scalability properties are demonstrated numerically in this paper.

For simplicity, and without any loss of generality, we consider here only the case where the subdomains in contact have matching discrete interfaces.

We also assume that all structural components undergo small deformations as well as small displacements and rotations. In Section 2, we review the FETI method in order to keep this paper as self-contained as possible. In Section 3, we overview the formulation of the frictionless contact problem. In Section 4, we highlight the suitability of the FETI framework to the solution of contact problems, and describe the FETI-C non-linear solver. In Section 5, we validate this solver using sample two-and three-dimensional frictionless contact problems and demonstrate numerically its numerical scalability with respect to both the problem size and the number of subdomains. Finally, we conclude this paper in Section 6.

AN OVERVIEW OF THE FETI METHOD

Domain decomposition with Lagrange multipliers

The FETI method [START_REF] Farhat | A Lagrange multiplier based divide and conquer finite element algorithm[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Farhat | The two-level FETI method for static and dynamic plate problems -Part I: an optimal iterative solver for biharmonic systems[END_REF][START_REF] Farhat | The two-level FETI method -Part II: extension to shell problems, parallel implementation and performance results[END_REF] is a non-overlapping domain decomposition or substructuring method with Lagrange multipliers. It was originally developed for solving efficiently on parallel processors large-scale linear or linearized systems of equations arising from the discretization of structural mechanics problems defined on a domain Ω. Such a system of equations can be written as

K g u g = f g (1)
where for a large class of applications, K g is a symmetric positive definite or semi-definite matrix constructed from the assembly of element level stiffness matrices by introducing a global numbering of the degrees of freedom (d.o.f.) indicated here by the subscript g, u g is the vector of generalized displacements, and f g the vector of prescribed forces. In the FETI method, the given domain Ω is partitioned into N s non-overlapping subdomains Ω (s) , the equilibrium equations are written at the subdomain level, and Lagrange multipliers are introduced at the subdomain interfaces to enforce the compatibility of the displacement field.

In such an approach, the initial mechanical problem is essentially reformulated as a saddle point problem whose Euler equations are given by [START_REF] Farhat | A Lagrange multiplier based divide and conquer finite element algorithm[END_REF]

δ v (s) ,µ L(v (s) , µ) = 0 (2)
where δ v (s) ,µ designates the variation with respect to the arguments v (s) and µ, and L is the Lagrangian of the problem and can be written as

L(v (s) , µ) = Ns s=1 1 2 v (s) T K (s) v (s) -v (s) T f (s) + µ T Ns s=1 B (s) v (s) (3)
Here, the superscript s designates a subdomain quantity, the superscript T denotes the transpose of a quantity, v (s) and µ are some admissible subdomain displacement and Lagrange multiplier fields, respectively, K (s) is the subdomain u (1) u (2) λ

Ω (1) Ω (2)
Figure 1: Interface assembly in the FETI formulation stiffness matrix, f (s) the subdomain vector of prescribed forces, and B (s) the signed Boolean matrix that extracts from a subdomain vector v (s) its signed (±) restriction to the subdomain interface boundary. The Euler equations (2) lead to the following constrained problem which is equivalent to problem (1)

K (s) u (s) = f (s) -B (s) T λ s = 1, • • • , N s (4)
Ns s=1

B (s) u (s) = 0 (5)
where u (s) and λ are respectively the subdomain displacements and Lagrange multipliers that are the stationary points of L. The mechanical interpretation of equations ( 4) and ( 5) is graphically depicted in Figure 1: each of equations ( 4) expresses the subdomain equilibrium under the action of the prescribed forces f (s) and the subdomain interface forces B (s) T λ, and equation (5) expresses the compatibility of the subdomain displacement fields at the subdomain interfaces.

There is one compatibility equation for each pair of d.o.f. that connect at a subdomain interface. Once λ is determined, the subdomain displacement fields can be recovered by solving concurrently the equilibrium equations to obtain

u (s) = K (s) + f (s) -B (s) T λ + R (s) α (s) (6)
where K (s) + denotes the inverse of K (s) if Ω (s) has sufficient Dirichlet boundary conditions to prevent K (s) from being singular, or a generalized inverse of K (s) if Ω (s) is a floating subdomain -that is, a subdomain without sufficient displacement boundary conditions to prevent K (s) from being singular. In the latter case, the columns of R (s) represent the rigid body modes of Ω (s) , i.e. R (s) = ker K (s) , and α (s) is the set of amplitudes that specifies the contribution of the null space R (s) to the solution u (s) . These coefficients can be determined by requiring that each subdomain problem be mathematically solvable -that is, each floating subdomain be self-equilibrated -which can be written as

R (s) T f (s) -B (s) T λ = 0 (7)
Substituting equation ( 6) into the compatibility equation (5) and exploiting the solvability condition [START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF] transforms the constrained problem (4,5) into the interface problem

F I -G I -G I T 0 λ α = d -e (8)
where

F I = Ns s=1 B (s) K (s) + B (s) T d = Ns s=1 B (s) K (s) + f (s) G I = B (1) R (1) • • • B (Ns) R (Ns) (9)
α = α (1) T • • • α (Ns) T T e = f (1) T R (1) • • • f (Ns) T R (Ns)
T Because λ is a dual variable to the primal variables u (s) , the interface problem ( 8) is called a dual interface problem. This interface problem is best solved by an iterative algorithm, the choice of which completes the description of the original FETI method.

Iterative solution of the dual interface problem

In the FETI method, the interface problem ( 8) is solved by a preconditioned conjugate projected gradient (PCPG) algorithm. More specifically, the indefinite interface problem ( 8) is transformed into a semi-definite system of equations by eliminating the self-equilibrium condition G T I λ = e using the splitting

λ = λ 0 + P (Q)λ ( 10
)
where λ 0 is a particular solution of G T I λ = e of the form

λ 0 = QG I (G T I QG I ) -1 e (11)
and P (Q) is a projection matrix defined for any given matrix Q by

P (Q) = I -QG I (G T I QG I ) -1 G T I ( 12
)
Note that for any matrix Q, G T I P (Q) = 0. Applying the splitting [START_REF] Dilintas | A conjugate projected gradient method with preconditioning for unilateral contact problems[END_REF] to the interface problem (8) leads to the alternative symmetric positive semi-definite interface problem

P (Q) T F I P (Q)λ = P (Q) T d -F I λ 0 (13)
which can be solved by a standard PCG algorithm. Two different preconditioners have been proposed: the mathematically optimal Dirichlet preconditioner F

D -1 I
introduced in [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF], and the computationally economical lumped preconditioner F

L -1 I
proposed in earlier works [START_REF] Farhat | A Lagrange multiplier based divide and conquer finite element algorithm[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. If each subdomain stiffness matrix is partitioned as

K (s) = K (s) ii K (s) ib K (s) T ib K (s) bb (14)
where the subscripts i and b designate the subdomain interior and interface boundary d.o.f., respectively, then the Dirichlet preconditioner can be written as

F D -1 I = W ( Ns s=1 B (s) 0 0 0 S (s) bb B (s) T )W ( 15
)
and the lumped preconditioner is given by

F L -1 I = W ( Ns s=1 B (s) 0 0 0 K (s) bb B (s) T )W (16)
In the above expressions of F

D -1 I and F L -1 I
, W is a diagonal matrix storing in each of its entries the inverse of the multiplicity of an interface d.o.f. [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Farhat | The two-level FETI method for static and dynamic plate problems -Part I: an optimal iterative solver for biharmonic systems[END_REF] -that is, the inverse of the number of subdomains to which an interface d.o.f. belongs -and S (s) bb is the following subdomain primal Schur complement

S (s) bb = K (s) bb -K (s) T ib K (s) -1 ii K (s) ib (17) Note that K (s)
ii is non-singular since it corresponds to the system matrix with all interface boundaries fixed and thus K (s) -1 ii exists. Both Dirichlet and lumped preconditioners have been recently extended in [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF] for addressing more efficiently heterogeneous problems. In particular, it was shown in [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF] that such extensions can be simply designed by redefining appropriately the scaling matrix W .

The iterative FETI solver with Q = I is summarized in Table 1. A more detailed description of this algorithm can be found in [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF].

Remarks

• For most problems, the simplest choice

Q = I is computationally the most effective one. Choosing Q = Q D = F D -1 I or Q = Q L = F L -1 I
is recommended for heterogeneous problems [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF].

Initialize λ 0 = GI (GI T GI ) -1 e r 0 = d -FI λ 0 Iterate k = 1, 2, • • • until convergence Project w k-1 = P T r k-1 Precondition z k-1 = F -1 I w k-1 Re-project y k-1 = P z k-1 Conjugate ζ k = y k-1 T w k-1 y k-2 T w k-2 (ζ 1 = 0) p k = y k-1 + ζ k p k-1 (p 1 = y 0 ) Minimize η k = p k-1 T w k-1 p k T FI p k Update λ k = λ k-1 + η k p k Compute residual r k = r k-1 -η k FI p k
• From Eqs. [START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF][START_REF]Proceedings of the Contact Mechanics International Symposium[END_REF][START_REF] Farhat | A Lagrange multiplier based divide and conquer finite element algorithm[END_REF] and Table 1, it follows that Ns s=1

B (s) u (s) = d -F I λ + G I α = r + G I α = P T r = w (18)
which shows that at each iteration k, the projected residual w k-1 is equal to the jump of the (k -1)-th iterate displacement field across the subdomain interfaces.

• At each FETI PCPG iteration, a matrix-vector product of the form P T r k-1 or P z k-1 incurs the solution of an auxiliary problem of the form

(G T I G I )x k-1 = b k-1 (19)
Hence, the matrix G T I G I (and more generally the matrix G T I QG I ) defines a coarse problem of size at most equal to 6N s , and whose solution couples all the subdomain computations, propagates the error globally, and accelerates convergence. It is because of this coarse problem that for second-order elasticity problems, when the Dirichlet preconditioner is used, the condition number of the FETI interface problem (13) can be bounded by a polylogarithmic function of H/h as follows

κ = O(1 + log 2 ( H h )) (20)
where H and h denote the subdomain and mesh sizes, respectively. The above conditioning result also holds for fourth-order plate and shell problems when the coarse problem ( 19) is augmented with the so-called subdomain "corner" modes [START_REF] Farhat | The two-level FETI method for static and dynamic plate problems -Part I: an optimal iterative solver for biharmonic systems[END_REF][START_REF] Farhat | The two-level FETI method -Part II: extension to shell problems, parallel implementation and performance results[END_REF][START_REF] Mandel | A scalable substructuring method by Lagrange multipliers for plate bending problems[END_REF]. We remind the reader that a bound such as [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF] proves the numerical scalability of the FETI method with respect to both the problem size and the number of subdomains [START_REF] Mandel | Convergence of a substructuring method with Lagrange multipliers[END_REF][START_REF] Mandel | A scalable substructuring method by Lagrange multipliers for plate bending problems[END_REF] -that is, its ability to solve larger problems using a larger number of subdomains in almost a constant number of iterations. The parallel scalability of the FETI method -that is, its suitability for massively parallel computing -has also been demonstrated on a large number of massively parallel processors for several realistic structural and structural dynamics problems [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Farhat | The two-level FETI method -Part II: extension to shell problems, parallel implementation and performance results[END_REF].

Relevance to contact problems

The subdomains Ω (s) associated with the FETI method can be viewed as substructures in a state of "perfect" contact. This perfect contact is defined and governed by the compatibility equation ( 5), which states that the jumps across the substructure interfaces of both the normal and tangential components of the displacement field must vanish. In the FETI method, this perfect contact condition is enforced by discrete Lagrange multipliers whose values are determined by the PCPG algorithm summarized in Table 1.

The above observations explain our motivation to extend the FETI method to the solution of contact problems whose governing laws are summarized next.

THE FRICTIONLESS CONTACT PROBLEM

Let Ω denote the initial state of a linear elastic body, ∂Ω C the part of its boundary that is a potential area of contact, n C the outward normal to ∂Ω C at a point M , ∂Ω ′ C the potential area of contact with ∂Ω C located on the surface of a foundation or another linear elastic body, and n 2). If u denotes the displacement field, the discontinuity of u in the normal direction at point M ∈ ∂Ω C can be written as

′ C = -n C the outward normal to ∂Ω ′ C at the point M ′ facing point M (see Figure
∆u(M ) = u(M ) • n C + u(M ′ ) • n ′ C = u(M ) • n C -u(M ′ ) • n C = (u(M ) -u(M ′ )) • n C (21)
and the non-penetration condition can be expressed as

∆u(M ) -c(M ) 0 ≤ 0 ( 22
)
where c(M ) 0 denotes the initial clearance at point M of the undeformed body Ω (see Figure 3), and is assumed to be small. Let σ denote the stress tensor, t = σ(M )n the traction vector at point M , and t n = t • n the normal component of t. If ∂Ω C and ∂Ω ′ C are in contact, the pressure at M is positive, and therefore t n ≤ 0 (Figure 4). The contact problem is defined by the unilateral conditions known as the Signorini-Fichera type of conditions

t n = 0 t n ≤ 0 M M' ∆u = c 0 ∆u ≤ c 0 n' c n c ∂Ω c ∂Ω' c Figure 2: A contact problem M M' n c ∂Ω c ∂Ω' c c 0 (M) Figure 3: Initial clearance M n c ∂Ω c t = σ(M)n t n = t.n
• if ∆u(M ) -c(M ) 0 = 0 then t n ≤ 0 and the contact is said to be active.

• if ∆u(M ) -c(M ) 0 < 0 then t n = 0 and the contact is said to be inactive.

(

) 23
The above relations are equivalent to stating that on ∂Ω C , the following holds

t n ≤ 0 ∆u(M ) -c(M ) 0 ≤ 0 (24) t n • (∆u(M ) -c(M ) 0 ) = 0

THE FETI-C METHOD

Our objective is to develop a computationally efficient domain decomposition method for the solution of contact problems that is based on FETI. In this paper, we limit however our attention to frictionless problems. For most domain decomposition methods including FETI, the optimal number of subdomains is dictated by a compromise between several issues related to, among others, computational complexity, memory efficiency, and parallel processing. For this reason, we consider the general case where the number of subdomains N s is larger than the given number of elastic bodies and foundations -that is, the case where each elastic body itself may be decomposed into several subdomains. This leads us to distinguishing between two types of subdomain interfaces: perfect contact interfaces (called perfect interfaces in the sequel), and potential contact interfaces. We construct the FETI-C method as a modified FETI algorithm that (a) distinguishes between perfect and potential contact interfaces, and (b) enforces the compatibility equation (5) on perfect interfaces and the non-penetration condition [START_REF] Klarbring | Mathematical programming and augmented lagrangian methods for frictional contact problems[END_REF] as well as the Signorini-Fichera conditions [START_REF] Mandel | Convergence of a substructuring method with Lagrange multipliers[END_REF] on contact interfaces, while (c) satisfying equilibrium and solvability in each subdomain Ω (s) .

Preliminaries

Let ∂Ω (s),p P = ∂Ω (s) ∩ ∂Ω (p) (25)
denote a perfect interface between subdomain Ω (s) and neighboring subdomain Ω (p) , ∂Ω

(s),q C = ∂Ω (s) ∩ ∂Ω (q) ( 26
)
denote a potential contact interface between Ω (s) and its neighbor Ω (q) , ∂Ω

P = p=l p=1 ∂Ω (s),p P (s)
denote the union of all the perfect interfaces attached to subdomain Ω (s) , ∂Ω

C = q=m q=1 ∂Ω (s),q C (28) (s)
denote the union of all the potential contact interfaces attached to subdomain Ω (s) ,

Γ P = s=Ns s=1 ∂Ω (s) P (29)
denote the set of all the perfect interfaces associated with the domain decomposition,

Γ C = s=Ns s=1 ∂Ω (s) C (30)
denote the set of all potential contact interfaces and is obtained from the initial configuration of the bodies, and

Γ I = s=Ns s=1 ∂Ω (s) = Γ P ∪ Γ C (31)
denote the global interface associated with the domain decomposition of all given bodies. Using this notation, the edge-wise partition of the interface boundary ∂Ω (s) of subdomain Ω (s) can be written as

∂Ω (s) = ∂Ω (s) P ∪ ∂Ω (s) C = { p=l p=1 ∂Ω (s),p P } ∪ { q=m q=1 ∂Ω (s),q C } ( 32
)
Let N I denote the total number of kinematic equations defined on Γ I for the given contact problem, N s,p denote the total number of displacement d.o.f. attached to the perfect interface ∂Ω (s),p P , and N s,q the total number of displacement d.o.f. attached to the potential contact interface ∂Ω (s),q C . To each perfect interface ∂Ω

where v (s) is a generic vector defined in subdomain Ω (s) , v (s),p P denotes its trace on the perfect interface ∂Ω (s),p P , and the sign (±) is defined by the orientation of the outward normal to ∂Ω (s),p P . Hence, b (s),p P is a signed Boolean matrix. For example, if all the interior d.o.f. of the subdomain Ω (s) are numbered first, and all the interface boundary d.o.f. of this subdomain are numbered last according to the edge-wise partition [START_REF] Wriggers | Finite element algorithms for contact problems[END_REF], the signed Boolean matrix B (s) introduced in Section 2.1 for the case without contact can be written as

B (s) = [0 B (s) I ] ( 34
)
where

B (s) I = [b (s),1 P • • • b (s),l P ] (35)
To each potential contact interface ∂Ω (s),q C

, we associate the matrix b

(s),q C ∈ R N I ×Ns,q defined by (b (s),q C v (s),q C ) M = (n (s),q T C ) M (v (s),q C ) M on ∂Ω (s),q C 0 elsewhere on Γ I ( 36
)
where n (s),q C

denotes the vector of outward normals to ∂Ω

(s),q C at each point M ∈ ∂Ω (s),q C
, and the subscript M designates the subvector or entry associated with node M . In Eq. ( 36) and throughout this paper, the subscript M also implies that the stated relation holds at each relevant point M . Note that b (s),q C is not a Boolean matrix, except when ∂Ω (s),q C is parallel to a reference axis. Hence, for contact problems, we generalize the matrix B (s) introduced in Eq. ( 3) to account for the zones of potential contact, and therefore transform Eq. ( 35) into

B (s) I = [b (s),1 P • • • b (s),l P ] [b (s),1 C • • • b (s),m C ] ( 37
)
At each node M ∈ ∂Ω (s),q C

, and for each interface ∂Ω (s),q ′ C that has the potential of getting in contact with ∂Ω (s),q C , we introduce one potential Lagrange multiplier along the normal direction n

(s),q C to ∂Ω (s),q C
. Such a Lagrange multiplier represents a contact force in the direction n (s),q C . In the case of contact, the governing subdomain equations (4,5) become

K (s) u (s) = f (s) -B (s) T λ s = 1, • • • , N s (38) Ns s=1 B (s) u (s) -c 0 ≥ 0 on Γ I ( 39
)
where B (s) differs from B (s) only when at least one zone of potential contact is inactive (see next section for further details), and c 0 denotes the extended vector of initial clearances -that is, c 0 is defined on both Γ P and Γ C but set to zero on the perfect interfaces. Note that this definition of c 0 is chosen so that Eq. (39) contain the expected compatibility relation on the perfect interfaces

Ns s=1 B (s) u (s) = 0 on Γ P (40)
Ω ( 4)

Ω (1)
Ω ( 2)

Ω (3) ∂Ω c (1),2 ∂Ω c (1),4 ∂Ω c (2),1 λ n c
(1),2

Ω (4) Ω (1)
Ω ( 2)

Ω (3)
∂Ω p (1),2

∂Ω p (1),4

∂Ω p (2),1 λ From Eq. (39) and Eq. ( 38), it follows that for contact problems, the jump w of the displacement field across the subdomain interfaces (see Eq. ( 18)) is

w = Ns s=1 (B (s) K (s) + f (s) -B (s) K (s) + B (s) T λ) + Gα -c 0 (41)
In summary, using the notation introduced in this section, we can re-write the contact conditions [START_REF] Mandel | A scalable substructuring method by Lagrange multipliers for plate bending problems[END_REF] on each potential contact interface ∂Ω (s),q C in terms of FETI variables as follows

λ M ≥ 0 w M ≤ 0 (42) λ M w M = 0

Remark

• At a crosspoint between four subdomains, six Lagrange multipliers are usually introduced by the standard FETI method to connect together what has been defined in this paper as perfect interfaces. However, from the presentation made above, it follows that at a crosspoint between four subdomains characterized by four potential contact interfaces, only four Lagrange multipliers can be introduced to represent the potential contact forces (Figure 5).

Next, we lay the foundations of the FETI-C method.

Discrete maximization problem

The standard FETI dual interface problem (8) can be reformulated as follows max

µ∈R N I , β∈R N R C(µ, β) = - 1 2 µ T F I µ + µ T d + β T (G T I µ -e) (43)
where C(µ, β) is a complementary energy-like functional, N I has the same meaning as in the previous subsection, and N R is the total number of subdomain rigid body modes. Indeed, differentiating C(µ, β) with respect to µ and β and computing (λ, α) as the solution of ( ∂C ∂µ = 0, ∂C ∂β = 0) leads to the same equations as in [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF].

From Eqs. (42), it follows that the standard FETI maximization problem (43) must be modified as follows in order to address contact problems

max e Γ C , µ∈R N I , µ M ≥0, β∈R N R C(µ, β) = - 1 2 µ T F I µ + µ T d + β T ( G I T µ -e) (44)
where

F I = Ns s=1 B (s) K (s) + B (s) T d = Ns s=1 B (s) K (s) + f (s) -c 0 (45)
G I = B (1) R (1) • • • B (Ns) R (Ns)
Γ C is the set of active contact interfaces -that is, the interfaces where µ M > 0 -and therefore

Γ C verifies Γ C ⊂ Γ C (46)
B (s) has the same pattern as in Eq. ( 34), B

I (s)
is given in Eq. (37),

B (s) = 0 B (s) I (47)
B (s) I = [b (s),1 P • • • b (s),l P ] [ b(s),1 C • • • b(s),m C ] (48)
and

( b(s),q C v (s),q C ) M = (b (s),q C v (s),q C ) M if M ∈ Γ C 0 otherwise ( 49
)
Let w be defined as follows

w = w on Γ P ∪ Γ C w M + on Γ C \ Γ C ( 50
)
where w is the jump of the displacement field across the subdomain interfaces as expressed in Eq. ( 41), and

a + = a if a > 0 0 otherwise (51)
From the second of Eqs. (42), it follows that if λ is solution of the maximization problem (44), and the corresponding subdomain displacement fields u (s) computed from the equilibrium equations (38) satisfy w = 0, then (λ, {u (s) } s=Ns s=1 ) is the domain decomposed solution of the given contact problem.

The FETI-C method we introduce next is a FETI-like method which finds the zone of active contact Γ C and computes the elastic state of each substructure, by solving the maximization problem (44) that embeds the equilibrium equations (38), while ensuring that w = 0.

Remarks

• The operator F I introduced in Eq. ( 45) is unsymmetric. However, from Eqs. (37,48,49), it follows that the set of admissible Lagrange multipliers on which F I operates is such that

B T λ = B T λ (52)
and therefore a PCG algorithm can be applied to the solution of the maximization problem (44).

• From Eqs. (9,45) and Eq. ( 52), it also follows that

G I (s) T λ = G (s) T I λ (53)

One-shot iterative procedure

As highlighted in the previous sections, the discrete maximization problem (44) associated with contact problems has many similarities with its counterpart (43) for standard structural mechanics problems. For this reason, a good starting point for solving problem (44) is the FETI PCPG algorithm summarized in Table 1. However, a major difference between the two maximization problems (44) and ( 43) is that the former one is non-linear, whereas the latter one is linear. Another difference pertains to the admissibility of the Lagrange multipliers.

For standard structural mechanics problems, the Lagrange multipliers must satisfy the self-equilibrium condition [START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF] in each floating subdomain. For contact problems, the Lagrange multipliers must also satisfy the additional positivity condition stated in the first of Eqs. (42). Hence, our idea here is to construct a non-linear solution algorithm for solving the discrete maximization problem (44) that mimics the FETI PCPG linear solver, while addressing the differences between perfect and contact interfaces. The resulting non-linear scheme, FETI-C, is a FETI-like iterative procedure where both the contact and equilibrium conditions are updated simultaneously. Hence, FETI-C does not involve any linearization and therefore unlike the standard FETI solver summarized in Table 1, it cannot by interpreted as a classical PCG algorithm. This new extension of the FETI method is presented in the following subsections in a step-by-step approach.

Initialization step.

Consider the FETI-PCPG algorithm summarized in Table 1. Suppose that the initialization step

λ 0 = G I (G T I G I ) -1 e (54)
is performed using

Γ 0 C = Γ C (55)
In that case, G I denotes an initialization of G I (45) where each subdomain matrix B (s) I (48) is constructed by assuming that all potential contact interfaces are initially active. In other words, G T I G I is based on Γ P and Γ C . If the given problem does not have any global rigid body mode, G T I G I is non-singular [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. On the other hand, if the given problem has one or several global rigid body modes, G T I G I is singular, but Eq. ( 54) admits at least one solution if the external loads are self-equilibrated. In any case, the initial value λ 0 as computed in Eq. ( 54) is not guaranteed to be admissible for the maximization problem (44) -that is, λ 0 is not guaranteed to satisfy λ 0 M > 0 on Γ C . For this reason, we propose to apply after the initialization step a "planing" procedure aimed at enforcing the positiveness of λ 0 in the sense defined in the first of Eqs. (42), without for that purpose destroying the self-equilibrium of any floating subdomain (see Eq. ( 7)). This planning step is represented by the following transformation λ ←-P(λ)

where P denotes the two-step iterative procedure described in Table 2. In each first step of the n-th iteration of this planing procedure, the positiveness of the iterate λ n is enforced simply by setting to zero the negative Lagrange multipliers acting on the zone of contact Γ C , and the active zone of contact Γ C and matrices B , respectively. In each second step, the self-equilibrium of the floating subdomains is restored by applying to λ n the following projector

P n = I -G I n ( G I n T G I n ) -1 G n T I (57)
and translating the result by the constant vector

G I n ( G I n T G I n
) -1 e. Indeed, from Eq. (57) and Eq. ( 53), it follows that

G I n T ( P n λ n + G I n ( G I n T G I n ) -1 e) = ( G I n T -G n T I )λ n + e = e ( 58
)
Note that the projector (57) is similar to the projector [START_REF] Farhat | A Lagrange multiplier based divide and conquer finite element algorithm[END_REF] equipped with Q = I, but differs from it in the fact that it accounts for the change in the

) n I , f GI n 1.(c) Test self-equilibrium: if f GI n T λ n = e end.
2.(a) Self-equilibrium:

λ n ←-(I -f GI n ( f GI n T f GI n ) -1 G n T I )λ n + f GI n ( f GI n T f GI n ) -1 e 2.(b) Test positiveness: if ∃M ∈ ΓC / λ n M < 0 go to 1.(a)
active zone of contact from Γ n-1 C to Γ n C . Experience shows that in general, the planing procedure described in Table 2 converges after a small number of planing subiterations (see Section 5).

Note also that Step 2.(a) of the planing procedure incorporates a coarse problem of the form

( G I n T G I n )x n = b n (59)
This coarse problem is called here a coarse "contact" problem because, as shown in Table 2, it influences the updating of the active zone of contact.

Evaluation of the projected residual.

As noted in Section 4.2, for contact problems, it is the modified projected residual w defined in Eq. ( 50) and not the projected residual w (18) that must be driven to zero. For this reason, at each iteration k, after w k-1 = P k-1 T r k-1 is evaluated, wk-1 is computed and the preconditioner is applied to wk-1 rather than w k-1 . For the same reason, wk-1 and wk-2 are used in FETI-C wherever w k-1 and w k-2 are used in FETI. The re-projection step is similar to that of the standard FETI PCPG solver and can be written as

y k-1 = P k-1 z k-1 , where z k-1 = F I -1
wk-1 , and F I -1 is the preconditioner discussed next.

Preconditioning step.

The straightforward extension to contact problems of the Dirichlet (15) preconditioner calls for substituting the matrices B (s) in Eqs. [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] by the matrices B (s) k associated with the k-th iterate active zone of contact Γ k C , and updating the definition of the interior (subscript i) and interface boundary (subscript b) d.o.f. accordingly. However, in order to avoid non-desirable oscillations in the enforcement of the non-penetration condition (second of Eqs. ( 42)), we propose to precondition only the equations corresponding to the d.o.f. where the modified projected residual w is non-zero. Hence, we propose the following "Dirichlet-like" preconditioner bb is computed once based on the initial configuration -that is, on Γ P ∪ Γ 0 C = Γ P ∪ Γ C -and is frozen throughout all iterations. This strategy is chosen to maximize computational efficiency.

F D -1 I k = W ( Ns s=1 B (s) k 0 0 0 S (s) bb B (s) k T )W (
Note that the fact that the above preconditioner changes at each iteration k is a direct consequence of the choice we have made to design a non-linear FETI-C solution algorithm that does not require any linearization of the non-linear maximization problem (44).

Conjugation and maximization steps.

The derivation and explanation of the conjugation and maximization (or minimization) steps that characterize a PCG algorithm are usually performed for a linear problem. In general, these steps are not required for a non-linear problem. Nevertheless, we enforce them in the FETI-C solver in order to avoid as much as possible numerical oscillations. However, in this case and even in the absence of round-off errors, the conjugation step enforces the orthogonality of the current search direction only to the previous one.

Planing step.

At each iteration k, the planing procedure (56) (or planing step) is applied to the updated Lagrange multipliers in order to enforce both the positiveness of the contact pressure and the self-equilibrium of the floating subdomains. Note that at the end of each k-th application of the planing procedure, a possibly new active zone of contact Γ k+1 C is predicted, and therefore is available for updating the matrices B (s) to become B (s) k , and the projector P and the precondititioner , respectively.

Evaluation of the residual.

In the sequel, we define a "status change" as the event where any node on the potential zone of contact Γ C changes status from active to inactive, or viceversa, during a FETI-C iteration. A status change affects the evaluation of the residual. Indeed, we first note that the residual can be computed only after the planing step, and that after each planing step the predicted active zone of contact is updated. If a status change occurs during this update, F I is modified, which prevents evaluating the residual as in a classical CG algorithm. Instead, if a status change occurs, this residual but must be re-computed as follows

r k = d -F I k+1 λ k (61)
where the notation F I k+1 means F I k+1 = F I ( Γ k+1 C ) and emphasizes the fact that F I depends on the predicted active zone of contact Γ C via the matrices B (s) (see Eqs. (45-49)). On the other hand, if no status change occurs during the updating of the active zone of contact, the residual can be computed as in the classical CG algorithm.

The FETI-C solver

The FETI-C solver for frictionless contact problems is summarized in Table 3. Convergence is checked by monitoring wk-1 and declared when

wk-1 ≤ ǫ ( Ns s=1 B (s) K (s) + f (s) + c 0 ) ( 62
)
where ǫ is a user-specified tolerance, and the triangular inequality has been applied to d (45) to avoid potential numerical difficulties.

NUMERICAL EXAMPLES

In this section, we assess the performance of the FETI-C solver and investigate its numerical scalability properties. For this purpose, we apply this algorithm, with and without preconditioning, to the solution of several two-and threedimensional frictionless contact problems. In all cases we monitor its convergence as described in (62) with ǫ = 10 -7 . Furthermore, we use the performance of the optimal FETI method as a reference to evaluate the numerical efficiency of the FETI-C solver, in the following manner. For each contact problem, after the FETI-C solver has converged and the active zone of contact has been determined, we "lock" the substructures in their final configuration -that is, we consider each determined active contact interface as a perfect interface in the normal direction -and perform a standard stress analysis of that configuration under the given external loads, using the standard FETI method as a solver. We compare the dependence on a few key parameters of the number of FETI iterations associated with this stress analysis to that of the number of FETI-C iterations performed during the solution of the underlying contact problem.

Initialize e Γ 0 C = ΓC λ0 = GI (GI T GI ) -1 e λ 0 = P( λ0 ) e P 0 = e P ( e Γ 1 C ) r 0 = d -e FI 1 λ 0 Iterate k = 1, 2, • • • until convergence Project w k-1 = e P k-1 T r k-1 wk-1 =  w k-1 on ΓP w k-1 M + on e Γ k C Precondition z k-1 = e F -1 I k wk-1 Re-project y k-1 = e P k-1 z k-1 Conjugate ζ k = y k-1 T wk-1 y k-2 T wk-2 (ζ 1 = 0) p k = y k-1 + ζ k p k-1 (p 1 = y 0 ) Maximize η k = p k-1 T wk-1 p k T e FI k p k Update λk = λ k-1 + η k p k Apply planing λ k = P( λk ) e P k = e P ( e Γ k+1 C ) Compute residual r k = ( r k-1 -η k e FI k p k if no satus change d -e FI k+1 λ k if status change x y -Fx -F(x -y)/2 -F(x + y)/2

Two-dimensional applications

Here, we consider a series of two-dimensional contact problems between 2 to 6 substructures. Each one of this problem features a different challenging aspect of contact. We assume that all substructures are made of the same linear elastic material characterized by a Young modulus E = 2.05 10 3 MPa, and a Poisson ratio ν = 0.3. We discretize all substructures by four-noded plane stress elements, design small meshes with 1,364 d.o.f. as well as larger ones with up to 175,692 d.o.f., and vary the total number of subdomains between 6 and 150.

A multi-block problem.

First, we consider the problem graphically depicted in Figure 6, with F = 10 4 N.

In this problem, 6 blocks sharing 7 contact interfaces are pushed against each other by a set of external forces. The two simple supports drawn in one of the blocks mean that the interface boundary nodes of that block can move only in the horizontal (x) direction. This condition ensures that the contact problem is well-posed and has a unique solution.

We construct three different domain decompositions with 6 (H = 1/2), 54 (H = 1/6), and 150 (H = 1/10) subdomains (Figure 7). In each subdomain, we generate a suite of (1/h) × (1/h) uniform meshes with 10 ≤ (H/h) ≤ 60. For each domain decomposition and each mesh size, we report in Table 4 the performance results of the FETI-C solver. These consist of the number of FETIlike iterations N itr , the total number of planing subiterations N pln , and the total number of status changes N stc . We denote by N dof the total number of d.o.f. associated with a given mesh, and by N lock itr the number of FETI iterations obtained for the solution of the locked configuration of this contact problem with Dirichlet preconditioner.

The performance results summarized in Table 4 show that

• The Dirichlet preconditioner reduces the number of FETI-C iterations by a factor ranging between 2 and 3.

subdomains

54 subdomains 150 subdomains • The FETI-C solver equipped with the Dirichlet preconditioner is numerically scalable with respect to both the problem size and the number of subdomains.

• For this contact problem, the FETI-C solver equipped with the Dirichlet preconditioner performs an amazingly small number of planing subiterations. Furthermore, for large-scale meshes and mesh partitions, it requires a number of iterations that is comparable to that of the FETI method applied to the standard stress analysis of the locked configuration of this problem.

The third observation formulated above can be explained by Figure 8 where the displacement solution of this contact problem using the coarsest mesh is displayed on the left, and the distribution of the contact Lagrange multipliers is graphically depicted using arrows on the right. In this figure, the reader can observe that almost all contact interfaces are active, which explains why few planing subiterations are required and N itr is close to N lock itr . The deformed shape shown on the left in Figure 8 also illustrates the reason why, at a crosspoint, no Lagrange multiplier should be introduced between two subdomains that share only one node in two-dimensions, or on edge in three-dimensions (see the remark in Section 4.1).

A contact problem with pre-stress.

Next, we consider a frictionless contact problem where a square elastic body is forced into the cavity of another U-shaped elastic body, as shown in Figure 9 with c 0 = -10 -6 m and F = 100 N. This problem has only 3 potential contact interfaces. However, it has the same topology as the previous problem. For this reason, we use the same meshes and same domain decompositions as for the previous problem (see Figure 7). We report in Table 5 and Figure 10 the performance results obtained for the FETI-C solver.

The reported values of N lock itr suggest that, even for a standard linear configuration, this problem is intrinsically "harder" to solve by a FETI method than the previous one. This explains why the FETI-C method performs more iterations for solving this problem than for solving the previous one. Nevertheless, the performance results summarized in Table 5 also show that the FETI-C solver equipped with the Dirichlet preconditioner: (a) is numerically scalable with respect to both the problem size and the number of subdomains, and (b) performs on average only 1 or 2 planing subiterations per iteration. The latter observation suggests that the coarse contact problem ( 59) is effective at reducing the oscillations that are usually observed when determining the active contact zone by an iterative procedure. Figure 11 reports on the evolution of the relative projected residual (62), number of planing subiterations, and number of nodes that experience a status change throughout the iterations, for the case with 54 subdomains, H/h = 40, and the Dirichlet preconditioner. The reader can observe that after a few iterations, no status change occurs and therefore the evaluation of the residual becomes as economical as in a classical CG algorithm. The solution of this contact problem with pre-stress is graphically depicted in Figure 12. The gap between the deformed interfaces of the two elastic bodies is due to the fact that the binding c 0 has been prescribed with respect to the

A three-dimensional contact problem

In order to test the behavior of the planing procedure (56) on two-dimensional interfaces, we finally consider the three-dimensional frictionless contact problem shown in Figure 13. In this problem -which is a three-dimensional version of the contact problem described in Section 5.1.1 -the translation in the z direction is restrained for all of substructures 4,5, and 6, and the translation in the y direction as well as the rotation around the x axis are restrained for substructure 5. The loading is identical to that of the multi-block problem shown in Figure 6, but only 80% of the potential contact nodes turn out to be active in this problem.

We generate two different meshes corresponding to discretizing each substructure by 4 × 4 × 4 and 12 × 12 × 12 8-noded brick elements. Hence, the first mesh contains 2, 250 d.o.f., and the second one 39, 546 d.o.f. We treat each substructure as one subdomain. We report in Table 6 the performance results of the FETI-C solver applied to this contact problem. Once again, these results show that the FETI-C solver equipped with the Dirichlet preconditioner performs on average about 2 planing subiterations per iteration, and is scalable with respect to both the size of the problem and the number of subdomains.

On the solution of the coarse contact problem

The implementation on sequential and parallel computing platforms of the FETI-C solver presented in this paper is similar to that of the FETI method. It will be discussed in details in a forthcoming paper, which will also report on CPU performance results. Here, we make only one comment on the solution of the coarse contact problem (59). When the number of potential contact interfaces can be expected to be a small fraction of the total number of interfaces -for example, when only a few bodies are in contact and a large number of subdomains is employed -we recommend that the contribution of the perfect interfaces to ( G I n T G I n ) be eliminated once for all by a static condensation procedure, before starting the FETI-C iterations. In that case, at each planing subiteration of a given FETI-C iteration, the cost of solving the coarse contact problem (59) becomes negligible compared to the other computational costs of the FETI-C method.

SUMMARY AND CONCLUSIONS

The FETI-C method is a FETI-based substructuring method with Lagrange multipliers aimed for the iterative solution of general contact problems. As a first step, it is developed in this paper for frictionless problems. It is organized around a single iteration loop that updates both the contact and equilibrium conditions. A distinctive feature of the FETI-C method is an iterative planing procedure which enforces at each FETI-C iteration the positiveness of the contact pressure while preserving the self-equilibrium of each floating subdomain. This planing procedure incorporates an auxiliary coarse problem which not only guides the prediction of the active zone of contact, but most importantly appears to be effective at reducing the oscillations that are usually observed when determining the active contact zone by an iterative procedure. The computational cost of each planing subiteration is a minor fraction of the computational cost of a FETI-C iteration. Furthermore, extensive numerical experiments reveal that the FETI-C method performs on average less than 2 planing subiterations per iteration. Extensive numerical experiments with two-and three-dimensional frictionles contact problems also suggest that the FETI-C method is numerically scalable with respect to both the number of subdomains, and the size of the problem. The extension of this method to contact problems with friction necessitates (a) modifying the admissibility condition on the Lagrange multipliers introduced at the contact interfaces to account for friction, and (b) superposing to the current residuals additional terms that account for the tangential components of the jumps of the displacement field. These issues will be dealt with in a forthcoming paper.

Figure 4 :

 4 Figure 4: Contact pressure: t n ≤ 0

  (s),p P , we associate the matrix b (s),p P ∈ R N I ×Ns,p defined by b

Figure 5 :

 5 Figure 5: Corner Lagrange multipliers: perfect interfaces (left), contact interfaces (right)

  and denoted by Γ n C and B (s) n I

  k is the restriction of B (s) k to the d.o.f. of Γ P and to the d.o.f. of Γ k C where the effective residual wk-1 is non-zero. For simplicity, we refer in the sequel to the Dirichlet-like preconditioner introduced above as the Dirichlet preconditioner. At each iteration k, the above Dirichlet preconditioner changes because the d.o.f. on Γ C where wk-1 = 0 change. However, only the matrices B (s) are updated. The Schur complement S (s)

F - 1 I

 1 to become P k and F

Figure 6 :

 6 Figure 6: A 6-block problem with 7 contact interfaces

Figure 7 :

 7 Figure 7: Domain decomposition

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Multi-block contact problem: deformed interfaces (left) and corresponding Lagrange multipliers (right, arrows)

Figure 11 :Figure 12 :

 1112 Figure 11: Convergence behavior of the FETI-C solver for the contact problem with a pre-stressed block

Figure 13 :

 13 Figure 13: A three-dimensional 6-block contact problem

Table 1 :

 1 The FETI PCPG algorithm

Table 2 :

 2 The planing procedure

	Iterate n = 1, 2, • • • until convergence 1.(a) Positiveness: λ n ←-		λ n λ n M +	on ΓP on ΓC
	1.(b) Active contact:	determine e Γ n C from λ n and update e B	(s

Table 3 :

 3 The FETI-C non-linear solution algorithm

Table 4 :

 4 Multi-block contact problem: performance results of the FETI-C solver

	H(Ns)	H/h	N dof	preconditioner Nitr N pln Nstc N lock itr
	1/2 (6)	10	1,452	I	23	0	0
	1/2 (6)	20	5,292	I	31	2	2
	1/2 (6)	40	19,200	I	43	4	3
	1/2 (6)	60	44,652	I	50	4	3
	1/6 (54)	10	11,532	I	52	3	2
	1/6 (54)	20	44,652	I	68	6	3
	1/6 (54)	40	175,692	I	88	11	5
	1/10 (150)	10	31,212	I	57	8	5
	1/2 (6)	10	1,452	Dirichlet	8	0	0	7
	1/2 (6)	20	5,292	Dirichlet	11	0	0	7
	1/2 (6)	40	19,200	Dirichlet	13	0	0	6
	1/2 (6)	60	44,652	Dirichlet	13	1	1	7
	1/6 (54)	10	11,532	Dirichlet	30	0	0	24
	1/6 (54)	20	44,652	Dirichlet	34	0	0	26
	1/6 (54)	40	175,692	Dirichlet	36	3	1	29
	1/10 (150)	10	31,212	Dirichlet	34	5	3	30

Table 5 :

 5 Contact problem with a pre-stressed block: performance results of the FETI-C solver

	H(Ns)	H/h	N dof	preconditioner Nitr N pln Nstc N lock itr
	1/2 (6)	10	1,364	I	53	21	7
	1/2 (6)	20	5,124	I	62	43	9
	1/2 (6)	40	19,844	I	84	32	11
	1/2 (6)	60	44,164	I	105	43	17
	1/6 (54)	10	11,284	I	63	14	4
	1/6 (54)	20	44,164	I	89	58	11
	1/6 (54)	40	174,724	I	121	80	14
	1/10 (150)	10	30,804	I	68	29	8
	1/2 (6)	10	1,364	Dirichlet	59	21	7	17
	1/2 (6)	20	5,124	Dirichlet	59	38	6	18
	1/2 (6)	40	19,844	Dirichlet	69	49	12	20
	1/2 (6)	60	44,164	Dirichlet	66	55	11	20
	1/6 (54)	10	11,284	Dirichlet	47	16	7	35
	1/6 (54)	20	44,164	Dirichlet	71	41	10	39
	1/6 (54)	40	174,724	Dirichlet	73	37	13	45
	1/10 (150)	10	30,804	Dirichlet	55	40	6	38

Table 6 :

 6 Three-dimensional 6-block contact problem: performance results of the FETI-C solver

	H(Ns) H/h	N dof	preconditioner Nitr N pln Nstc N lock itr
	1/2 (6)	4	2,250	I	18	11	3
	1/2 (6)	12	39,546	I	15	18	13
	1/2 (6)	4	2,250	Dirichlet	11	6	3	6
	1/2 (6)	12	39,546	Dirichlet	8	17	7	5

Acknowledgements The first author acknowledges partial support by the Service des Relations Internationales of the ENS Cachan. The second author acknowledges the support by the Sandia National Laboratories under Contract No. BD-2435.


Lev Barinov

Sergey Ryabichko

William Bialek

Thomas Gregor

Transcription-dependent spatial organization of a gene locus

There is growing appreciation that gene function is connected to the dynamic structure of the chromosome. Here we explore the interplay between three-dimensional structure and transcriptional activity at the single cell level. We show that inactive loci are spatially more compact than active ones, and that within active loci the enhancer driving transcription is closest to the promoter. On the other hand, even this shortest distance is too long to support direct physical contact between the enhancer-promoter pair when the locus is transcriptionally active. Artificial manipulation of genomic separations between enhancers and the promoter produces changes in physical distance and transcriptional activity, recapitulating the correlation seen in wild-type embryos, but disruption of topological domain boundaries has no effect. Our results suggest a complex interdependence between transcription and the spatial organization of cis-regulatory elements.

It has long been appreciated that developing embryos provide accessible laboratories in which to study the basic mechanisms of gene regulation [1][2][3]. More recently it has become clear that embryos provide the opportunity for much more quantitative experiments, illuminating the physics of regulation and the resulting flow of information through the relevant genetic networks [4][5][6][7][8][9]. Here we exploit the striped expression patterns of the gene even-skipped (eve) in the Drosophila embryo to explore connections between the spatial organization of the chromosome and the functional activation of transcription.

In higher eukaryotes, transcription is controlled by interactions among DNA segments that are widely separated along the chromosome, specifically, cis-regulatory enhancer sequences and their target promoters [10,11]. These interactions occur within complex genetic loci containing multiple enhancers and in a spatially crowded nuclear environment [12][13][14]. It is unclear whether the three-dimensional spatial architecture of such complex loci is connected to their functional output, transcriptional activity [15]. Are there causal links between chromatin conformation and transcription?

Recent studies using imaging and chromatin capture methods provide conflicting observations regarding locus architecture. In particular, it is not clear whether physical enhancer-promoter proximity is required for transcription [16][17][18][19][20] or whether topological domain boundaries are needed to organize the spatial structure of a transcriptionally active locus [21,22]. It seems likely that we are missing something fundamental about the spatial components of transcriptional regulation.

The genomic locus of the Drosophila patterning gene eve provides an ideal system in which to address these issues [23,24]. Eve is expressed in the early embryo in seven spatially distinct stripe domains (Fig 1A -B). Each stripe is controlled by one of five enhancers that are all located within an 18 kb topologically associated domain (TAD), which is flanked by two insulator boundaries, homie and nhomie [25][26][27][28]. Further, eve expression occurs in a syncytium prior to the embryo separating into discrete cells. Thus, by observing chromosome conformations in nuclei from different stripe and interstripe domains we have an opportunity to connect spatial structure to transcriptional activity, with all other variables set by the common cytoplasm.

Qualitatively, our strategy is to mark promoter and enhancer sites throughout the eve locus with fluorescent probes and measure the distances among these sites in large numbers of nuclei across multiple stripes and interstripe domains. In detail, we use an adapted multicolor oligopaint approach [29,30]. We visualize multiple sites within the eve locus simultaneously using pools of single-stranded DNA (ssDNA) oligonucleotides, each targeting 1-2.5 kb long regions, and Concurrent visualization of DNA and nascent mRNA was performed with a catalyzed reporter deposition scheme. We achieved subdiffraction resolution with spectrally-based localization microscopy using a conventional confocal microscope to capture optical sections of the surface of blastoderm embryos (Fig 1B). We find that with careful calibration we can localize sites with ∼ 25 nm precision in the imaging plane and ∼ 56 nm along the optical axis; our respective optical resolution is ∼ 7 nm, and ∼ 55 nm. As an example, we have tagged two adjacent 1 kb segments with probes in two colors, and find these to be 111 ± 2 nm apart.

To establish the spatial scale of the eve locus, we measured distances between seven pairs of elements within the locus, across thousands of nuclei in 15-23 embryos. We can summarize the results as distributions of pairwise distances (Fig 1C ), which have median separations in the range 140 -220 nm and an overall mean of 173 ± 10 nm. The smallest median distance is between the stripe 3+7 (E3+7) and 5 (E5) enhancers, which are separated by 11 kb along the chromosome, while and the largest median distance is between the eve promoter and the 3' homie boundary, 9 kb downstream. These results support a compact globular structure for this 18 Another way to conceptualize these measurements is to define the probability of a promoter encountering an enhancer element at some distance R from itself. If we normalize this probability by the volume of a corresponding sphere with radius R, we can define an effective density or concentration, e.g. for enhancer sites as seen from the promoter (Fig 1D). As expected, if we make R very large this effective concentration decreases, but there is a plateau below R ∼ 200 nm. Inside this plateau, it is as if the promoter were in a solution of enhancer sites at a concentration in the range of 20 -100 µm -3 . This is equivalent to 30 -150 nM, comparable to the concentration of transcription factors in the surrounding solution [4][5][6]31]. Thus the promoter appears to be no more in contact with enhancers than it would be with the typical unbound transcription factor. The spatial proximity of the enhancer and promoter thus does not, by itself, provide a mechanism for transcription factor binding to control transcription.

One might worry that there is a small population of close enhancer-promoter contacts which are invisible to us because of our limited resolution. As a control we analyze the colocalization of two-color labels targeting a ∼1 kb sequence Fig 1D and see that the effective concentration is a factor of 30 to 100 times larger than for any of the enhancer-promoter pairs, plateauing just below R ∼ 20 nm. This shows that fewer than 1% of the enhancer-promoter pairs could be in close contact, even when the promoter is active, making the possibility of a proximal sub-population unlikely.

Our measurements on fixed tissue are snapshots of a dynamic, fluctuating structure. We get a more detailed picture of this structure by measuring simultaneously the positions of three sites within the eve locus, e.g. the enhancers for stripes 3+7 and stripe 5, along with the promoter. In all cases we observe significant crosscorrelations among the two pairwise distances, which suggests that the locus forms a fluctuating globule across the whole range of conditions sampled in these experiments.

When we sort nuclei based on whether the eve locus is transcriptionally active or inactive, we find that the distributions of enhancer-promoter distances are subtly different (Fig 2A). In particular, median distances are larger by ∆R = 17 ± 5 nm at active loci, independent of which enhancer-promoter pair is tagged. This ∼ 10% expansion of distances corresponds to a ∼ 30% expansion of the volume of the locus when it is transcriptionally active, qualitatively in agreement with previous work [19]. We can find inactive loci in the interstripe regions or by searching for alleles within eve-expressing nuclei that happen to have been captured in a silent state, but the degree of compaction relative to the active state is the same (Fig 2B). Thus the spatial rearrangements that we see are linked to transcriptional activity, even across alleles in individual nuclei.

The fact that different eve stripes are activated by different enhancers allows us to ask whether enhancerpromoter distances are specific to the interactions that drive transcription. We see that E5 indeed is on average significantly closer to the promoter in stripe 5 (Fig 2C). Similarly, E4+6 is closer to the promoter in stripes 4 and 6 (Fig 2D). Interestingly, enhancer-promoter distances are the same in stripes 4 and 6, which are activated by the same enhancer, and similarly in stripes 3 and 7. These differences in distance are similar to the differences seen between transcriptionally active and inactive loci and the qualitative relationship does not depend on the chosen distance metric. If we think of the individual enhancers as being active or inactive depending on the position of the nucleus, these results can be summarized by saying that enhancers and promoters are close either when both are active or when both are inactive, and farther apart when their states disagree, in qualitative agreement with recent theoretical predictions [32]. Intriguingly, as an exception to this rule, the E3+7 distance to the promoter appears to be constitutively short, which might be due to additional intra-TAD substructures between the eve promoter and the upstream insulator element nhomie, as evidenced through DNA FISH and Hi-C. Is spatial structure causally related to transcriptional activity? We should be able to change the spatial structure of the eve locus by inserting inactive segments of DNA, expanding the genomic distances. We have done this, adding 10 kb of exogenous DNA downstream of the eve promoter (Fig 3A). This increases the separation between the promoter and the two enhancers E4+6 and E5, while not affecting the separation between the promoter and E3+7 (Fig 1A). In embryos carrying this modification, we find that the fraction of transcription- We performed a second manipulation to the locus, deleting a 368 bp DNA segment containing the 3' homie insulator element. In embryos carrying this deletion the fraction of active nuclei in the different stripe domains closely matches wild-type levels (Fig 3B). In addition, no major shifts are observed between the WT and ∆Homie spatial configurations when tagging E5 and the promoter (Fig 3C). Thus, the flanking topological boundary has only minor impact on transcriptional activity and spatial organization of the active eve locus.

While none of the genomic manipulations impact the E3+7-promoter distance, the additional 10 kb genomic separation between E5 and the promoter leads to a large decrease in transcriptional activity and an increase in physical enhancer-promoter distance, with a corresponding overall decrease in effective enhancer concentration (Fig 3D). Together these results suggest that a disrup-tion of genomic organization can have a direct impact on both the spatial configuration and the activity level of intervening genes. In particular, there seems to be a strong modulation of transcriptional activity with enhancerpromoter distance at the 150 -250 nm scale (Fig 2C ), encompassing the outer edge of the globular locus structure mapped out in Fig 1D .   In summary, our measurements of enhancer-promoter distances conditioned on transcriptional state solidify a relationship between the spatial organization of the Drosophila eve locus and transcriptional activity. We find evidence for a globular locus structure with linear dimensions ∼ 200 nm, in line with recent observations of microenvironments [33], nuclear hubs [34], or liquid condensates [20]. Most strikingly, even transcriptionally active enhancer-promoter pairs are rarely closer than ∼ 150 nm, making sustained direct molecular contact highly unlikely, and calling for a novel mechanism for enhancer-to-promoter information transfer [32]. The plateau we observe in enhancer density for distances below ∼ 200 nm suggests that enhancer-promoter pairs are insensitive to distance at this length scale, and thus distance cannot be the parameter driving functional control.

Exploiting the multiplicity of the active eve domains, each driven by an independent enhancer [35], and the high-throughput nature of internally controlled side-byside experiments, we provide insights into the mechanisms governing both the spatial and transcriptional regulation of the locus at an even finer spatial scale. In agreement with previous measurements [19], we observe a general decompaction of the locus when it is transcriptionally active. In we observe enhancer-specific displacements between enhancer-promoter pairs in the domains in which the respective enhancer is driving activity. These results may help reconcile multiple conflicting observations regarding enhancer-promoter proximity as a function of activity [16,18,19,36]. It is especially interesting that enhancers driving transcription are typically the same distance from an active promoter as an inactive enhancer is from an inactive promoter.

Finally, we demonstrated that the boundaries flanking the eve TAD are not necessary for enhancer-promoter proximity or for transcriptional activity. This suggests that these boundaries may be primarily insulating, rather than mediating a necessary spatial precursor to transcription. This is supported by a clear relationship among the enhancer-promoter genomic distances, their spatial separation, and the transcriptional output. This suggests that the spatial separation of enhancer-promoter pairs is subtly but crucially coupled to the regulation of transcription, while flanking architectural mediators may not be sufficient to overcome genomic distances on the order of 10 kb within the domains they flank.

  kb locus with a spatial scale of ∼ 200 nm (Fig 1D, inset).

FIG. 1 :

 1 FIG. 1: Spatial organization of cis-regulatory elements in the Drosophila eve locus. (A) Schematic sequence of the 18 kb eve locus with 5 genomically distinct enhancer regions (E1, E2, E3+7, E4+6, and E5), flanked by insulators Homie (H) and nHomie (nH). (B) B, Sample image of a 2.5 h old Drosophila embryo with fluorescent oligonucleotides labeling the eve promoter in red, the stripe-5 enhancer E5 in yellow, and nascent eve RNA of the seven-striped pattern ( S1-S7) in green), scale bar 10 µm; inset shows close-up of small white rectangular region, scale bar 10 µm. (C) Distributions of distances between pairs of sites, eight combinations of tagged elements: Colocalization control (blue, n = 6 embryos, N = 4583 loci), E5 to eve promoter (n = 23, N = 6637), E3+7 to eve promoter (n = 23, N=6119), E4+6 to eve promoter (n = 11, N = 4222), Homie to nHomie (n = 21, N = 5251), E5 to E3+7 (n = 15, N = 3612), nHomie to eve promoter (n = 17, N = 5484), and Homie to eve promoter (yellow, n = 18, N = 5053) from blue to yellow, respectively. Inset shows box and whisker plots for all distributions with 25th, 75th, and median (EC50) values; average median value is 173 ± 10 nm (dashed line). (D) Effective density or concentration of sites as a function of the averaging volume (sphere of radius R), colors as in (C). Inset shows schematic of eve locus, with all measured points on average at a fixed distance apart, enclosing a droplet-like region (brown shade) with a spatial scale of ∼ 200 nm.

FIG. 2 :

 2 FIG. 2: Transcription-dependent and enhancer-specific spatial locus organization. (A) Cumulative pairwise enhancerpromoter distance distributions for active (green) and inactive (black) alleles. Distributions shown for distances measured between the eve promoter and three enhancers: E5 (n = 23 embryos, Non = 2471 loci, N off = 4166 loci), E4+6 (n = 11, Non = 697, N off = 3525) and E3+7 (n = 23, Non = 2571, N off = 3548). Distances are measured as differences form the median for all inactive loci. Inset shows median distances for inactive (black; average= 0±0.4 nm) and active (green; 17±5 nm) loci. (B) Pooled median distances of (eve promoter to E5 and E4+6 ) for active (green, 188 ± 9 nm, n = 34, N = 204) and inactive loci; inactive loci are split within-domain (red; 164 ± 10 nm, n = 34, N = 209) and inter-domain (gray; 164 ± 8 nm, n = 34, N = 1699). (C) Distances between eve promoter and E5, in different stripe domains (median±st. error, n = 23 embryos). Grey band shows inter-stripe distance (163 ± 4 nm, N = 1056). (D) As in (C), but for E4+6 enhancer (n=11 embryos; gray band 164 ± 11 nm, N=583).

FIG. 3 :

 3 FIG. 3: Mutual interaction between spatial locus organization and transcriptional activity. (A) Sample images of 2.5 h old embryos with fluorescently labeled eve transcripts: wild-type (left), an embryo carrying a 10 kb insertion of exogenous DNA (middle), and ∆Homie (right), with their respective genomic landscapes (top). (B) Fraction of active nuclei per stripe domain in WT (blue), ∆Homie (red), and the 10 kb insertion (green), seen in (A) (for each genetic background domain fractions are normalized to average fraction of S3 and S7). (C) Normalized fraction of active nuclei (from B) versus median distance between the promoter and the relevant domain-driving enhancer the for the three genetic constructs shown in A (same color code as in (B); black dashed line is a linear fit). (D) Effective density (for E5 and E3+7) as a function of radial distance (see Fig 1D) for embryos of all three genetic backgrounds shown in A.

We thank MS Levine, M Levo, M Nollmann, P Schedl, G Tkačik, and EF Wieschaus for helpful discussions. This work was supported in part by the US National Science Foundation, through the Center for the Physics of Biological Function (PHY-1734030) and by National Institutes of Health Grants R01GM097275, U01DA047730, and U01DK127429.