# Comparison of multi-level approaches in domain decomposition for structural analysis

Pierre Ladevèze, David Dureisseix

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Ladevèze, David Dureisseix. Comparison of multi-level approaches in domain decomposition for structural analysis. 2nd Euro-Conference on Parallel and Distributed Computing for Computational Mechanics - EURO-CM-PAR-98, Apr 1998, Sintra, Portugal. pp.55-63. hal-00322448

HAL Id: hal-00322448

https://hal.science/hal-00322448v1

Submitted on 24 Jul 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Comparison of multi-level approaches in domain decomposition for structural analysis

P. Ladevèze and D. Dureisseix

## Abstract

Domain decomposition methods are well suited to multiprocessor computers with either shared or distributed memory, which are now the most powerful computers known, [12].

When increasing the number of subdomains, the maximum efficiency expected from domain decomposition methods is not attained. Use of a global mechanism to propagate information among all of the substructures can overcome this drawback [1]. Such a mechanism has now been implemented on several domain decomposition-like algorithms, such as the FETI method [5] or the “Balancing Domain Decomposition method” [11].

The approach used in this work is based upon the LArge Time INcrement (LATIN) method, coupled with a substructuring technique, in order to solve the problem concurrently [7]. The initial structure is decomposed into substructures and interfaces. Since each possesses its own behaviour and equations, the unknowns in the problem are both the displacement and the efforts on the interfaces. The resulting approach then is a “mixed” domain decomposition approach. In initial extension to the LATIN approach, which takes into account the two scales arising from the substructuring, has been applied in order to improve performance [4]. A more complete extension pertains to the homogenisation techniques between a macro level and a micro level, which are also derived from the substructuring. This approach is described for the case of linear elasticity and perfect interfaces. The macro level is chosen to be the “homogenised” structure (in order to represent the large-scale effects), to propagate the information globally. The choice of this macro-level problem then allows us to express the operators that transfer information between the micro and macro levels throughout the iterations of the algorithm.

The displacement within each substructure, as well as the displacement field on the interfaces, are additively split according to the definition of the two levels. The stress field and effort on the interfaces are then split by duality when separating the energy on the two levels.

As a first step, we focus herein on linear elasticity and perfect interfaces; however, the resulting multi-level approach naturally takes place within the LATIN framework, which has formerly been designed for non-linear evolution problems. Last, a comparison of the mixed domain decomposition methods is proposed.

## 1 Introduction

Domain decomposition methods have proved to be useful iterative algorithms for performing structural analyse on parallel

computers. In order to improve the efficiency of these methods, a global mechanism to exchange the information globally among all of the subdomains at each iteration is required, [1].

Using the LArge Time INcrement method, a mixed domain decomposition method can be built, [9, 2], and an initial extension, which takes into account two scales arising from the substructuring, has been applied in order to improve performance. Its principles have been presented in [3] and its feasibility has been shown for linear elasticity and perfect interfaces in [4]. Its performance is similar to that obtained from various versions of the FETI method.

The resulting algorithm is numerically scalable, i.e. the number of iterations needed to reach convergence is independent of the number of substructures. At each iteration, the algorithm leads to solving a large-scale problem on the whole structure in order to propagate information.

The approach proposed herein is related to the homogenisation techniques between a macro level and a micro level, which have also been derived from the substructuring.

In Section 2, we specify the notations for the reference problem and its substructuring. Then, two scales are defined for both the kinematic and static quantities in Section 3. The LATIN method is used to build an iterative algorithm in Section 4, and Section 5 presents the first results for a 2D elasticity analysis. The parallel computation is performed with an initial version of such a micro-macro 2-level approach.

## 2 Problem to be solved

The reference problem is related to the quasi-static behaviour of a structure denoted by $\Omega$, over a time interval $[0, T]$, for small perturbations and isothermic evolutions. Loadings are a prescribed displacement $U_d$ on one part of the boundary $\partial_1\Omega$, a prescribed traction force $F_d$ on a complementary part of the boundary $\partial_2\Omega$, and an eventual body force $f_d$ onto $\Omega$, see Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Reference problem</figcaption></figure>
In order to simplify the presentation, the reference problem is written for the case of linear elasticity (the non-linear case is treated with a full LArge Time INcrement approach [8]). Time is no longer taken into account because only the final configuration is of interest herein.

The displacement field at each point M of Ω is U(M); the associated space is U(U₀ will denote the related virtual field space). ε is the strain field and the current state of the structure is given by the stress field σ.

The reference problem is then to find s = (U; σ) which satisfies:

* kinematic equations:

$$ \epsilon = \epsilon(U), \quad U|_{\partial_1 \Omega} = U_d, \quad U \in U \quad (1) $$

* equilibrium equations:

$$ \forall \underline{U}^* \in \mathcal{U}_0, \quad \int_{\Omega} \mathrm{Tr}[\sigma \epsilon(\underline{U}^*)] d\Omega = \\ = \int_{\Omega} f_d \cdot \underline{U}^* d\Omega + \int_{\partial_2 \Omega} \underline{F}_d \cdot \underline{U}^* dS \quad (2) $$

* the constitutive relation (K is Hooke's tensor):

$$ \sigma = \mathbf{K}\epsilon \quad (3) $$

In order to obtain an accurate solution to the previous problem when several structural effects occur — such as high-stress gradient areas — or when a precise model of the material's behaviour is needed, one would expect to use a large number of degrees of freedom in a finite element approach. This problem will be called the micro-scale problem. Obtaining a direct solution to the fine-scale problem can be very costly. One way to solve it efficiently is to use a substructuring approach, thereby allowing us to perform the necessary treatment on multi-processor computers [12].

The originality of the substructuring approach used herein is to decompose the former structure into substructures and interfaces, with each of them being a mechanical entity in its own right: each possesses its own behaviour and equation [7]. Such substructuring leads to introducing additional quantities; the state of the liaison Γ<sup>EE'</sup> between substructures Ω<sup>E</sup> and Ω<sup>E'</sup> is defined by efforts and displacements on the interface, (F<sup>E</sup>; W<sup>E</sup>) and (F<sup>E'</sup>; W<sup>E'</sup>). They represent the action of Γ<sup>EE'</sup> on Ω<sup>E</sup> and Ω<sup>E'</sup>, respectively, as shown in Figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2: Substructure and interface</figcaption></figure>

As the size of the reference problem increases, domain decomposition methods become more efficient than direct methods. However, efficiency is not as high as expected when

the number of substructures increases; in particular, numerical scalability is not attained. An additional effort is required in order to improve the performance of such an algorithm. In classical terms, a global problem has to be used for the whole structure [1]. For instance, the FETI method [5], which is a dual Schur complement problem solved with a projected conjugate gradient, must balance the efforts field on its boundary globally for each subdomain. At each iteration, it has to solve a global problem whose unknowns are the rigid body movements of the subdomains. When using a primal Schur complement method [10], the Balancing Domain Decomposition method [11] makes use of a similar feature during the preconditioning step, in which a Neumann problem is solved for each subdomain at each iteration.

In the case of the LATIN method, an initial multi-level extension to the algorithm has been proposed in [3]. It employs a coarse discretisation of the whole problem at the large-scale in order to propagate the information globally among all of the substructures. We are proposing herein another large-scale problem generated from a homogenisation-like approach, [8], which is expected to be more efficient.

# 3 Using two levels

The solution is considered to possess two parts, related to the micro scale m and to the macro scale M. Such an approach uses homogenisation techniques similar to those which can be developed with multigrid approaches [6]. The state of the structure is then described with two parts, each defined on one scale.

In order to simplify the notations, further substructure-related quantities will no longer be denoted by a superscript E.

## 3.1 Kinematic splitting

According to the previously considered splitting, the displacement field is

$$ \underline{U} = \underline{U}^m + \underline{U}^M $$

U<sup>M</sup> is regular over the whole structure, except on any eventual discontinuity surfaces (unilateral contact, for instance): U<sup>M</sup> ∈ U<sup>M</sup>. This constitute the large variation length displacement. U<sup>m</sup> is not required to be continuous throughout the interfaces between substructures (though the converged solution will be): U<sup>m</sup> ∈ U<sub>E</sub><sup>m</sup>. This is the additive micro correction displacement.

Let us consider a classical displacement-oriented discretisation of the problem, and let us focus on the case where both the macro-scale and the micro-scale unknowns are searched in embedded discretised spaces, see Figure 3. The macro-scale unknowns are hierarchically projected and the splitting on macro and additive micro spaces is performed via a hierarchical h-basis splitting, [15, 14].
<figure><img src="image_5.png" /><figcaption>Figure 3: Hierarchical discretisation of the two levels</figcaption></figure>

The strain field is then:

$$\varepsilon(\underline{U}) = \varepsilon(\underline{U}^m) + \varepsilon(\underline{U}^M) = \varepsilon^m + \varepsilon^M$$

From a kinematic point of view, the displacement on the interfaces is also split on two scales:

$$\underline{W} = \underline{W}^m + \underline{W}^M$$

with the corresponding spaces: $\underline{W}^m \in W_E^m$, and $\underline{W}^M \in W_E^M$.

## 3.2 Dual static splitting

Both stresses and efforts can be derived from previous representations, with an energy splitting, which is a key point in the proposed approach:

$$\begin{aligned} (\underline{F}, \underline{W}) &= (\underline{F}, \underline{W}^m) + (\underline{F}, \underline{W}^M) = \\ &= (\underline{F}^m, \underline{W}^m) + (\underline{F}^M, \underline{W}^M) \end{aligned} \quad (4)$$

where $(\cdot, \cdot)$ is the symmetric form:

$$(\underline{F}, \underline{W}) = \sum_E \int_{\partial\Omega^E} \underline{F}^E \cdot \underline{W}^E d\Omega$$

and:

$$\begin{aligned} \int_{\Omega} \mathrm{Tr}[\sigma\varepsilon]d\Omega &= \int_{\Omega} (\mathrm{Tr}[\sigma\varepsilon^m] + \mathrm{Tr}[\sigma\varepsilon^M])d\Omega = \\ &= \int_{\Omega} (\mathrm{Tr}[\sigma^m\varepsilon^m] + \mathrm{Tr}[\sigma^M\varepsilon^M])d\Omega \end{aligned} \quad (5)$$

These relations, (4) and (5), allow us to define the micro and macro dual unknowns: $\sigma^M$ and $\underline{F}^M$ are the macro stress and efforts, and $\sigma^m$ and $\underline{F}^m$ are the micro correction stress and the micro correction efforts.

## 4 A LATIN approach

Within the linear elasticity case, the LATIN framework is transformed since time no longer plays a role. The usual duality is now related to energy, i.e. it occurs between displacements and efforts, and no longer between celerity and efforts, as this is the case with the dissipation duality for non-linearities or for dynamic problems.

Since we are dealing with perfect interfaces, let us express their behaviour: the displacement field on the interfaces must be continuous, $\underline{W} = \underline{W}'$, and the efforts must be balanced,

($\underline{F} + \underline{F}', \underline{W}^*) = 0$. With $\underline{W}^* = \underline{W}^{m*} + \underline{W}^{M*}$, this equilibrium can be rewritten in the micro-macro form:

$$\forall \underline{W}^{m*} \in W_{E,0}^m, \quad (\underline{F}^m + \underline{F}^{lM}, \underline{W}^{m*}) = 0 \quad (6)$$

$$\forall \underline{W}^{M*} \in W_{E,0}^M, \quad (\underline{F}^M + \underline{F}^{lM}, \underline{W}^{M*}) = 0 \quad (7)$$

while the displacement continuity is:

$$\underline{W}^m = \underline{W}^{lM} \quad \text{and} \quad \underline{W}^M = \underline{W}^{lM}$$

In order to break the globality of the reference problem, (1) and (2), equations for the substructured problem are separated into two groups. We intend to keep the macro-scale problem global over the whole structure, while treating the micro-scale equilibrium on each substructure separately.

For the version presented herein, the problem now consists of finding $\underline{s} = (\underline{U}, \underline{W}; \sigma, \underline{F})$ which satisfies:

* an initial group of equations ($\Gamma$) with possibly non-linear, but local-in-space, equations:

- for each substructure $\Omega^E$, Hooke's law (3), and

- for each interface $\Gamma^{EE'}$, the corresponding behaviour for the micro scale: equilibrium (6) and continuity of displacement $\underline{W}^m = \underline{W}^{lM}$.

* a second group ($\underline{A}_d$) with possibly global, but linear, equations, in which the splitting into the two levels had been previously performed: for each substructure,

- kinematic admissibility:

$$\begin{aligned} \varepsilon^m &= \varepsilon(\underline{U}^m), & \underline{U}^m|_{\partial\Omega^E} &= \underline{W}^m, \\ \underline{U}^m &\in \underline{U}_E^m, & \underline{W}^m &\in \underline{W}_E^m \end{aligned} \quad (8)$$

$$\begin{aligned} \varepsilon^M &= \varepsilon(\underline{U}^M), & \underline{U}^M|_{\partial\Omega^E} &= \underline{W}^M, \\ \underline{U}^M &\in \underline{U}^M, & \underline{W}^M &\in \underline{W}^M \end{aligned} \quad (9)$$

- static admissibility:

$$\begin{aligned} \forall \underline{U}^{m*} \in \underline{U}_{E,0}^m, & \quad \int_{\Omega^E} \mathrm{Tr}[\sigma^m \varepsilon(\underline{U}^{m*})] d\Omega = \\ &= \int_{\Omega^E} \underline{f}_d \cdot \underline{U}^{m*} d\Omega + \int_{\partial\Omega^E} \underline{F}^m \cdot \underline{W}^{m*} dS \end{aligned} \quad (10)$$

$$\begin{aligned} \forall \underline{U}^{M*} \in \underline{U}_0^M, & \quad \sum_E \int_{\Omega^E} \mathrm{Tr}[\sigma^M \varepsilon(\underline{U}^{M*})] d\Omega = \\ &= \sum_E \int_{\Omega^E} \underline{f}_d \cdot \underline{U}^{M*} d\Omega + \int_{\partial\Omega^E} \underline{F}^M \cdot \underline{W}^{M*} dS \end{aligned} \quad (11)$$

and for each interface, the macro-scale behaviour: equilibrium (7) and displacement continuity $\underline{W}^M = \underline{W}^{lM}$.

A two-stage iterative algorithm can now be built [7], see Figure 4. Its main feature is to successively produce an element of $\Gamma$ and an element of $\underline{A}_d$ at each iteration. Each stage involves a search direction, $\underline{E}^+$ and $\underline{E}^-$. These are the parameters of the method.
<figure><img src="image_6.png" /><figcaption>Figure 4: 2-stage iterative LATIN scheme</figcaption></figure>

One can notice that if the continuity of additive micro dis-placements throughout interfaces had been placed into the A<sub>d</sub> group, while keeping the additive micro equilibrium on inter-face in Γ, this would have constituted a primal approach (each solution **s** ∈ A<sub>d</sub> satisfying the continuity of displacement). In contrast, the equilibrium in A<sub>d</sub> and the continuity in Γ would have produced a dual approach. Keeping both equilibrium and continuity in Γ as an interface behaviour leads to a mixed sub-structuring approach.

4.1 Local stage

Once **s** is known, the problem then is to find **\hat{s}** ∈ Γ such that **\hat{s}** − **s** belongs to a search direction **E**<sup>+</sup>, i.e. for linear elasticity, we use the Hooke’s operator:

$$
(\hat{\sigma} - \sigma) + \mathbf{K}(\hat{\varepsilon} - \varepsilon) = 0 \quad (12)
$$

for each substructure, and:

$$
\begin{equation}
\begin{aligned}
& \forall \underline{W}^{m\star} \in \mathcal{W}_E^m, \\
& \quad (\underline{\hat{F}}^m - \underline{F}^m, \underline{W}^{m\star}) - (\mathbf{k}^m(\underline{\hat{W}}^m - \underline{W}^m), \underline{W}^{m\star}) = 0, \\
& \quad (\underline{\hat{F}}'^m - \underline{F}'^m, \underline{W}^{m\star}) - (\mathbf{k}^m(\underline{\hat{W}}'^m - \underline{W}'^m), \underline{W}^{m\star}) = 0
\end{aligned}
\tag{13}
\end{equation}
$$

for each interface. Of course, **\hat{s}** must also satisfy Hooke's law for each substructure (3) and the micro-scale behaviour of each interface.

It can be noticed that micro and macro quantities are not separated in the search direction (12) in order to ensure cou- pling between the two scales. For the boundary fields, the search direction (13) is expressed separately on the micro scale. One characteristic of the parameter **k**<sup>m</sup>, when using a multi-level approach, is its relationship with an interface char- acteristic [3] and no longer with a structural characteristic, as it was the case for the former mono-level approach, [9].

Since the search direction **E**<sup>+</sup> is local and linear, the re-
sulting local stage problem is also local in terms of the space
variable. It can be easily parallelised once the interfaces have
been mapped onto the available processors. For instance, with
the interface behaviour and the search direction (13), the addi-
tive micro displacement and efforts can be explicitly expressed
as:

$$
\begin{align*}
(k^m \underline{\hat{W}}^m, \underline{W}^{m\star}) &= (k^m \underline{\hat{W}}'^m, \underline{W}^{m\star}) = \\
&= \frac{1}{2} \left[ (k^m (\underline{W}^m + \underline{W}'^m), \underline{W}^{m\star}) - (\underline{F}^m + \underline{F}'^m, \underline{W}^{m\star}) \right]
\end{align*}
$$

$$
\begin{align*}
(\underline{\hat{F}}^m, \underline{W}^{m\star}) &= -(\underline{\hat{F}}'^m, \underline{W}^{m\star}) = \\
&= \frac{1}{2} \left[ (\underline{F}^m - \underline{F}'^m, \underline{W}^{m\star}) - (k^m(\underline{W}^m - \underline{W}'^m), \underline{W}^{m\star}) \right]
\end{align*}
$$

4.2 Linear stage

Once $\hat{s}$ is known, the problem then is to find $\mathbf{s} \in \mathbf{A}_d$, in such a
way that $\mathbf{s} - \hat{s}$ belongs to the search direction $\mathbf{E}^{-}$, i.e. for each
substructure:

$$
(\sigma - \hat{\sigma}) - \mathbf{K}(\varepsilon - \hat{\varepsilon}) = 0 \quad (14)
$$

and for each interface:

$$
\begin{equation}
\begin{array}{@{}l@{}}
\forall \underline{W}^{m\star} \in \mathcal{W}_E^m, \\
\qquad (\underline{F}^m - \underline{\hat{F}}^m, \underline{W}^{m\star}) + (\underline{k}^m(\underline{W}^m - \underline{\hat{W}}^m), \underline{W}^{m\star}) = 0
\end{array}
\tag{15}
\end{equation}
$$

As previously noted, **K** does not separate out micro and macro quantities in (14), although **k**<sup>m</sup> does so in (15).

Since **s** ∈ **A<sub>d</sub>**, it is said to be admissible. It is subjected to both kinematic (8 and 9) and static (10 and 11) admissibilities. Such a problem is both linear and global per each substruc- ture only. It can also be parallelised when the substructures get distributed among the processors.

We are now going to describe more precisely the kind of
problems to be solved during the linear stage. Taking into ac-
count the properties that ˆs has to satisfy, the search direction
(14) leads to σ = Kε, i.e.:

$$
\forall \varepsilon^{*} = \varepsilon^{m*} + \varepsilon^{M*},
\quad
\int_{\Omega^{E}} \mathrm{Tr}[\sigma \varepsilon^{*}] d\Omega = \int_{\Omega^{E}} \mathrm{Tr}[\varepsilon \mathbf{K} \varepsilon^{*}] d\Omega \quad (16)
$$

4.2.1 Micro-scale problem

Let's look more closely now at the micro-scale problem, in
which all macro quantities are supposed to be known. One
must find **s**<sup>*m*</sup> = (**U**<sup>*m*</sup>, **W**<sup>*m*</sup>; **σ**<sup>*m*</sup>, **F**<sup>*m*</sup>) which satisfies, for each sub-
structure, the kinematic and static admissibilities (8) and (10),
as well as the search direction (14). The later then leads to:

$$
\int_{\Omega^E} \mathrm{Tr}[\sigma^m \varepsilon^{m\star}] d\Omega = \int_{\Omega^E} \mathrm{Tr}[(\mathbf{K}\varepsilon^m + \mathbf{K}\varepsilon^M)\varepsilon^{m\star}] d\Omega
$$

Using the static admissibility and the search direction, the
resulting displacement-oriented formulation consists of finding
on each substructure **U**<sup>*m*</sup> and **W**<sup>*m*</sup> such that:

$$
\begin{align*}
\varepsilon^m &= \varepsilon(\underline{U}^m), & \underline{U}^m|_{\partial\Omega^E} &= \underline{W}^m, \\
\underline{U}^m &\in \mathcal{U}_E^m, & \underline{W}^m &\in \mathcal{W}_E^m
\end{align*}
$$

and:

$$
\begin{align*}
& \forall \underline{U}^{m\star} \in \mathcal{U}_{E,0}^{m}, \\
& \qquad \int_{\Omega^E} \mathrm{Tr}[\varepsilon(\underline{U}^m) \mathbf{K} \varepsilon(\underline{U}^{m\star})] + \int_{\partial\Omega^E} \underline{U}^m \cdot \underline{k}^m \underline{U}^{m\star} dS = \\
& = \int_{\Omega^E} \underline{f}_d \cdot \underline{U}^{m\star} d\Omega + \int_{\partial\Omega^E} (\underline{\hat{F}}^m + \underline{k}^m \underline{\hat{W}}^m) \cdot \underline{U}^{m\star} dS + \\
& \phantom{=} - \int_{\Omega^E} \mathrm{Tr}[\varepsilon^M \mathbf{K} \varepsilon(\underline{U}^{m\star})] d\Omega \quad (17)
\end{align*}
$$

The discretised problem is then obtained when expressing
the displacement with the classical set of basis functions N:
**U** = N[U], where [U] is the column of unknowns. If [U<sup>m</sup>] and
$[U^M]$ denote degrees of freedom for the additive micro and for the macro quantities, the discretised problem can be written as:

$$([K^m] + [k^m])[U^m] = [f_d^m] + [\hat{F}^m] - [K^{MM}][U^M] \quad (18)$$

$[K^m]$ is a part of the classical rigidity matrix generated from the constitutive relation; $[k^m]$ is a boundary rigidity arising from the search direction. Both are constant throughout the iterations, and are therefore factorised once during an initialisation phase. $[f_d^m]$ is the prescribed loading; $[\hat{F}^m]$ is a boundary loading due to $\hat{F}^m + k^m \hat{W}^m$, which depends on the iteration; and $[K^{MM}][U^M]$ is a body loading generated from the macro-scale displacement $U^M$. $[K^{MM}][U^M]$ is an out-of-diagonal block in the classical rigidity matrix, expressed on the previous h-basis split, when splitting the degrees of freedom in micro and macro parts:

$$[K] = \begin{bmatrix} K^m & K^{MM} \\ K^{Mm} & K^M \end{bmatrix}$$

In order to express the macro-scale problem, in the next section, the structure of the micro-scale problem must be detailed. Since the problem is linear, the solution can be expressed as:

$$\underline{U}^m = \hat{U}_d + \tilde{U}$$

$\tilde{U}$ is the part of the solution for which $\varepsilon^M$ is the only loading term on the right-hand side of (17). Let $\tilde{U} \in U^M$ be such that:

$$\forall \underline{U}^{M*} \in \mathcal{U}^M, \quad \int_{\Omega^E} \mathrm{Tr}[\varepsilon(\tilde{\underline{U}}) \mathbf{K} \varepsilon(\underline{U}^{M*})] d\Omega = \\ = \int_{\Omega^E} \mathrm{Tr}[\varepsilon(\tilde{\underline{U}}) \mathbf{K} \varepsilon(\underline{U}^{M*})] d\Omega$$

i.e. from a discretised point of view:

$$[K^M][\tilde{\underline{U}}] = [K^{Mm}][\tilde{\underline{U}}]$$

Since $\tilde{\underline{U}} \in U_E^m$ and $\tilde{\underline{U}} \in U^M$, the later is energetically equivalent to the response in the micro scale to a macro-scale loading. Finally, $\tilde{\underline{U}}$ is derived from $\underline{U}^M$ using $\tilde{\underline{U}}$. It can be formally written as:

$$\mathbf{K}\varepsilon(\tilde{\underline{U}}) = \mathbf{L}\mathbf{K}\varepsilon(\underline{U}^M) \quad (19)$$

or:

$$[K^M][\tilde{\underline{U}}] = -[K^{Mm}][K^m] + [k^m]^{-1}[K^{MM}][U^M]$$

### 4.2.2 Macro-scale problem

This problem is treated similarly with the search direction (16), which then leads to:

$$\int_{\Omega^E} \mathrm{Tr}[\sigma^M \varepsilon^{M*}] d\Omega = \int_{\Omega^E} \mathrm{Tr}[(\mathbf{K}\varepsilon^M + \mathbf{K}\varepsilon^m)\varepsilon^{M*}] d\Omega$$

Using expression (19) in the macro-scale equilibrium, the displacement-oriented formulation now consists of finding $\underline{U}^M$ and $\underline{W}^M$, such that:

$$\varepsilon^M = \varepsilon(\underline{U}^M), \quad \underline{U}^M|_{\partial\Omega^E} = \underline{W}^M, \\ \underline{U}^M \in \mathcal{U}^M, \quad \underline{W}^M \in \mathcal{W}_E^M$$

and:

$$\begin{aligned} \forall \underline{U}^{M*} \in \mathcal{U}_0^M, \\ \sum_E \int_{\Omega^E} \mathrm{Tr}[\varepsilon(\underline{U}^M)(1+\mathbf{L})\mathbf{K}\varepsilon(\underline{U}^{M*})] d\Omega = & \\ = \sum_E \int_{\Omega^E} \underline{f}_d \cdot \underline{U}^{M*} d\Omega + & \\ & \int_{\partial\Omega^E} \underline{F}^M \cdot \underline{U}^{M*} dS + \\ & - \sum_E \int_{\Omega^E} \mathrm{Tr}[\varepsilon(\hat{\underline{U}}_d)\mathbf{K}\varepsilon(\underline{U}^{M*})] d\Omega \quad (20) \end{aligned}$$

The regularity constraints for $\underline{U}^M$ lead to a global problem over the whole structure for the macro scale. Since only perfect interfaces are used in this example, the interface loading in the global macro-scale equilibrium (20) is then:

$$\begin{aligned} \sum_E \int_{\partial\Omega^E} \underline{F}^M \cdot \underline{U}^{M*} dS = & \\ = \sum_{\Gamma^{EE'}} \int_{\Gamma^{EE'}} (\underline{F}^M + \underline{F}'^M) \cdot \underline{U}^{M*} dS = 0 \end{aligned}$$

for the version used herein, as an initial implementation of the micro-macro scheme.

This macro-scale problem is used as the mechanism to propagate the information globally at each iteration.

When discretised, the macro-scale problem becomes:

$$[\bar{K}][U^M] = [\bar{F}] \quad (21)$$

with the condensed rigidity and right-hand side:

$$[\bar{K}] = [K^M] - [K^{Mm}][[K^m] + [k^m]]^{-1}[K^{MM}] \\ [f_d^M] - [K^{Mm}][\hat{U}_d]$$

$(1+\mathbf{L})\mathbf{K}\varepsilon(\underline{U}^M)$ is the homogenised stress and $(1+\mathbf{L})\mathbf{K}$ is the homogenised Hooke's tensor. $\underline{U}^M$ and $\sigma^M$ are then the effective displacement and stress fields.

The solution $\underline{U}^M$ now allows us to solve the micro-scale problem (17) as a second problem.

## 5 Example

The reported results are intended to enable comparing the proposed micro-macro approach with the initial version of the multi-level scheme [3], and with the initial mono-level algorithm [9, 2]. In particular, the initial multi-level approach does not take into account the homogenised part $\mathbf{L}$, but it does retain the coupling terms in both the micro and the macro-scale equilibrium, and solves the problem in a fixed-point manner together with the LATIN iterations. In contrast, the proposed micro-macro approach uses a homogenisation technique performed with a static condensation on the macro-scale problem.

The general layout of the micro-macro approach is described in Table 1, and has been implemented in the industrial-type code CASTEM 2000 [13], developed at the CEA in Saclay, France. All of the computations have been performed on an SGI ORIGIN 2000 with 32 processors and 8 Gb of shared memory.

In order to check the performance of the proposed approach, let us consider the model problem of a slendered bidimensional structure submitted to a parabolic bending loading
<table><caption>Table 1: Micro-macro LATIN algorithm</caption><thead><tr><th>Macro scale<br>dedicated to 1 processor</th><th>Micro scale<br>dedicated to n processors</th></tr></thead><tbody><tr><td><strong>Initialisation</strong><br>receiving [K&#773;<sub>E</sub>]<br>assembling the substructure contributions<br>factorisation of [K&#773;]</td><td>Initialisation<br>compute elementary rigidities [K<sub>E</sub>] and [k<sub>E</sub><sup>m</sup>]<br>assembling the element contributions<br>factorisation of [K<sub>E</sub>] + [k<sub>E</sub><sup>m</sup>] as well as<br>condensation to obtain [K&#773;<sub>E</sub>]<br>sending [K&#773;<sub>E</sub>]<br>initialisation of s&#770; = 0</td></tr><tr><td><strong>Loop over iterations</strong><br>linear macro stage<br>receiving [F&#773;<sub>E</sub>]<br>assembling the substructure contributions<br>forward-backward on (21): global problem<br>sending [U<sup>M</sup>]</td><td>Loop over iterations<br>linear micro stage<br>compute elementary loading<br>assembling the element contributions<br>forward substitution on (18) as well as<br>condensation to obtain [F&#773;<sub>E</sub>]<br>sending [F&#773;<sub>E</sub>]<br>receiving [U<sup>M</sup>]<br>compute coupling term [K<sub>E</sub><sup>mM</sup>][U<sub>E</sub><sup>M</sup>]<br>backward substitution to obtain [U<sub>E</sub><sup>m</sup>]<br>local micro stage<br>convergence check</td></tr></tbody></table>

(see Figure 5). For this example, the finite element meshes of the two scales are embedded. The reference $(\underline{U}; \sigma)_{\text{ref}}$ here is the direct finite element solution without decomposition. This allows us to define the convergence rate in the energy norm:

$$ \tau = -\log \frac{e_{n+1}}{e_n} $$

where:

$$ e_n^2 = \frac{\frac{1}{2} \int_{\Omega} \text{Tr}[(\sigma_n - \sigma_{\text{ref}})\mathbf{K}^{-1}(\sigma_n - \sigma_{\text{ref}})]d\Omega}{\frac{1}{2} \int_{\Omega} \text{Tr}[\sigma_{\text{ref}}\mathbf{K}^{-1}\sigma_{\text{ref}}]d\Omega} $$

<figure><img src="image_7.png" /><figcaption>Figure 5: Model problem and its decompositions into 4, 16 and 64 substructures and interfaces</figcaption></figure>

For the following computations, the size of the problem local to one substructure is kept constant. The size of the refer-

ence problem increases as well as the number of substructures. Table 2 shows the various characteristics of the solved problems.

Figure 6 presents the averaged convergence rate (up to a high level of convergence: $e_n \le 0.1\%$) versus the number of substructures. The relative convergence rate is obtained by normalising it with the convergence rate for the case of 4 substructures and the mono-level LATIN algorithm. Table 3 reports the numerical values. It illustrates a well-known behavioural characteristic of domain decomposition methods without a multi-level scheme: slowing the convergence rate when increasing the number of subdomains [1].

<figure><img src="image_8.png" /><figcaption>Figure 6: Relative convergence rate versus number of substructures</figcaption></figure>

The evolution in the stress and displacement fields throughout the iterations are compared in Figure 7 for the mono-level
<table><caption>Table 2: Characteristics of the problems involved (micro scale is meshed with 6-node triangles, macro scale with 3-node triangles and 2 elements per substructure)</caption><thead><tr><th colspan="3">substructuring (total)</th><th colspan="2">macro scale problem</th><th colspan="3">direct computation</th></tr><tr><th>nb sstr</th><th>nb dof</th><th>memory requirement for factorised matrix</th><th>nb dof</th><th>memory requirement for factorised matrix</th><th>nb dof</th><th>nb elements</th><th>memory requirement for factorised matrix</th></tr></thead><tbody><tr><td>4</td><td>8 712</td><td>3.3 Mb</td><td>20</td><td>0.7 Kb</td><td>8 514</td><td>2048</td><td>4.7 Mb</td></tr><tr><td>16</td><td>34 848</td><td>13 Mb</td><td>54</td><td>3 Kb</td><td>33 410</td><td>8192</td><td>39.7 Mb</td></tr><tr><td>64</td><td>139 392</td><td>52 Mb</td><td>170</td><td>14 Kb</td><td>132 354</td><td>32768</td><td>324 Mb</td></tr></tbody></table>

<table><caption>Table 3: Values of relative convergence rates</caption><thead><tr><th rowspan="2">number of substructures</th><th colspan="3">relative convergence rate</th></tr><tr><th>mono-level</th><th>2-level</th><th>micro-macro</th></tr></thead><tbody><tr><td>4</td><td>1</td><td>1.41</td><td>2.15</td></tr><tr><td>16</td><td>0.25</td><td>1.10</td><td>2.02</td></tr><tr><td>64</td><td>0.10</td><td>1.20</td><td>2.05</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 7: Comparison of solutions (Mises stress field and displacement field) at iterations 1,2 and 15 for: a) on the left: a mono-level approach; and b) on the right: a multi-level approach.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: a) on the left: the micro scale solution (with a 10-magnification for displacement); b) on the right: the macro scale solution.</figcaption></figure>
and multi-level approaches (here, the initial version), at itera-
tions 1, 2 and 15. For the latter approach, both the micro and
macro-scale parts of the solution are shown in Figure 8.

Last, it should be pointed out that for the mono-level ap-
proach, the search direction is **k** = *E*/*L*₀, where *E* is Young's
modulus and *L*₀ has an optimum value of twice the length of
the structure, [9]. When using the multi-level or micro-macro
approach, the value of **k**<sup>*m*</sup> = *E*/*L*₀ has been set such that *L*₀ is
equal to 0.25 times the length of an interface. It is no longer
characterised by the global behaviour of the structure, [4].

With a global macro-scale problem, the algorithm also be-
comes numerically scalable: the convergence rate is quasi-
independent of the number of substructures.

The micro-macro approach deals with a structural ho-
mogenisation: even for a problem with a homogeneous ma-
terial, the multi-level approach using the “homogenised” op-
erator (1 + **L**) **K** for the large scale problem has an improved
convergence rate when compared to a classical coarse discreti-
sation.

6 Conclusion

The objective herein is first to build a mechanical and parallel
approach related to domain decomposition methods, in order to
use parallel architecture computers for problems with a large
number of degrees of freedom. It deals with a mixed substruc-
turing method since the unknowns are both the displacement
and the efforts on the interfaces.

Using a macro-scale problem to propagate the informa-
tion throughout the substructures globally improves the con-
vergence rate; as this rate becomes quasi-independent of the
number of substructures, the resulting algorithms are numeri-
cally scalable. Moreover, the “optimal” search direction is then
related to an interface characteristic, and does not depend any
longer on the behaviour of the whole structure.

The proposed technique belongs to the field of structural
analysis with homogenisation. It combines both micro and
macro representations of the problem into the LATIN frame-
work. Since the macro scale problem remains global over
the whole structure, this leads to an information exchange
mechanism for all of the substructures. This approach does
not require any specific treatment for boundary-located areas
and concerns structures that are not required to be a periodic
medium.

Even with homogeneous material, the use of this structural-
homogenised like operator for the global macro-scale, im-
proves the convergence rate. As this operator is relatively easy
to build, the additive cost is expected to be small enough to
have an improved algorithm in term of CPU cost as well. Of
course, improved performances are to be expected when deal-
ing with highly heterogeneous structures as this is the case of
reinforced or composite structures.

References

[1] J. H. Bramble, J. E. Pasciak, and A. H. Schatz. The con-
struction of preconditioners for elliptic problems by sub-
structuring. *I. Math. Comp.*, 47(175):103–134, 1986.

[2] L. Champaney, J.-Y. Cognard, D. Dureisseix, and P. Ladevèze. Large scale applications on parallel computers of a mixed domain decomposition method. *Computational Mechanics*, (19):253–263, 1997.

[3] D. Dureisseix and P. Ladevèze. Parallel and multi-level strategies for structural analysis. In J.-A. Désidéri, editor, *Proceedings of the Second European Conference on Numerical Methods in Engineering*, pages 599–604. Wiley, September 1996.

[4] D. Dureisseix and P. Ladevèze. A multi-level and par-
allel strategy for structural analysis. In D. R. J. Owen,
E. Onate, and E. Hinton, editors, *Proceedings of the
Fifth Int. Conf. on Computational Plasticity*, pages 1790–
1797, Barcelona, March 1997. CIMNE.

[5] C. Farhat and F.-X. Roux. A method of finite element tearing and interconnecting and its parallel solution algorithm. *International Journal for Numerical Methods in Engineering*, 32:1205–1227, 1991.

[6] J. Fish and V. Belsky. Multigrid method for periodic het-
erogeneous media. Part 1 and part 2. *Computer Methods
in Applied Mechanics and Engineering*, 126:1–38, 1995.

[7] P. Ladevèze. *Mécanique non-linéaire des structures — Nouvelle approche et méthodes de calcul non incrémentales*. Springer Verlag, 1998. English version.

[8] P. Ladevèze and D. Dureisseix. Une nouvelle stratégie de calcul parallèle et micro / macro en mécanique non-linéaire. Rapport Interne n° 188, Laboratoire de Mécanique et Technologie, Cachan, July 1997.

[9] P. Ladevèze and P. Lorong. A large time increment approach with domain decomposition technique for mechanical non linear problems. In R. Glowinski, editor, *Comput. Meths. Appl. Sc. Engng.*, pages 569–578, New York, 1992. INRIA, Nova Science.

[10] P. Le Tallec. Domain decomposition methods in com-
putational mechanics. In *Computational Mechanics Advances*, volume 1. North-Holland, 1994.

[11] J. Mandel. Balancing domain decomposition. *Communications in Applied Numerical Methods*, 9:233–241, 1993.

[12] A. K. Noor. New computing systems and future high performance computing environment and their impact on structural analysis and design. *Computers & Structures*, 64(1–4):1–30, 1997.

[13] P. Verpeaux, T. Charras, and A. Millard. CASTEM 2000 : une approche moderne du calcul des structures. In J.-M. Fouet, P. Ladevèze, and R. Ohayon, editors, *Calcul des Structures et Intelligence Artificielle*, volume 2, pages 261–271. Pluralis, 1988.

[14] H. Yserentant. On the multi-level splitting of finite ele-
ment spaces. *Num. Math.*, 49:379–412, 1986.

[15] O. C. Zienkiewicz, J. P. De S. R. Gago, and D. W. Kelly. The hierarchical concept in finite element analysis. *Computers & Structures*, 16(1–4):53–65, 1983.