A Parret-Fréaud

C Rey

P Gosselet

F Feyel

Fast estimation of discretization error for FE problems solved by domain decomposition

Keywords: verification, error in constitutive relation, non overlapping domain decomposition, FETI, BDD

This paper presents a strategy for a posteriori error estimation for substructured problems solved by non-overlapping domain decomposition methods. We focus on global estimates of the discretization error obtained through the error in constitutive relation for linear mechanical problems. Our method allows to compute error estimate in a fully parallel way for both primal (BDD) and dual (FETI) approaches of non-overlapping domain decomposition whatever the state (converged or not) of the associated iterative solver. Results obtained on an academic problem show that the strategy we propose is efficient in the sense that correct estimation is obtained with fully parallel computations; they also indicate that the estimation of the discretization error reaches sufficient precision in very few iterations of the domain decomposition solver, which enables to consider highly effective adaptive computational strategies.

Introduction

The setting-up of robust numerical methods to solve complex systems of partial differential equations has become a key issue in applied mathematics and engineering, driven by the increasing use of numerical simulation in both research and industry. Among the latter, virtual testing has become a short term aim, with the objective to replace expensive experimental studies and validations by numerical simulations, even in order to certify large structures as planes and bridges.

Thus, one key point of the numerical methods to develop is the verification of computations which enables to warranty that the computed solution is sufficiently close to the original continuum mechanics model. This topic of numerical analysis has been the subject of many studies for the last decades. Three main classes of error estimator have been developed, based either on equilibrium residuals [START_REF] Babuska | Error estimates for adaptative finite element computation[END_REF], flux projection [START_REF] Zienkiewicz | A simple error estimate and adaptative procedure for practical engineering analysis[END_REF] or error in constitutive law [START_REF] Ladevèze | Comparaison de modèles de milieux continus[END_REF]. An overview of those various methods can be found in [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF].

Another key point of numerical methods is their ability to quickly provide solutions to large (nonlinear) systems. The most classical answer to this issue is to use domain decomposition methods in order to take advantage of the parallel hardware architecture of recent clusters and grids. In engineering, non-overlapping domain decomposition methods are mostly employed, such as the well known FETI [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] or BDD [START_REF] Mandel | Balancing domain decomposition[END_REF]. An overview of the main approaches related to non-overlapping domain decomposition can be found in [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF].

We aim to provide fully integrated adaptive strategies to compute large structural mechanics problems with certified quality. To do that, our current approach is to explore some ways of making bidirectional interactions between domain decomposition and a posteriori error estimation. Our developments are based both on the error in constitutive relation to measure the quality of our results and to forecast mesh refinement, and on a generic vision of non-overlapping domain decomposition methods which enables to do high-performance computing.

This paper focuses on the estimation of the global error in constitutive relation in order (among others) to study how it is influenced by the error in the convergence of the domain decomposition solver which is linked to the non-satisfaction of interface equations (continuity of displacements and balance of forces). To do so we propose a strategy to build, in parallel and during the iterations, displacement and stress fields which are kinematically admissible (KA) and statically admissible (SA) on the whole structure. We face two main difficulties. First, since before convergence interface fields do not possess the classical properties of discretized fields (continuity of displacements and weak equilibrium), the recovery of admissible displacements and stresses requires some preprocessing. Second, the computation of statically admissible fields being an operation which can not be conducted indepen-dently on each element (in some methods it can even be a large bandwidth operation), classical recovery methods [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Ladevèze | New advances on a posteriori error on constitutive relation in f.e. analysis[END_REF][START_REF] Parés | Subdomain-based flux-free a posteriori error estimators[END_REF][START_REF] Gallimard | A constitutive relation error estimator based on tractionfree recovery of the equilibrated stress[END_REF][START_REF] De Almeida | Recovery of equilibrium on star patches using a partition of unity technique[END_REF][START_REF] Ladevèze | A new non-intrusive technique for the construction of admissible stress fields in model verification[END_REF] would require inter-subdomain communications.

Our generic method to build continuous displacement and balanced traction fields for both primal and dual approaches of non-overlapping domain decomposition is presented through this paper. It will be shown that the properties of the preconditioners involved in domain decomposition solvers make this reconstruction costless, and that an error estimator can then be computed in a fully parallel way.

This paper is organized as follows. Section 2 recalls the general framework related to our upcoming developments, mainly the estimation of the error in constitutive equation and the use of domain decomposition method. Section 3 shows how the problem of error estimation in a substructured context can be brought back to the computation of nodal displacement and traction fields which are admissible in a discrete sense. Sections 4 and 5 describes how to obtain these fields without inter-subdomains exchanges when using classical primal (BDD) and dual (FETI) domain decomposition methods with good preconditioners. Section 6 presents numerical assessments, first to validate the parallel recovery procedure, then to prove that a good estimation can be obtained far earlier than the solver converged (in the sense of domain decomposition iterative solver). Finally, Section 7 concludes this paper.

Framework of the study

Reference mechanical problem

Let us consider the static equilibrium of a structure which occupies the open domain Ω ⊂ R d and which is submitted to given body forces f , to given traction forces g on ∂ f Ω and to given displacements u 0 on the complementary part ∂ u Ω = ∅. We assume the structure undergoes small perturbations and that the material is linear elastic, characterized by the Hooke's tensor H. Let u be the unknown displacement field, ε(u) the symmetric part of the gradient, σ the Cauchy stress tensor.

Let 1). We introduce two affine subspaces and one positive form: where tr is the trace operator.

ω ⊂ Ω be an open subset of Ω, ∂ f ω = ∂ω ∩ ∂ f Ω, ∂ u ω = ∂ω ∩ ∂ u Ω and Γ = ∂ω \ (∂ u ω ∪ ∂ f ω) (see Figure
• Subspace of kinematically admissible fields KA(ω) = u ∈ H 1 (ω) d , tr(u) = u 0 on ∂ u ω (1)
• Subspace of statically admissible fields

SA(ω) = τ ∈ L 2 (ω) d×d , τ symmetric, ∀u * ∈ KA 00 (ω) , ω τ : ε(u * )dω = ω f.u * dω + ∂ f ω g.u * dS (2)
where KA 00 (ω

) = u ∈ H 1 (ω) d , tr(u) = 0 on ∂ u ω ∪ Γ
• Measure of the non-verification of the constitutive equation [START_REF] Ladevèze | Comparaison de modèles de milieux continus[END_REF] e

CR(ω) (u, σ) = σ -H : ε (u) H -1 ,ω (3)
where

x H -1 ,ω = ω (x : H -1 : x) dω
The mechanical problem set on Ω can be formulated as:

Find (u ex , σ ex ) ∈ KA(Ω) × SA(Ω) such that e CR(Ω) (u ex , σ ex ) = 0

Finite element approximation for the global problem

Let Ω h be a tessellation of Ω to which we associate a finite dimensional subspace KA h (Ω) of KA(Ω). The classical finite element displacement approximation consists in searching

u h ∈ KA h (Ω) σ h = H : ε(u h ) Ω σ h : ε(u * h )dΩ = Ω f.u * h dΩ + ∂ f Ω g.u * h dS, ∀u * h ∈ KA 00 h (Ω) (4)
After introducing the d × N dof matrix ϕ h of shape functions which form a basis of KA h (Ω) and the vector of nodal unknowns u (of size N dof , number of degrees of freedom) so that u h = ϕ h u, the classical finite element method leads to the well-known linear system:

Ku = f ( 5
)
where K is the (symmetric positive definite) stiffness matrix of domain Ω h and f is the vector of generalized forces.

A posteriori error estimator

The finite element approximation (u h , σ h ) satisfies u h ∈ KA(Ω) and e CR(Ω) (u h , σ h ) = 0 but σ h / ∈ SA(Ω). The error in constitutive relation consists in deducing from (u h , σ h ) an admissible displacement-stress pair ( u h , σ h ) ∈ KA(Ω) × SA(Ω) in order to measure the residual on the constitutive equation (3) e CR(Ω) ( u h , σ h ) 0. Using the well-known Prager-Synge theorem it can be proved that

ε(u ex ) -ε( u h ) 2 H,Ω + σ ex -σ h 2 H -1 ,Ω = e CR(Ω) ( u h , σ h ) 2
Hence, the evaluation of the error in constitutive relation e CR(Ω) ( u h , σ h ) for any admissible pair ( u h , σ h ) provides a guaranteed upper bound of the global error

ε(u ex ) -ε( u h ) H,Ω e CR(Ω) ( u h , σ h ) (6)
KA h (Ω) being a subspace of KA(Ω), the construction of an admissible displacement field u h is straightforward since it can be taken equal to u h . On the other hand, as σ h is not statically admissible, the construction of an admissible stress field σ h ∈ SA(Ω) is a crucial point which has already been widely studied in the literature. A first solution is to use a dual formulation of the reference problem [START_REF] Beckers | 3-d error estimation and mesh adaptation using improved r.e.p. method[END_REF] to compute σ from scratch. Unfortunately building a subspace of SA(Ω) is a complex task and most people prefer to post-process a statically admissible field from Field σ h obtained by a displacement formulation. Classical methods are the element equilibration techniques [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Ladevèze | New advances on a posteriori error on constitutive relation in f.e. analysis[END_REF], which have been improved by the use of the concept of partition of unity which lead to [START_REF] Gallimard | A constitutive relation error estimator based on tractionfree recovery of the equilibrated stress[END_REF][START_REF] De Almeida | Recovery of equilibrium on star patches using a partition of unity technique[END_REF][START_REF] Ladevèze | A new non-intrusive technique for the construction of admissible stress fields in model verification[END_REF] and the flux-free method [START_REF] Parés | Subdomain-based flux-free a posteriori error estimators[END_REF]. In most cases they involve the computation of efforts on "star-patches" which are the set of elements sharing one node, for each node of the mesh. Though rather simple these computations are in great number and thus expensive.

In the following, we note by F h the algorithm which has been chosen to build an admissible stress field σh . Whatever the choice, the algorithm takes as input not only the finite element stress field σ h but also the continuous representation of the imposed forces (f, g).

σh = F h (σ h , f, g) ∈ SA(Ω)
The algorithm we have used for our applications is the one proposed in [START_REF] Ladevèze | New advances on a posteriori error on constitutive relation in f.e. analysis[END_REF] using a three degrees higher polynomial basis when solving the local problems on elements [START_REF] Babuska | Validation of a posteriori error estimators by numerical approach[END_REF].

Substructured formulation

Let us consider a decomposition of domain Ω in open subsets (Ω (s) ) 1 s N sd (N sd is the number of subdomains) so that Ω (s) ∩ Ω (s ′ ) = ∅ for s = s ′ and Ω = ∪ s Ω(s) . Let u = (u (s) ) s , we define the global assembling operator A:

u = A(u ) ⇔ u |Ω (s) = u (s) (7)
In order to reformulate the mechanical problem on the substructured configuration, we need to specify the conditions that should be satisfied at the boundary between subdomains Γ (ss ′ ) = ∂Ω (s) ∩ ∂Ω (s ′ ) . We have the fundamental properties:

A(u ) ∈ KA(Ω) ⇔ u (s) ∈ KA Ω (s) , ∀s tr(u (s) ) = tr(u (s ′ ) ) on Γ (ss ′ ) , ∀(s, s ′ ) (8) A(σ ) ∈ SA(Ω) ⇔ σ (s) ∈ SA Ω (s) , ∀s σ (s) .n (s) + σ (s ′ ) .n (s ′ ) = 0 on Γ (ss ′ ) , ∀(s, s ′ ) (9)
In other words, in order to be admissible on the whole domain Ω, not only fields need to be admissible in a local sense (independently on each Ω (s) ), but they also need to satisfy interface conditions, namely displacements continuity and tractions balance (action-reaction principle).

Finite element approximation for the substructured problem

We assume that the tessellation of Ω and the substructuring are conforming so that (i) each element only belongs to one subdomain and (ii) nodes are matching on the interfaces. Each degree of freedom is either located inside a subdomain (subscript i) or on its boundary Γ

(s) = ∪ s ′ Γ (ss ′ ) (subscript b)
where it is shared with at least one neighboring subdomain. Let λ (s) b be the vector of unknown efforts imposed on the interface of subdomain Ω (s)

h by its neighbors. The finite element problem ( 5) can be written highlighting the contributions of subdomains:

∀s, K (s) u (s) = f (s) + t (s) T λ (s) b with    s A (s) λ (s) b = 0 s A (s) u (s) b = 0 (10)
where t (s) is the discrete trace operator (u

(s) b = t (s) u (s)
) and where A (s) and A (s) are assembling operators so that A (s) enables to formulate the mechanical equilibrium of interfaces ( 9) and A (s) enables to formulate the continuity of displacements (8) (in the case of two subdomains, we have

s A (s) λ b = λ (1) b + λ (2) b = 0 and s A (s) u (s) b = u (1) b -u (2)
b = 0, see Fig. 7 for less trivial example and [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF] for more an extensive description of all operators). One fundamental property of assembling operators is their orthogonality:

s A (s) A (s) T = 0 (11)
Note that the equilibrium of subdomain Ω (s) also writes:

K (s) ii K (s) ib K (s) bi K (s) bb u (s) i u (s) b = f (s) i f (s) b + 0 (s) i λ (s) b (12)
or in an equivalent condensed form:

S (s) u (s) b = b (s) p + λ (s) b (13)
with

S (s) = K (s) bb -K (s) bi K (s) ii -1 K (s) ib b (s) = f (s) b -K (s) bi K (s) ii -1 f (s) i ( 14
)
where S (s) is the Schur complement and b (s) is the condensed right-hand side.

A posteriori error estimator in substructured context

The key point for the efficient evaluation of the error in constitutive relation in a substructured context (without overlapping) is to define admissible pairs ( u

(s) h , σ (s) h ) ∈ KA Ω (s)
×SA Ω (s) on each subdomain so that the associated assembled pair is admissible for the reference problem (A( u h ), A( σ h )) ∈ KA(Ω) × SA(Ω). Due to the absence of overlap, the additive structure of the associated error in constitutive relation leads to a fully parallel evaluation of the a posteriori error estimator:

e CR(Ω) A( u h ), A( σ h ) 2 = s e CR(Ω (s) ) u (s) h , σ (s) h 2
The application of a classical recovery strategy to compute admissible fields raises two difficulties in a substructured context. First, the star-patches can not be employed on the boundary nodes without assuming communication between subdomains. Though these exchanges would remain limited, we propose an alternate strategy to achieve full parallelism without impairing the properties of the error in constitutive relation. Second, in order to solve the substructured problem [START_REF] Parés | Subdomain-based flux-free a posteriori error estimators[END_REF] parallel strategies consist in using iterative solvers which are based on the loosening of at least one of the interface conditions which is only verified (up to a certain precision) once the solver converged. Thus recovering strategies need to be adapted so that the local fields ( u

(s) h , σ (s)
h ) satisfy the interface conditions. The aim of this section is to prove that the determination of the admissible pair (A( u h ), A( σ h )) can be brought back to the determination of nodal interface fields ( u

(s) b , λ (s)
b ) s which satisfy specific interface conditions. The construction of these nodal fields depends on the chosen domain decomposition strategy and is discussed in the following sections.

Kinematically admissible fields

In order to ensure interface Condition (8) when building u

(s) h ∈ KA Ω (s) so that A( u h ) ∈ KA(Ω), we introduce continuous interface displacement fields û(s)
bh from which we shall deduce internal displacement fields:

û(s) bh = û(s ′ ) bh , ∀(s, s ′ ) û(s) h|Γ (ss ′ ) = û(s) bh , ∀s
Since discretizations are matching on the interface, the first condition can directly be imposed on finite element nodal quantities:

u (s) b = u (s ′ ) b , ∀(s, s ′ )
In order to deduce the internal fields, one finite element problem is solved independently on each subdomain with imposed Dirichlet conditions on the interface:

u (s) i = K (s) ii -1 f (s) i -K (s) ib A (s) T u (s) b û(s) h = ϕ (s) h u (s) = ϕ h (s) i ϕ h (s) b u (s) i u (s) b u = A( u h ) ∈ KA(Ω)

Statically admissible fields

In order to ensure interface Condition (9) when building σ

(s)
h ∈ SA Ω (s) so that A( σ h ) ∈ SA(Ω), we introduce for each subdomain the continuous balanced interface traction fields F (s) bh defined on Γ (s) which satisfy:

σ (s) h .n (s) = F (s) bh on Γ (s) F (s) bh + F (s ′ ) bh = 0 on Γ (ss ′ ) Ω (s) f.ρdΩ + ∂ f Ω (s)
g (s) .ρdS +

Γ (s) F (s)
bh .ρdS = 0 ∀ρ ∈ RKA 0 (Ω (s) ) [START_REF] Babuska | Validation of a posteriori error estimators by numerical approach[END_REF] where RKA 0 (Ω (s) ) is the set of rigid body motions which are compatible with Dirichlet conditions imposed on ∂ u Ω (s) :

RKA 0 (Ω (s) ) = ρ ∈ H 1 (Ω (s) ), ρ = 0 on ∂ u Ω (s) , ε(ρ) = 0,
The last condition of ( 15) is the translation of Fredholm's alternative in order to ensure the well-posedness of the static problem on domain Ω (s) . To build these traction fields in a simple way, we associate them with the finite element nodal reaction field λ (s) b :

Γ (ss ′ ) F (s) bh • ϕ (s) j |Γ (ss ′ ) dS = λ (s) b,j (16)
where j denote a node of the interface, ϕ

A (s) λ (s) b = 0 R (s) T t (s) T λ (s) b + f (s) = 0 (17)
where R (s) is a basis of ker(K (s) ). As said earlier, the first equation corresponds to the equilibrium between subdomains. The second equation corresponds to the balance of the subdomain with respect to virtual rigid body motions (since this kind of displacement field is exactly represented in the finite element approximation, the discrete condition is equivalent to the continuous one).

As a first approach, we define F (s) bh as:

F (s) bh = ϕ h (s) |Γ (s) F (s) b (18)
where

F (s) b is the vector of nodal values of F (s)
bh and ϕ h

|Γ (s) refers to the vector of the trace on Γ (s) of finite element shape functions. Vector F (s) b is then obtained by the inversion of the (small) "mass" matrix of the interface of each subdomain. In the following, we denote by G h the previous procedure which associates a continuous balanced interface force F

F (s) bh = G h ( λ (s) b ) (19)
The traction field F

ũ(s) = K (s) + f (s) + t (s) T λ (s) b σ (s) h = F h H : ε(ϕ (s) h ũ(s) ), f (s) , g (s) , G h ( λ (s) b ) σ h = A( σ h ) ∈ SA(Ω) (20)
The use of the pseudo-inverse K (s) + is due to the potential lack of Dirichlet boundary conditions on the substructure. Displacement field ũ(s) is defined up to a rigid body motion which needs not to be determined since only the symmetric gradient of the associated displacement field is required.

It has to be noted that the fully parallel procedure G h proposed above leads to a different admissible traction field as would have been obtained using standard patch-technique [START_REF] Ladevèze | New advances on a posteriori error on constitutive relation in f.e. analysis[END_REF] (referred in the sequel as the sequential approach). Thus the use of G h implies that the parallel error estimation is different from the standard sequential one even when discrete interface conditions are satisfied. For now there are no theoretical results on the quality of the resulting fields, examples (as given in Section 6) show that sequential estimator and parallel estimator (when interface conditions have sufficiently converged, which happens very quickly) can not be distinguished.

Recovery of admissible fields in BDD

In the Balancing domain decomposition [START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Tallec | Domain-decomposition methods in computational mechanics[END_REF], a unique interface displacement unknown u b is introduced so that continuity is always insured:

u (s) b = A (s) T u b =⇒ s A (s) u (s) b = 0 (21)
Other quantities can be deduced from u b and equations (12,13):

u (s) i = K (s) ii -1 f (s) i -K (s) ib A (s) T u b λ (s) b = S (s) u b -b (s) (22)
The BDD solver consists in iteratively finding the interface displacement u b which insure global equilibrium ( s A (s) λ

(s) b = 0), 0 = s A (s) λ (s) b = s A (s) S (s) A (s) T u b - s A (s) b (s) (23)

Recovery of KA fields

In the BDD solver, kinematic interface conditions are satisfied anytime and using u

(s) b = u b enables to build u (s)
h so that u h = A u h ∈ KA(Ω). Note that all associated computations are realized during the standard resolution process so that no extra operation is required.

Recovery of SA fields

For a given interface displacement u b , we note:

λ b = s A (s) λ (s) b = s A (s) S (s) u b -b (s)
Obviously λ b is zero if and only if u b is the solution to [START_REF]Z-set user manual[END_REF]. We then define:

λ (s) b = λ (s) b -Ã(s) T λ b (24)
where ( Ã(s) ) s are scaled assembling operators so that s A (s) Ã(s) T = I.

The multiplicity scaling is a typical example of such operator Ã(s) :

Ã(s) T = A (s) T j A (j) A (j) T -1
which, in the case of two subdomains, gives Ã(s) T λ b = 1 2 λ b . In the case of heterogeneous structures, other scaled assembly operators which take the heterogeneity into account are used [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF][START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: connections and new results[END_REF][START_REF] Gosselet | A domain decomposition method for quasi incompressible formulations with discontinuous pressure field[END_REF]. b also satisfies Fredholm's alternative, we note that since R (s) is a basis of ker(K (s) ) and

It is clear that by definition, λ

K (s) ii is invertible, we have S (s) R (s) b = 0 and R (s) i = -K (s) ii -1 K (s) ib R (s)
b . The condition then writes in an equivalent condensed form:

R (s) b T λ (s) b + b (s) = 0 R (s) b T S (s) u b -b (s) + Ã(s) T λ b + b (s) = 0
Using the symmetry of S (s) (inherited from the symmetry of K (s) ) to nullify R (s) b T S (s) , the condition writes:

Ã(s) R (s) b T λ b = 0 (25)
which is exactly the balancing condition [START_REF] Mandel | Balancing domain decomposition[END_REF] of the iterative BDD solver: the residual of the BDD iterative solver [START_REF]Z-set user manual[END_REF] has to be orthogonal to all local weighted rigid body motions so that preconditioning step is well posed.

λ b = s A (s) λ (s) b (
Then we have constructed a pair of interface nodal Vectors ( u b , λ b ) which satisfy all required conditions to build admissible fields.

Note that all the involved operations are already realized during classical steps of the primal domain decomposition approach with a Neumann-Neumann preconditioner and the associated coarse problem, so that all finite element quantities (even the internal ones) are available at no cost; the only extra operations are due to the use of Algorithms G h (to compute F bh ) and F h (to compute σ h ).

Recovery of admissible fields in FETI

In the Finite Element Tearing and Interconnecting domain decomposition [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF], a unique interface effort unknown λ b is introduced so that interface equilibrium is always insured:

λ (s) b = A (s) T λ b =⇒ s A (s) λ (s) b = 0 (26)
Displacements can be deduced from λ b if it satisfies Fredholm's alternative on each substructure:

u (s) = K (s) + f (s) + t (s) T A (s) T λ b + R (s) α (s) 0 = R (s) T f (s) + t (s) T A (s) T λ b (27)
where α (s) is the unknown magnitude of rigid body motions. The FETI solver consists in iteratively finding an interface effort λ b , under the previous constraint, which insures the continuity of interface displacement: and apply algorithms G h and F h to compute σ

0 = s A (s) u (s) b = s A (s) t (s) K (s) + t (s) T A (s) T λ b + s A (s) t (s) K (s) + f (s) + s A (s) t (s) R (s) α (s) (28)
(s) h ∈ SA Ω (s) with σ h = A( σ h ) ∈ SA(Ω).

Recovery of KA fields

For a given balanced nodal interface traction λ b , we introduce, in agreement with (27), the gap of the interface displacement :

u b = s A (s) u (s) b
and we define u

(s) b = u (s) b - Ã(s) T u b
where ( Ã(s) ) s are scaled assembling operators so that s A (s) Ã(s) T = I.

Similarly to the BDD case, a typical example of such operator Ã(s) is the multiplicity scaling:

Ã(s) T = A (s) T j A (j) A (j) T -1 (29)
Note that in the case of two subdomains, we have:

Ã(s) T u b = 1 2 u b .
The connection between FETI and BDD scaling operators (even in the heterogeneous case) is given in [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF].

It is clear that by construction

s A (s) u (s) b = 0 Hence nodal interface displacement u (s)
b can be used to deduce an admissible displacement field u (s)

h so that u h = A u h ∈ KA(Ω).
Then we have constructed a pair of interface nodal Vectors ( u b , λ b ) which satisfies all required conditions to build admissible fields. Note that all the involved operations are already realized during classical steps of the dual domain decomposition approach (with built-in coarse problem) with Dirichlet's preconditioner, so that all finite element quantities (even the internal ones) are available at no cost: the quantity u b is directly available during the classical solution procedure (without computing any α (j) ) which is based on an initialization/projection algorithm [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF], and the displacement field u (s) can be defined up to an element of the kernel (a rigid body motion) since only its symmetric gradient is used during the computation of the error. The only extra operations are due to the use of algorithms G h (to compute F bh ) and F h (to compute σ h ).

Numerical assessment

In order to assess the performance of our parallel error estimator, we consider the 2D toy problem of the Γ-shape structure of Figure 2(a) which has been used in other papers like [START_REF] Chamoin | A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems[END_REF]. Plane stresses are assumed. The material behavior is isotropic, linear and elastic, with Young modulus E = 2000 MPa and Poisson's ratio ν = 0.3. The structure is clamped on its basis (whose length is denoted L) and it is submitted to traction and shear on its upperright side, while all the remaining boundaries are traction-free. Several regular meshes have been generated, constituted by triangular elements of characteristic size h = L m with m = 2, 4, 8, 16, 32. For each mesh, a sequential (mono-domain) computation is driven, followed by domain decomposition computations obtained by an automatic splitting of the mesh in an increasing number N sd of subdomains ( N sd = 2, 4, 8 when m 4 and N sd = 2, 4, 8, 16, 32 when m 8). All the computations are driven in the ZeBuLoN finite element code [START_REF]Z-set user manual[END_REF], using elements of polynomial degree p = 1. Both BDD and FETI algorithms are used to solve the substructured problems, respectively used together with Neumann-Neumann and Dirichlet preconditioners. Beside, the convergence criterion of the solver, which stands here for the interface traction gap (resp. displacement gap) in the primal (resp. dual) approach, is set to 10 -6 .

On each case, in addition to the new parallel error estimator e ddm CR , we compute the standard sequential e seq CR and the true error e h obtained using a reference field u ex computed on a very fine mesh:

e seq CR = e CR(Ω) ( u h , σ h ) e ddm CR = s e CR(Ω (s) ) u (s) h , σ (s) h 2 e h = ε (u ex -ûh ) H,Ω = ε (u ex ) 2 H,Ω -ε (û h ) 2
H,Ω

Quality of the parallel error estimator

We first study the quality of the parallel error estimator e ddm CR for computations when convergence of the domain decomposition solver is reached. As said earlier, the proposed technique does not lead to the same statically admissible field because of the special treatment of the interface traction [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: connections and new results[END_REF]. Our estimator might then be sensitive to the substructuring, we thus compare the estimations obtained with meshes of characteristic size h and decomposition into N sd subdomains. Results are given in Figure 3 and Table 1.

We observe that:

• The results obtained by FETI and BDD can not be distinguished (which is why only FETI results are given in Table 1).

• e ddm CR barely depends on the substructuring; the results are quite similar whether they are conducted on a single domain ("sequential" curve) or on N sd subdomains. Only a slight rise of the estimation can be observed when the number of interface degrees of freedom is not small compared to the number of internal degrees of freedom, which is logical since the description of interface traction fields is coarser in parallel than in sequential.

As a conclusion, the parallel error estimator e ddm CR enables to recover the same efficiency factor as the standard sequential one, while the CPU-time is divided by N sd .

Convergence of the parallel estimator along DD-solver iterations

Previous results enabled to analyse the quality of the parallel estimator when interface quantities had converged. A new feature associated to the use of an iterative solver for the domain decomposition (DD) problem is that the discretization error estimation can be conducted before DD convergence is reached, that is in presence of displacement or traction discontinuity at the interface as explained in Sections 4 and 5

We then compute the parallel error estimator e

r n = u n b Γ u 0 b Γ or r n = λ n b Γ λ 0 b Γ (30)
Classical stopping criterion for the of convergence of DD solver is this residual being below 10 -6 . Because of the similarity between the curves, the only shown cases correspond to h = L/8 and h = L/16. The curves show a rapid convergence of the parallel error estimator along iterations of the solver, so that e ddm CR can be considered as converged when FETI residual reaches an order of magnitude of 5.10 -3 or BDD residual reaches 5.10 -1 , which corresponds to at most 5 iterations whereas the solver convergence is achieved in 10 to 20 iterations.

Actually, e ddm CR is driven by both the discretization error and the convergence of the solver (interface error). The "L"-shaped curves show that the impact of residual of the DD solver is preponderant only at the first iterations (when interface fields are very poorly estimated), after e ddm CR stagnates at a value very close to e seq CR which is only associated to the discretization error. Then, it seems possible to stop the iterations of the solver far before convergence while still obtaining an accurate global estimate for the discretization error.

Figures 5 and6 show maps of the elementary contributions e estimator highlights both discretization errors (around the re-entrant angle) and lack of convergence of the solver (along the interfaces), whereas very quickly the solver (that is the interfaces) does not contribute any more to the estimator.

The various examples show that the convergence of the global estimator is due to the convergence of elementary Contributions e ddm,E CR , which means that when willing to carry out remeshing procedures, the maps obtained after few iterations of the solver are sufficient to define correct refinement instructions.

Conclusions

In this paper, we presented a new approach to handle robust model verification based on constitutive relation error in a domain decomposition context.

The method relies on the construction of fields that are kinematically and statically admissible on the whole structure. We showed that a fully parallel construction is possible even when starting from fields which do not satisfy interface conditions. The construction is a three-step procedure: first displacement and traction nodal fields are built-up so that discrete admissibility conditions are satisfied, second continuous admissible traction fields are deduced, third these fields are used as input by any classical recovery procedure. The first step is implicitly done when good preconditioners are employed within the domain decomposition methods and the second step corresponds to the inversion of small and sparse "mass" matrices.

Our first results show that not only the estimation error does not suffer from the approximation that are made at the interface in order to achieve full parallelism, but that even roughly estimated interface fields enable to obtain a good estimation of the discretization error and correct maps of elementary contributions which are required by mesh adaptation procedures. Thus not only the computational cost are divided by the number of processors but the prior obtainment of the finite element solution can be accelerated since a coarse solution is sufficient (which corresponded to 3 to 5 times less iterations in our case).

Future studies will deal with parallel mesh adaptation.

       , A (2) =        1 0 0 0 0 -1 0 0 0 0 1 0 0 -1 0 0 0 0        , A (3) =        0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 0 0 -1 0       

Figure 1 :

 1 Figure 1: Domain Ω, subdomain ω and boundaries

  the corresponding nodal component of λ (s) b . This equation then imposes that the discrete field λ (s) b and the continuous field F (s) bh develop the same virtual work in any finite element displacement field. The conditions on F (s) bh have these discrete counterparts on λ (s) b : s

  bh allows to satisfy the interface conditions associated to the static admissibility. The next step is to build internal finite element stress fields which match the associated nodal boundary field λ (s) b . This is done by solving one finite element problem on each subdomain with imposed Neumann conditions on the interface.

=

  A (s) T λ b are by construction always balanced at the interface (26) and associated to wellposed discrete Neumann problems on each substructure (27). Hence, we can directly set λ

8 , N sd = 8 )Figure 2 :

 882 Figure 2: Γ-shape structure

Figure 2 (

 2 b) shows such a decomposition for N sd = 8 and m = 8.

Figure 4 :

 4 Figure 4: Convergence of estimator vs. DD residual

Figure 5 :

 5 Figure 5: Maps of e ddm,E CR

Figure 6 :

 6 Figure 6: Maps of e ddm,E CR

Figure 7 :

 7 Figure 7: Local numberings, interface numberings, trace and assembly operators

  Convergence of Error e h and Estimators e seq CR and e ddm CR (for various N sd ) vs. element size h

				Exact error and estimate as function of mesh size
		1					primal : 2 sd
							dual : 2 sd
							primal : 4 sd
		0.5					dual : 4 sd primal : 8 sd
							dual : 8 sd
							primal : 16 sd
	log(error)	-0.5 0					dual : 16 sd primal : 32 sd dual : 32 sd sequential exact
		-1
		-1.5
		0.2	0.4	0.6	0.8	1	1.2	1.4	1.6
					log(1/h)
	Figure 3:

Table 1 :

 1 Error e h and Estimators e seq CR and e ddm CR (for various N sd ) vs. element size h the FETI (resp. BDD) residual, defined (for Iteration n) as the normalized displacement (resp. traction) gap at the interface:

	ddm CR at each iteration of the

  ddm CR at different steps of the convergence, for N sd = 8 with h e = L/2 or h e = L/8. At the first iterations it can be seen that the

		0.55	2 subdomains			0.55	2 subdomains
		0.5	4 subdomains 8 subdomains			0.5	4 subdomains 8 subdomains
		0.45	16 subdomains 32 subdomains			0.45	16 subdomains 32 subdomains
			sequential			sequential
		0.4				0.4
	e ddm CR	0.25 0.3 0.35			e ddm CR	0.25 0.3 0.35
		0.2				0.2
		0.15				0.15
	1e-05 FETI normalized displacement jump 0.0001 0.001 0.01 FETI Case: h = L/8 1e-05 0.0001 0.001 0.01 0.1 1 BDD normalized tractions jump 2 subdomains 4 subdomains 8 subdomains 16 subdomains 32 subdomains sequential parallel error estimator e 0.1 0.1 0.1 0.15 0.2 0.4 0.45 0.5 0.55 10 0.25 0.3 0.35 e ddm CR	1e-06 1e-06	e ddm	0.1 0.1 0.15 0.2 0.4 0.45 0.5 0.55 0.25 0.3 0.35	ddm,E CR 1e-05 FETI normalized displacement jump 0.0001 0.001 0.01 FETI Case: h = L/16 0.1 1e-05 to the 1e-06 1e-06 0.0001 0.001 0.01 0.1 1 10 BDD normalized tractions jump 2 subdomains 4 subdomains 8 subdomains 16 subdomains 32 subdomains sequential

CR

BDD Case

: h = L/8 BDD Case: h = L/16


P Gosselet
email: gosselet@lmt.ens-cachan.fr

D J Rixen

C Rey

A domain decomposition strategy to efficiently solve structures containing repeated patterns

Keywords: repeated patterns, quasi-cylic structures, domain decomposition methods, block Krylov solvers, FETI

This paper presents a strategy for the computation of structures with repeated patterns based on domain decomposition and block Krylov solvers. It can be seen as a special variant of the FETI method. We propose using the presence of repeated domains in the problem to compute the solution by minimizing the interface error on several directions simultaneously. The method not only drastically decreases the size of the problems to solve but also accelerates the convergence of interface problem for nearly no additional computational cost and minimizes expensive memory accesses. The numerical performances are illustrated on some thermal and elastic academic problems.

INTRODUCTION

Industrial structures often exhibit repeated patterns. Indeed the use of repetitions simplifies the design and enables to reduce costs in assembled structures because repeated components can be mass-produced; in addition symmetries enable to equilibrate inertia and to balance loadings. The existence of repetitions is naturally used when defining the CAD model of the structure, using "copy-paste" and rigid transformations. Unfortunately, taking into account the repetitions during the numerical simulation of the structure is not as an easy task. Strategies based on Fourier expansions [9] or homogenization techniques [2] are only valid in a limited range of application, most of them require either true geometrical periodicity, repetition of one single pattern or large number of repetitions; their application to non-periodic geometries, loadings and boundary conditions is not straightforward.

The crankcase shown in figure 1 and typically used in power stations of EDF (Electricité de France) is a good example of the structures we dedicate our method to: it is composed of repeated cooling winglets and disks, and a complex clamping system. The winglets and disks are repeated patterns and we will assume that the mesh of every occurrence of a given pattern has the same discretization (as is often the case in practice to simplify the mesh generation). The remainder of the structure is not a repeated pattern. The winglets and the disk have different orientations and thus the structure in its entirety cannot be considered as quasicyclic. The classical method to simulate the carter would be to obtain a complete mesh of the structure and to solve this large system (possibly using a domain decomposition algorithm with automatic substructuring software). In the present work we propose a numerical strategy based on domain decomposition methods [5, 11, 4] (see [6] for a review) and block-Krylov solvers [10, 1, 8]. The main idea behind the strategy developed here is to redefine the structure as a collection of occurrences of patterns as described in figure 2: this structure is made of two patterns, one of which occurs three times. Our method exhibits the following properties compared to classical domain decomposition methods:

• only the patterns need to be meshed (not their occurrences),

• the computational operations are realized more efficiently, • the iteration schemes takes advantage of the numerical information associated to all occurrences to converge faster. The paper is organized as follow. First, basics on domain decomposition and Krylov solvers are briefly recalled. Then the method is presented on the most simple case: the linear thermal analysis of a fully periodic structure (one repeated pattern). The method is then generalized to non-periodic structures, to elastic problems and to decompositions involving floating substructures. Eventually, conclusion and prospects are provided.

Basics on domain decomposition methods and Krylov solvers

The method we propose is a variant of the most classical non-overlapping domain decomposition methods [6]. For the sake of simplicity, we restrict our presentation to the dual domain decomposition method (a.k.a. FETI method [5, 7]) in the context of the finite element analysis of a linear elasticity problem.

Let us denote K Ω the stiffness matrix, u Ω the displacement field and f Ω the generalized forces set on domain Ω. The so called "global"(or "assembled") problem to be solved is

K Ω u Ω = f Ω (1)
We consider a partition of Ω into n non-overlapping subdomains denoted Ω (s) , and introduce interface internal forces λ b (s) imposed on subdomain Ω (s) by its neighbors. We call t (s) the trace operator which extracts boundary degrees of freedom from subdomain Ω (s) . Equilibrium of subdomains, equilibrium of interfaces and connectivity of submeshes read:

K (s) u (s) = f (s) -t (s) T λ (s) b for s = 1 . . . n n s=1 A (s) t (s) u (s) = 0 λ (s) b = A (s) T λ for s = 1 . . . n (2)
where A (s) is the signed boolean assembly operator which connects pairwise degrees of freedom on the interface. We have chosen to express the reactions on the interface of subdomains from one unique interface stress field λ insuring automatically the equilibrium of the interface reactions. The second equation corresponds to the equality of interface displacements. In the following subscript b stands for boundary and subscript i for internal, so that t (s) u (s) = u (s)

b . In case the Dirichlet boundary conditions of a subdomain are not sufficient to fix it, we call the subdomain a "floating subdomain"; the kernel of the associated stiffness matrix K (s) is denoted by N (s) , a pseudo-inverse of K (s) is denoted by K (s) + , it verifies K (s) K (s) + y = y for any y in range(K (s) ). Eliminating the displacements u (s) from the formulation (2), one finds the dual interface problem as appearing in the Finite Element Tearing and Interconnecting (FETI) method [5]:

  F G G T 0     λ g   =   d e   (3)
where we introduced the following definitions: s) . . .

F = s A (s) t (s) K (s) + t (s) T A (s) T G = . . . A (s) t (s) N (s) . . . d = s A (s) t (s) K (s) + f (s) e T = . . . f (s) T N (

We will also denote local Schur complement by S

(s) = K (s) bb -K (s) bi K (s) ii -1 K (s)
ib . It can easily be proven that S (s) + = t (s) K (s) + t (s) T and thus F = s A (s) S (s) + A (s) T . In order to take advantage of the additive structure of matrix F , a Krylov iterative solver is chosen to solve the dual interface problem (3). Our strategy is independent of the Krylov solver. For simplicity reasons we choose, as in the basic FETI methods, to use the conjugate gradient. We hence suppose that F is symmetric positive definite (properties which are directly inherited from matrix K Ω ).

One key point of coupling domain decomposition methods and Krylov iterative solvers is the preconditioner F -1 . Various versions have been proposed for the FETI method: the so called "Dirichlet", "lumped", "superlumped" and "identity" preconditioners (in increasing order of efficiency and computational burden):

F -1 D = s ∆ (s) A (s) S (s) A (s) T ∆ (s) F -1 L = s ∆ (s) A (s) K (s) bb A (s) T ∆ (s) (4) F -1 SL = s ∆ (s) A (s) diag(K (s) bb )A (s) T ∆ (s) F -1 I = I (5)
where ∆ (s) is a diagonal scaling operator ( s ∆ (s) = I).

Another keypoint in solving the dual interface problem is the coarse problem imposing at every iteration the self-equilibrium constraint of the interface forces, the second set of equations in (3). The coarse problem is solved using a projector Q and initialization λ 0 , defined by

Q = I -F -1 G G T F -1 G -1 G T (6)
λ 0 = F -1 G G T F -1 G -1 e (7)
where F -1 can be any of the preconditioners introduced above.

The standard FETI solver is summarized in the left column of table I. Conjugate gradient consists in building a F -orthogonal basis (w i ) 0 i m of Krylov subspace

K m z 0 , Q F -1 Q T F = span z 0 , . . . , Q F -1 Q T F m-1
z 0 and finding approximation λ m ∈ λ 0 + K m z 0 , Q F -1 Q T F so that at each iteration the error λ -λ m F is minimized (which is equivalent to making the residual r m = (d -F λ m ) orthogonal to K m ). Because of the good conjugation properties of basis (w i ) the optimization is decoupled and only one scalar coefficient α m is sought for at each iteration so that λ -λ m-1 -α m w m F is minimized.

Modification of the FETI method for periodic structures

For pedagogical reasons the method is first presented in the simplest case, namely a structure made of a single pattern repeated in a periodic manner; more complex situations are considered in the section following this one.

Principles of the method

In this toy problem, we consider a periodic structure. In order to introduce concepts without excessive technical difficulties we first focus on scalar problems so that no local frame is

2 Compute λ 0 = F -1 G(G T F -1 G) -1 e + Qλ 00 3 Compute r 0 = Q T (d -F λ 0 ) 4 z 0 = Q F -1 r 0 set w 0 = z 0 z 0 = Q F -1 r 0 set w 0 = z 0 obtain w 0 through permutations 5 for j = 0, . . . , m 5.1 p j = Q T F w j p j = Q T F w j 5.2 α j = (w T j r j )/(w T j p j ) α j = (w j T p j ) -1 (w j T r j ) 5.3 λ j+1 = λ j + α j w j λ j+1 = λ j + w j α j 5.4 r j+1 = r j -α j p j r j+1 = r j -p j α j 5.5 z j+1 = Q F -1 r j+1 z j+1 = Q F -1 r j+1 obtain z j+1 through permutations 5.6 β j = -(p T j z j+1 )/(w T j p j ) β j = -(w j T p j ) -1 (p T j z j+1 ) 5.7 w j+1 = z j+1 + β j w j w j+1 = z j+1 + w j β j required.
All substructures thus have the same local Schur complement S (s) = S () . To illustrate our discussion, one can think of the toy problem of determining the temperature field on a ring whose internal face is submitted to a given temperature (so that no zero energy mode exist in our problem), the external face being submitted to non-periodic heating. To fix the ideas, we assume the pattern to be a third of the structure, see figure 3.

∂ n Ω Ω (1)
Ω ( 2)

Ω (3) ∂ (1,2) ∂ (2,1) ∂ (2,3) ∂ (3,2) ∂ (3,1) ∂ (1,3) ∂ d Ω Figure 3: "Donut" toy problem
The main assumption underlying the method is that the numerical information generated at the interface of one occurrence is pertinent for all identical occurrences in the system. So, given a search direction w i and instead of "simply" finding optimal length α, a block of search direction w i is generated by permutation of the information. Referring to figure 3 for the definition of interfaces,

w i =     w (1,2) w (2,3) w (3,1)     -→ w i =     w (1,2) w (2,3) w (3,1)
w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF]

w (1,2) w (2,3)     (8)
The undertilde notation is used to denote the block of permutations. As can be seen three different directions of descent have been constructed by simple permutation of the interface partitions of the preconditioned residual w i . Then optimization can be realized simultaneously on these three search directions and an optimal linear combination of these search directions (represented by vector α) can be deduced so that x i+1 = x i + w i α i . Obviously the search directions w i+1 for the next iteration have to be made orthogonal to previous. Depending on the loading and the structure, the method may converge much faster than classical methods. In any case the convergence can not be deteriorated because the initial information is still present in the multivector basis. The algorithm is presented in the right column of table I. As can be seen it is quite similar to the classical conjugate gradient since only a few more operations are involved. The following section discusses how an efficient implementation can be obtained from the pattern-based framework.

It is to be noted that the algorithm presented in the right column of table I applies independently of the way the multivector is generated and thus, as shown in section 4, it can be applied to non-fully-periodic structures. In the case of truly-periodic structures (namely where the interface trace operators A (s) are cyclic) computational costs can be minimized when performing products and when orthogonalizing the new research direction with respect to previous ones: indeed, because F and F -1 themselves exhibit periodicity properties, the multivectors p j and z j+1 can also be deduced from one single vector. For instance, it is sufficient to compute (line 5.1) p j = Q T F w j then to deduce p j and w j ; it is also sufficient to compute (line 5.5) z j+1 = F -1 r j+1 then to make z j+1 orthogonal to previous multivector of search directions w j ((line 5.6) β j = -(w j T p j ) -1 (p T j z j+1 ) and (line 5.7) w j+1 = z j+1 + w j β j ) and to finally deduce w j+1 through permutations.

In the following section we discuss how the computations in the multivector FETI algorithm can be organized efficiently for periodic structures.

Efficient implementation of the blocked algorithm

It is first to be noticed that our method introduces nearly no additional computations compared to the basic method. The main differences are that dot-products are replaced by matrixproducts p j T w j and that the resulting square matrix needs to be factorized. Moreover, giving the pattern a predominant role leads to the rewriting of all the operations in a computationally more efficient way: the main idea is that operations will be realized simultaneously on blocks of vectors (instead of on single vectors). The principle is that the complete interface can be observed from the point of view of the pattern with each occurrence of the pattern describing one part of the complete interface. Classically the interface vector w is stored as subdomain contributions w (s) themselves split in interface contributions w (s,i) . Here however the complete interface is stored as a matrix, named W , of local contributions (refer to figure 3):

w =     w (1,2)
w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF]

w (3,1)     w (1) =   w (1,2) w (1,3)   W =   w (1,2) w (2,3) w (3,1)
w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF]

w (3,2)   (9)
where we define w (i,j) = -w (j,i) . Clearly W stores the same information as w. The redundancy inside W , i.e. storing w (i,j) and w (j,i) corresponds simply to storing the information on different subdomains as is usually done anyway in parallel implementations when organizing the communications. In the following capital letters refer to quantities stored in this specific "block" format. Let S () be the Schur complement of the pattern which, in the case of thermal problems, is identical for all its occurrences. Computing a product F w = s A (s) S () + A (s) T w can be realized by computing S () + W and using adapted assembly operations. The interest of such an implementation is that the solutions for all the right hand sides (i.e. all columns of W ) are computed simultaneously, which enables to minimize the cost of memory transfers which represent a significant contribution to the overall cost in sparse matrix operations.

Let us now analyse the extra cost associated to the minimization of the residue on a set of research directions. The issue of applying F and F -1 on multivectors was discussed above. One must however note that to put the result of the operation in the multivector format as described for W in (9) a slightly modified implementation of the assembly procedure is required. The matrix operations such as p T w can directly be deduced from the computation of P T W (using again a modified assembly procedure). In this way this operation is hardly more expensive than simply computing the vector dot product p T w. Finally the only significant extra-cost is caused by the factorization of matrix p T w ; since, in practical problems, the number of occurrences of a pattern is not extremely high the dimension of this matrix is expected to remain small. The issue related to the possibility of this minimization matrix to become singular is discussed in following section.

Numerical issues

As said earlier the proposed algorithm optimizes many computations associated to domain decomposition methods, and only introduces one additional operation: the factorization of matrix (p T w ) required during the computation of optimal descent coefficients α j . It might occur that for some iteration this matrix became non-invertible or at least very poorly conditioned. The bad conditioning of matrix (p T w ) can be traced back the non full-columnrankness of the block of search directions w which is itself associated to a redundant information between (at least) two interfaces. In that case the use of a pseudo inverse is sufficient to avoid breakdowns.

It has been observed that this algorithm, due to the fact that it handles several directions simultaneously, is more sensitive to roundoff errors than classical conjugate gradient schemes. For that reason full re-orthogonalization of the directions of descent is often used in order to make the iteration procedure more robust. The F -normalization of search directions is then a simplifying feature; after the computation of η = (p T w ), one substitutes w ← w η -1 2 and p ← p η -1 2 such that (p T w ) = I. As a result no system of equation needs to be solved during the full-reorthogonalization and residue minimization steps. Obviously the normalization procedure described above is the right place to handle redundancy of directions of descent in the multivector. As can be seen, using several directions of descent per iteration can reduce the number of iteration by about 35% (respectively 55%) for the five-repetition example (respectively for the nine-repetition test), reminding that the cost per iteration is nearly identical. As can be observed in almost every cases, the norms of the residual during the first iterations of the block and non-block methods are almost coincident. Indeed, the first iterations of conjugate gradient explore the higher part of the active spectrum of the problem which corresponds to large wavelength phenomena which are almost identical in any column of the multivector. Once these contributions are found, the enrichment of the search space offered by multivectors becomes effective and leads to significant improvement of the convergence.

Five repetitions Nine repetitions

Donut tests

In order to evaluate the performance results on more realistic cases the conductivity matrix of a piece of donut is obtained from FreeFem++ (www.freefem.org) and imported in Scilab. The pattern in these examples has about 2 000 degrees of freedom (dof) of which about 100 are on its interface. Here again the loading is random and non-periodic. The Dirichlet preconditioner is used. Figure 5 presents the number of iterations of the conjugate gradient to reach a given precision for a classical dual domain decomposition method and for the block dual domain decomposition method. The algorithms are assessed on the five and nine-part donuts. The decrease of the number of iterations is less spectacular than in previous cases, although with 9 repetitions the block algorithm converges in about 50% less iterations.

Pattern (Five-part donut)

Five-part donut Nine-part donut It is also interesting to estimate the computational costs of the classical FETI solver and the block variant for systems with repeated patterns (and applied to periodic structure in this section). Here the algorithm is programmed in scilab (an interpreted programming language) and cpu time are only indicative. Nevertheless it is informative to compare the CPU time for three different variants of the FETI solver:

• A basic ("classical") FETI method where the fact that the structure is made of repeated patterns is not at all accounted for. • The same classical FETI algorithm but where the periodicity is utilized to compute efficiently F w = s A (s) S () + A (s) T w. Since every subdomain has the same local operator all subdomain contributions in this interface problem are computed simultaneously, namely applying the forward and backward substitution at once for several right hand sides. Clearly this algorithm will exhibit the same convergence history but the cost related to the factorization of S () is paid only once, and the cost per iteration is significantly reduced. This version of the FETI method will be called the classical FETI-mrhs. • Finally we apply the multivector FETI algorithm as proposed in this paper where several directions of descent are considered per iteration whereas the cost per iteration is similar to the FETI method denoted above as the classical FETI-mrhs.

The cpu times are listed in Table II. Obviously the performance of the classical FETI is significantly penalized by the fact that no use is made of the property that local operators are identical; nevertheless if implemented on a multiprocessor computer and assuming perfect parallelism (i.e. dividing the cpu time for the classical FETI by the number of domains) the CPU time would be very close to the one of the classical FETI-mrhs. Hence in a parallel computing setting using the fact that local operators are identical as in the classical FETImrhs approach is not very advantageous. When applying the multivector FETI however faster computing will be achieved thanks to the fact that the entire information generated at every iteration is exploited in the multiple search directions, thereby leading to a faster convergence with a similar cost per iteration. It amounts to a gain of 50% in computing time in our test for the 9-part donut case. When a structure can not be decomposed into several occurrences of a single pattern, it cannot be treated with the procedure described in the previous session dealing with periodic structures. It is to be noticed that boundary conditions and/or connection conditions between subdomains can also alter the periodicity. Figure 6 presents different cases where the strategy developed in the previous section needs to be slightly modified to still be able to make efficient use of the presence of patterns.

Ω (1) Ω ( 2)

Ω (3)
Ω [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] (a) Donut with one stand

Ω (1)
Ω ( 2)

Ω (3)
Ω (4) Ω [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] (b) Donut with two stands

Ω (1)
Ω ( 2) There are thus two patterns, one of which being instantiated three times. The topology of interfaces is modified and the system cannot be seen as periodic anymore. Nevertheless one can still build a multivector that makes use of the information computed on the interfaces between repeating subdomains. The simplest multivector one can consider is:

Ω (3) ∂ n Ω ∂ d Ω (c) Non periodic Dirichlet bc's
w =        w (1,2)
w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF])

       w =        w (1,2) w (2,3) w (3,1)
w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF])

       (10)
In the present work the method is assessed on academical matrices (randomly generated as in the examples of section 3.4.1). Each interface is of dimension 20 and either 5 or 9 repetitions are considered. In figure 7 we plot the convergence history of the classical and the multivector FETI. It is observed that the exploitation of the repetition leads to high efficiency, the number of iterations being divided by a factor greater than 2. It is also remarkable that the multivector method is independent on the number of repetitions. It thus seems that the proposed multivector FETI approach leads to convergence efficiency that is scalable with the number of repetitions. A possible interpretation of this result is that the exploitation of the repetition has the same algorithmic effect as the use of a coarse grid problem (typically associated to rigid body motions of floating substructures): at each iteration the residual is made orthogonal to a subspace augmented by the permutations (that is to say a space that is larger than the single vector of classical conjugate gradient). The permutations enable to transmit instantaneously a global information on the whole structure which is precisely the aim of coarse grid problems.

Donut with two stands

In the case of two identical stands (figure 6.b) new partial periodicity appears. It is possible to consider the following enriched multivector:

w =           w (1,2)
w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF]

          w =          
w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF][START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Dostál | Total feti -an easier implementable variant of the feti method for numerical solution of elliptic pde[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Bergman | Les méthodes de l'homogénéisation -Théorie et applications en physique[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] w [START_REF] Aliaga | A lanczos-type method for multiple starting vectors[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]

          (11)
Of course one danger of this philosophy is an exponential increase of components in the multivector as the number of patterns and repeated interfaces grows. Proper strategies to choose the permutations are yet to be developed in that case.

Dirichlet bc's

Let us now consider the case of one piece of the donut having part of its outer frontier submitted to Dirichlet conditions (imposed temperature, figure 6.c); in that case the Schur complement of the subdomain should be modified. A first option would be to consider that the structure is made out of two patterns (the two free pieces versus the clamped one). Another option, following the Total FETI method proposed in [3], is to use Lagrange multipliers to impose the boundary conditions and thus the Schur complements are always associated to free substructures so that all subdomain can be represented by the same pattern.

The preconditioning problem

In any of the last two cases we showed that non-periodic interface conditions (extended to the boundary of the subdomain) could be integrated into our acceleration framework. One key point which gives precedence to the dual approach is that the dual Schur complement involves the resolution of Neumann problems associated to matrix K (s) + which is independent on the location of the imposed fluxes, so that a block treatment of the F w product is always possible.

The problem is different when dealing with the Dirichlet preconditioner since it is strongly dependent on the location of the Dirichlet conditions (which modifies the splitting between interior and boundary degrees of freedom). Thus the theoretically optimal preconditioner is made out of contributions which are potentially different even for subdomains associated to the same pattern. As an example, consider again the structure with one stand (figure 6.a): the Dirichlet operator of subdomain 3 is different from the one of subdomains 1 and 2 because it should integrate the action of the stand. Thus a block treatment of the optimal preconditioner is complex to manage; anyhow lumped strategies are still available for block treatment and should lead to good performance results.

Linear elasticity problems

For now we have considered thermal problems because the research of a scalar unknown leads to simpler writing. To use previous strategy for elastic problems one has to take into account the rotations which are necessary to transform the pattern to one of its occurrences. These rotations modify the interpretation of the degrees of freedom of the interface nodes. Nevertheless this additional local rotation step does not modify the strategy proposed here (see for instance [9] for a detailed formulation of periodic elastic problems in the dual setting).

More precisely two modifications are required:

• The assembly operator needs to incorporate rotation from local to global frame. Since Schur complements are defined in the pattern's frame, not only the assembly operator selects degrees of freedom associated to the interfaces of specific occurence but it expresses them in the global frame: A (s) T goes from the structure's frame to the pattern's frame and A (s) sends back data defined on the pattern to data defined on the structure. • The permutations which enable to obtain multivectors (w to w ) also need to take into account the rotation. Data defined on interface ∂ 1,2 in the structure's frame are not directly pertinent for interface ∂ 2,3 , they need to be transformed to be defined in the same frame.

These operations are easy to implement and do barely incur additional computational costs.

As an example we consider the "fully-periodic" donut as a bidimensional elasticity problem (plane strain). The internal ring is clamped and loading is the result of a random process. The pattern has about 2000 degrees of freedom (1000 nodes) of which about 100 are on its interface. The acceleration due to the use of multivectors is very similar to the one observed for thermal problems, performance are summed up in table III. Here again the use of multivectors seems to provide numerical scalability, CPU gain settles around 25% (which could be improved if some redundant operations where suppressed in the yet non-optimal implementation of the multivector approach). Table III: Performance of the multivector approach for elastic donut problems

Floating substructures

Part of the treatment of rigid body motions is simplified in our framework because it can be realized at the scale of the pattern (and not at the level of its occurrences). For instance in the fully-periodic elastic donut problem, assuming the pattern is only fixed by a hinge at the central node of its internal ring, each occurence has one rigid body motion (rotation around the hinge) and the structure has enough Dirichlet boundary conditions (see figure 8).

G =     G (1,2) G (2,1) 0 0 G (2,3) G (3,2) G (1,3) 0 G (3,1)     =     O (1) a N a -O (2) b N b 0 0 O (2) a N a -O (3) b N b -O (1) b N b 0 O (3) a N a     (12)
where matrix O

(1) a

is the rotation matrix which transforms interface a of the pattern into interface a of occurence (1).

The computation of the coarse grid operator (G T F -1 G) sums up to the evaluation of the reaction of the pattern to traces of rigid body modes on its interfaces. So one just has to introduce matrix G block defined in the frame of the pattern:

G block =   N a N a 0 N b 0 N b   (13)
and deduce (G T F -1 G) from the computation of G T block (∆ S() ∆ T )G block (∆ is the block adaptation of the scaling matrix ∆ () , S() is the chosen approximation of the primal Schur complement of the pattern). Hence the coarse grid problem can be built by computing the following responses of the pattern: its response when submitted to a rigid body mode on its entire interface and its response to the trace of the rigid modes on each of its interface segment successively.

Results presented in table IV (obtained with S() = I) show that rigid body motions have no influences on the performance of the multivector technique. Such a result is in good agreement with previous analysis: rigid body motions are associated to a large wavelength information whereas permutations deal with small wavelength phenomena so that the two effects reinforce one another.

Conclusion and outlook

This paper presents the first developments around a new approach to compute structures made out of repeated patterns. The solution algorithm is based on a nonoverlapping domain decomposition method (in order to extract the patterns from the mesh) and on block Krylov solvers (in order to distribute time consuming operations and optimize the solution procedure on a larger subspace). The method always converges at least as fast as classical method with a cost per iteration only slightly higher. We have illustrated the efficiency of the proposed strategy on simple examples.

Few drawbacks need to be investigated yet, mostly the problem of the generation of redundant information (which leads to handling too much data and to numerical difficulties).

The extension of the method to more realistic problems is a rather complex task which at first implies an object oriented description of the mesh. The entire system is then based on the patterns and connections between their different instantiations. This requires some adaptation of usual codes and probably a better cooperation with CAD modelling tools (where the information related to the existence of the patterns and their occurences is directly available). Moreover, it is important to note that the same way rotations are used to transport information in the correct frames, non-boolean assembly operators (as obtained with mortar techniques) could be employed to deal with non-matching meshes and transport information from one discretization to another, which would avoid to impose too complex contraints on the meshing process (patterns could be meshed independently from one another, and faces of the same pattern would not need to have the same discretization).

Nevertheless we believe that handling such a description would lead to high computational benefits: only patterns would be meshed (saving memory) and would communicate with one another larger amount of data at one time (saving message passing time), operations would always be realized on blocks and minimization would be realized on larger subspaces.

Figure 1 :

 1 Figure 1: EDF carter (courtesy of O. Boiteau)

Figure 2 :

 2 Figure 2: Schematic view of structures with occurrences of a pattern

1 .

 1 Academical tests Here the Schur complement of the pattern is represented by a random dense matrix made symmetric definite positive. The Dirichlet's preconditioner is used. The pattern is repeated either five or nine times. The loading is random and non-periodic. Figure4presents the number of iterations of the conjugate gradient to reach a given precision for a classical dual domain decomposition method and for the block (multivector) version (implemented in Scilab, www.scilab.org).

Figure 4 :

 4 Figure 4: Comparison between new (dashed) and classical (solid) algorithms

Figure 5 :

 5 Figure 5: Comparison between new (dashed) and classical (solid) algorithms

Figure 6 : Various extensions 4 . 1 . 1 .

 6411 Figure 6: Various extensions

Figure 7 :

 7 Figure 7: Comparison between new (dashed) and classical (solid) algorithms

Figure 8 

 8 Figure 8: 3-piece structure with hinges

Table I :

 I Dual Schur complement with conjugate gradient

	Classical version	Multivector version
	1	Initialize with arbitrary λ 00

Copyright c 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1-6 Prepared using nmeauth.cls

Acknowledgements P. Gosselet wishes to thank Erasmus European Exchange Program for partly funding his visit to Delft University of Technology where part of this work was realized.