 Research article
 Open Access
 Published:
Noninvasive global–local coupling as a Schwarz domain decomposition method: acceleration and generalization
Advanced Modeling and Simulation in Engineering Sciences volume 5, Article number: 4 (2018)
Abstract
The noninvasive global–local coupling algorithm is revisited and shown to realize a simple implementation of the optimized nonoverlapping Schwarz domain decomposition method. This connection is used to propose and compare several acceleration techniques, and to extend the approach to non conforming meshes.
Introduction
The noninvasive local–global coupling technique proposed by Allix [1] is an iterative method which aims at making accurate the well known submodeling technique [2,3,4]. It is strongly related to many reanalysis techniques [5,6,7] and domain decomposition methods [8].
The aim of this technique is to evaluate the effect of local modifications inside a computational model (geometry, material and load) without requiring heavy developments. More precisely the objective is to use an industrial model with a given commercial software and to simulate the presence of local alterations by iteratively spawning computations with only extra traction loads inside the model. Moreover, the alterations can be computed on any chosen software including dedicated research codes.
This philosophy was successfully applied in many different contexts like: the introduction of local plasticity and geometrical refinements [1], the computation of the propagation of cracks in a sound model [9], the evaluation of stochastic effects with deterministic computations [10], the taking into account of the exact geometry of connectors in an assembly of plates [11]. In [12] the method was used in order to implement a nonlinear domain decomposition method [13,14,15,16] in a noninvasive manner with Code_aster. The extension of the approach to explicit dynamics was proposed in [17], improved in [18] and applied to the prediction of delamination under impact in [19]. Alternative noninvasive strategies can be derived from the extended finite element method [20, 21].
After a description of the method (“Derivation of the noninvasive algorithm” section), this paper provides several contributions. First the noninvasive coupling algorithm is proved to realize a simple implementation of the optimized nonoverlapping Schwarz domain decomposition method (“Connexion with alternate nonoverlapping Schwarz method” section). Several accelerations techniques are proposed (“Analysis and acceleration of the global/local algorithm” section), some are classical but the linear and nonlinear conjugate gradient is new in this framework. The algorithms are described in a very programmerfriendly manner. Last, an overlapping version of the method is proposed (“Overlapping version” section) which can be used to handle fully nonconforming meshes.
Derivation of the noninvasive algorithm
The algorithm we study is very general and applies to the study of many PDEs; in order to fix the ideas, we consider problems of nonlinear quasistatic structure mechanics under the small strain hypothesis. We note u the displacement field, \(\varepsilon \) the symmetric part of the gradient, \(\sigma \) the Cauchy stress tensor. The domain \(\Omega \) is submitted to given body force f, Dirichlet condition \(u_d\) on the part \(\partial _d \Omega \) of the boundary and Neumann condition g on the complement part \(\partial _n \Omega \). In order to manage viscous materials, the study is conducted over a time interval \(\mathcal {T}=[0,T]\), and the following equations are meant to be satisfied at any time \(t\in \mathcal {T}\), which we omit to write except when necessary.
Let \(V(\Omega )= \left\{ v\in H^1(\Omega ), v=u_d \text { on }\partial _d\Omega \right\} \) be the affine space of admissible displacement fields and \(V^0(\Omega )\) the associated vector space. The conservation of momentum can be written as:
The notation \(\sigma (u)\) stands for local or non linear constitutive laws defined under the following functional expression:
This modeling of the mechanical behavior is typically suited for elastoviscoplastic materials. For most models an alternative description by internal variables summarizing the effect of the past history can be found.
The mechanical problem above takes the following classical form:
where l is a continuous linear form, and a is a continuous coercive form, linear in the second variable, note that a may be nonlinear in the first variable.
In the following we handle several space subdomains and models, when any quantity is specifically attached to one model, a superscript mentions it. The domains are illustrated on Fig. 1.
The Reference problem (superscript \({}^R\)) is set on the domain \(\Omega {}^R\) which is the assembly of two nonoverlapping subdomains: the zone of interest where a Fine model is required for a reliable simulation (superscript \({}^F\)), and a Complement zone (superscript \({}^C\)) where a simpler model is sufficient (and which in general covers most of the structure). The interface is \(\Gamma =\partial \Omega {}^C\cap \partial \Omega {}^F\), it is thus immersed in \(\Omega {}^R\). Note that using several zones of interest presents no difficulty as long as they do not overlap. The Reference problem is defined as the assembly of the Fine and Complement subproblems; it is never formed in practice:
Note that \(l{}^F\) is associated with the restriction of body force to \(\Omega {}^F\) and of Neumann traction to \(\partial _n\Omega {}^R\cap \partial \Omega {}^F\).
We assume that we have another representation of the zone of interest, named Auxiliary representation (superscript \({}^A\)) which shares the same characteristics as the Complement zone, and which is thus coarser than the Fine representation. Typically if \(\Omega {}^F\) was a zone where material coefficients had strong variations, the Fine representation would follow the exact distribution whereas the Auxiliary representation could use a homogenized behavior. An application is the case where the Fine model is stochastic whereas the Auxiliary model is deterministic [10]. The load could also be simplified: \(l{}^A\) is associated with body force applied to \(\Omega {}^A\) and traction applied to \(\partial \Omega {}^A\setminus \Gamma \). We insert the Auxiliary representation of the zone of interest in the Reference problem:
The Global problem, (superscript \({}^G\)), is the assembly of the Complement zone with the Auxiliary (coarse) representation of the zone of interest, this problem is in practice assembled and dealt with by commercial software.
From the previous equation, we could derive the following stationary iteration:
which would correspond to a fixed point of the Reference problem preconditioned by the coarse Global system. Not only convergence would be slow but also the righthandside terms would not be easy to compute in practice. Moreover, this iteration needs the Auxiliary domain \(\Omega {}^A\) to be coincident with the Fine domain \(\Omega {}^F\) which is a limitation we want to get rid of. In the following, we only assume that the interface is on the boundary of the Auxiliary domain \(\Gamma \subset \partial \Omega {}^A\), so that \(\Gamma \) is on the boundary of all subdomains. Note that no condition on \(\Gamma \) is imposed to the fields belonging to \(V(\Omega {}^X)\) and \(V^0(\Omega {}^X)\) (\(X\in \{A,F,C\}\)), the Dirichlet conditions on the subdomains are imposed with Lagrange multipliers. We note \(V_\Gamma \) the trace space of displacements on \(\Gamma \) and \(V_\Gamma ^*\) its dual space of interface tractions, \(\langle \lambda ,v\rangle _\Gamma \) is the associated duality bracket with \(\lambda \in V_\Gamma ^*\) and \(v\in V_\Gamma \).
We thus choose to associate the right hand side of (5) with the evaluation of local problems. Starting from \(p_0=0\) (the mathematical space for \(p_n\) is discussed later), the basic global/local iteration is then the following:
In words, the Global problem is the coarse problem with extra load p, the Fine and Auxiliary systems are resolutions on the domain of interest with imposed Dirichlet conditions on \(\Gamma \). We chose a Lagrangian formulation for these problems in order to make appear the reaction forces \(\lambda {}^F\) and \(\lambda {}^A\). The update is simply the equivalent of (6) with fields issuing from the local solves instead of the global one.
Remark 1
Of course, the Lagrange multipliers are equal to the normal stress:
where \(X\in \{A,F\}\) and \(n^X\) is the outward normal vector. \(\square \)
We have the following properties:

Assuming the fine and auxiliary problems were solved exactly, we have:
$$\begin{aligned} p_{n+1} = \left( \lambda {}^A_{n}\lambda {}^F_{n} \right) \in V_\Gamma ^* \end{aligned}$$(9)the corrective load p is then an immersed surface traction. In the following, we always assume the exactness of the computations; note that using inexact solvers was investigated in [22] where the method is identified with a localized multigrid iteration.

Because the Auxiliary problem corresponds to the restriction of the Global problem on the zone of interest with global displacement imposed, we directly have:
$$\begin{aligned} u{}^A_{n}= u{}^G_{n\Omega {}^A} \end{aligned}$$(10)The introduction of the Auxiliary problem is thus not mandatory, it is just a workaround in case of software unable to compute the reaction in an immersed surface. Of course, the Auxiliary problem can be solved in parallel with the Fine problem.

We can also define the reaction from the Complement zone for a given \(u{}^G_{n}\):
$$\begin{aligned} \langle \lambda {}^C_{n} , v \rangle _\Gamma = a{}^C(u{}^G_{n},v)  l{}^C(v), \quad \forall v \in V(\Omega {}^C) \end{aligned}$$(11)Then we see that:
$$\begin{aligned} \lambda {}^C_{n} + \lambda {}^A_{n} = p_{n} \end{aligned}$$(12)The surface traction \(p_{n}\) generates a discontinuity in the normal stress of the Global problem.

If we replace the auxiliary reaction by the complement one, we have:
$$\begin{aligned} p_{n+1} = p_n + r_{n} \text { with } r_{n+1} = \left( \lambda {}^F_{n}+\lambda {}^C_{n}\right) \end{aligned}$$(13)in words, the correction brought to \(p_{n+1}\) corresponds to the lack of balance between the Complement zone and the Fine representation of the zone of interest. This lack of balance is the residual r of the algorithm. The algorithm converges when the two representations are in equilibrium (\(r=0\), in which case the extra load p shall not evolve anymore).

The algorithm makes no use of domain integrals to communicate between subdomains; only interface data (on \(\Gamma \)) are exchanged, namely the displacement \(u{}^G\) and the reactions \(\lambda {}^F\) and \(\lambda {}^A\) (or \(\lambda {}^C\)). As long as the interface \(\Gamma \) is well represented in all models, it is not necessary to use the exact Fine domain \(\Omega {}^F\) in the Auxiliary problem, any coarser representation is possible (\(\Omega {}^A\)). Typically microperforations or micro cracks need not be represented in the Auxiliary problem. Of course modifying the representation of the zone of interest may have consequences on the convergence of the algorithm (but not on its limit which is the reference solution).
Connexion with alternate nonoverlapping Schwarz method
The question of linking the noninvasive global–local coupling method to the many variants of domain decomposition and associated algorithms, like chimera, was studied in other publications like [8]. Here we propose to connect the method with the iterations of a nonoverlapping optimized Schwarz method. The theoretical framework of Schwarz method will allow us natural extensions to the method, in particular the use of overlaps to treat mesh incompatibilities.
We consider two nonoverlapping subdomains, \(\Omega {}^C\) and \(\Omega {}^F\), connected by the interface \(\Gamma \). The decomposed problem to solve can be written as:
In words, subdomains must be in mechanical equilibrium while displacements shall be equal on the interface and force fluxes shall be balanced. The optimized Schwarz method consists in using Robin conditions at the interface. The Robin conditions are materialized by operators called interface impedances (or interface stiffnesses): \(Q{}^C\) and \(Q{}^F\) from \(V_\Gamma \) to \(V_\Gamma ^*\). The interface conditions are rewritten as:
where we need in particular \((Q{}^C+Q{}^F)\) to be injective for the equivalence with initial conditions to hold. In general \(Q{}^C\) and \(Q{}^F\) are chosen to be such that each associated form \(V_\Gamma ^2 \ni (u,v) \mapsto \langle Q{}^X(u), v \rangle _{\Gamma }\) is bilinear symmetric continuous coercive.
The new conditions can be combined with the equilibrium:
Hence the alternate optimized Schwarz stationary iterations (\(\lambda {}^F_0=0\), \(u{}^F_0=0\)):
It is well known that the optimal value for one subdomain’s impedance is the DirichlettoNeumann operator of the other subdomain which we note \(\mathcal {S}{}^X\). Typically for \(\Omega {}^C\), we have:
Note that for linear problems, \(\mathcal {S}{}^X\) is an affine operator (and not just a linear operator) since it also takes into account the effect of the load. The existence of this operator is conditioned to the wellposedness of the Dirichlet problem over subdomains. There are many contexts where this wellposedness can be proved, at least locally, see [23] for details. The global existence of the operator can be proved in the case of coercive continuous monotone operators; see [24, 25] for an analysis at the level of the variational formulation and [26] for the analysis of the finite element approximation. Mechanically this case is associated with positive hardening behaviors and certain contact laws, in small strains [27, 28]. Moreover, the DirichlettoNeumann operator inherits properties from the initial problem (typically monotonicity, coercivity and continuity; see [29] and associated bibliography).
The global–local algorithm corresponds to the choice \(Q{}^C=\mathcal {S}{}^A\) and formally \(Q{}^F=\infty \) (the Dirichlet condition being seen as the limit case of an infinite interface impedance). The choice \(Q{}^C=\mathcal {S}{}^A\) is extremely strong because we can expect \(\mathcal {S}{}^A\) to be a good approximation of \(\mathcal {S}{}^F\), not only in term of stiffness (\(a{}^A\) vs \(a{}^F\)) but also in term of load (\(l{}^A\) vs \(l{}^F\)) which corresponds to providing a good initialization to the algorithm.
The framework of Schwarz method enables us to recover the following features:

Krylov acceleration: replacing stationary iterations by Krylov solvers is classical in Schwarz methods [30]. The Dirichlet condition \(Q{}^F=\infty \) preserves some symmetry so that we can derive a conjugate gradient algorithm, see “Conjugate gradient” section.

Mixed approach: the condition \(Q{}^F=\infty \) is a poor approximation of the optimal choice. In [31] a twoscale approximation of \(\mathcal {S}{}^C\) was proposed for the global–local coupling.

Parallel processing: the global–local method corresponds to the alternate version of the optimized Schwarz method. The parallel version could be tried in the noninvasive context. Note that this would only make sense in the presence of multiple Fine zones with finite Fine impedance \(Q{}^F<\infty \).

Nonlinearity: stationary iterations can directly be transferred to nonlinear problems, in particular the ones with monotone operators (positive hardening) [32, 33]. The local–global method was successfully applied in many nonlinear problems like plasticity or fracture [1, 22]

Overlapping version: optimized Schwarz methods also exist with overlaps. In [34], the overlap was used as a buffer zone to dampen edge effects in plate/3D coupling. In “Overlapping version” section, we present another application, the handling of nonmatching meshes.
Analysis and acceleration of the global/local algorithm
Notations
In order to further analyze the algorithm and be more practical, we now consider the finite element discretization of the problem. We use the following notations: \(\mathbf {f}\) for the generalized forces, \(\mathbf {u}\) for the nodal displacement and \(\varvec{\lambda }\) for the nodal reactions and \(\mathbf {p}\) for the nodal component on the immersed surface effort. When indexing degrees of freedom, \(F,\ A,\ C\) stand for the internal degrees of freedom whereas \(\Gamma \) stands for nodes on the interface (whose description is identical in all models). We tried to use minimal notations, but sometimes a quantity defined on the interface is issued from one side specifically, in which case we make it clear by an extra superscript. In the linear(ized) case notation \(\mathbf {K}\) is used for the stiffness matrices.
Remark 2
We recall that the nodal reaction is not the discretization of the Lagrange multiplier. Indeed for a boundary degree of freedom i associated with shape function \(\phi _i\), we have:
where \(X\in \{C,A,F\}\), \(n^X\) is the outward normal vector and \(\sigma _h\) is the stress tensor obtained from the finite element computation. Thus the nodal reactions \(\varvec{\lambda }^X\) can be computed either by using a Lagrangian formulation for the Dirichlet condition (which is fairly common in commercial software) or by using the formula above to postprocess it from the finite element stress (which may be complex to implement in legacy software); hence the use of the Auxiliary model to compute reactions on the immersed interface. \(\square \)
If we assume that the Reference and the Global problems are wellposed, then Dirichlet problems are well posed on all subdomains, at least locally near the solution. We can then define the following nonlinear discrete DirichlettoNeumann operators \(\mathcal {S}{}^X\) (we use the same notation as in the continuous case) which compute the reactions \(\varvec{\lambda }{}^X\) from a given interface displacement \(\mathbf {u}_\Gamma \):
Because of the nonlinearity, the effects of given loads appear as parameter of the method, they will be omitted in the absence of ambiguity.
Remark 3
In the case of linear problems, it is possible to give an explicit formula for the DirichlettoNeumann operators. As an illustration, the equilibrium of the Fine problem can be written as:
which can be condensed as:
In that case, \(\mathcal {S}{}^F\) is an affine operator: \(\mathbf {S}{}^F\) is the well known Schur complement of the Fine domain on the interface. Linearity allows to set apart the contrition of the given load, with \(\mathbf {b}{}^F\) the condensed righthand side. Note that the internal displacement in the Fine domain was implicitly computed as:
\(\square \)
Because of the additivity of integral with respect to the domain, the Global operator verifies the following decomposition \(\mathcal {S}{}^G=\mathcal {S}{}^C+\mathcal {S}{}^A\) and the Reference operator writes \(\mathcal {S}{}^R=\mathcal {S}{}^C+\mathcal {S}{}^F\), and we can rephrase the Global and Reference problems in a condensed manner:
Note that each time one of the condensed operators is employed, the displacement inside the subdomains is implicitly computed: for instance, \(\mathbf {u}{}^G\) is a byproduct of (24) and \(\mathbf {u}{}^F\) in a byproduct of (20). To make it clearer, we will use notations \(\mathcal {S}{}^X\) when analyzing the methods, whereas we will use the following functional notations when describing the algorithms:

\([\mathbf {u}{}^G]=\mathtt {SolveGlobal}(\mathbf {p};\mathbf {f}{}^G)\), \(\mathbf {u}{}^G\) is defined on the whole Global model and in particular we have \(\mathbf {u}{}^G_\Gamma = \mathcal {S}{}^{G^{1}}(\mathbf {p};\mathbf {f}{}^G)\).

\([\mathbf {u}{}^F,\varvec{\lambda }{}^F]=\mathtt {SolveFine}(\mathbf {u}{}^G;\mathbf {f}{}^F)\), \(\mathbf {u}{}^F\) is defined in the Fine model and we have \(\varvec{\lambda }{}^F=\mathcal {S}{}^F(\mathbf {u}{}^G_\Gamma ;\mathbf {f}{}^F)\).

\([\varvec{\lambda }{}^A]=\mathtt {SolveAux}(\mathbf {u}{}^G;\mathbf {f}{}^A)\), which in corresponds to \(\varvec{\lambda }{}^A=\mathcal {S}{}^A(\mathbf {u}{}^G_\Gamma ;\mathbf {f}{}^A)\). When authorized by the software, it can be replaced by the postprocessing of the stress (19).
The Fine and Auxiliary solves are in general gathered in one line because the computations can be run in parallel.
In order to keep notations simple, we assume that the coarse and fine meshes are conforming at the interface. In most cases the zone of interest is deduced from an initial coarse computation and it is defined as a subset of coarse elements. Then the interface \(\Gamma \) lies on faces of coarse elements and it defines the boundary of the Fine domain. In that context, even if the Fine discretization is chosen to be finer than the Global one on \(\Gamma \), a simple Globalmaster–Fineslave strategy may give satisfying results: a transfer matrix \(\mathbf {T}\) is computed (for instance using interpolation or Mortar techniques) such that the interface conditions can be written as:
so that the mechanical work is preserved. The algorithms presented below are unchanged, one just need to consider the Global interface as the master interface. Note that “Overlapping version” section presents a technique to handle noncoincident interfaces.
Stationary iterations
The global/local coupling iterations of Eq. (7) can formally be written as:
One recognizes fixed point iterations. The convergence is controlled by the contraction property of the operator \(I\mathcal {S}{}^R\circ \mathcal {S}{}^{G^{1}}=(\mathcal {S}{}^A\mathcal {S}{}^F)\circ (\mathcal {S}{}^A+\mathcal {S}{}^C)^{1}\).
Remark 4
In the linear case, the convergence condition can be written in term of spectral radius as:
where we use the ordering of the quadratic form associated with the (symmetric nonnegative) Schur complements. A trivial sufficient condition for the operator to be a contraction is thus \(\mathbf {S}{}^A\geqslant \mathbf {S}{}^F\). Mechanically speaking this means that the Auxiliary model shall be stiffer than the Fine one; this is usually the case when the Fine model has a refined mesh or holes, but this might not be the case if Fine reinforcements are omitted in the Auxiliary model (as in loaded rubbers). Moreover, we can expect the Auxiliary model to be a good approximation of the Fine model leading to \((\mathbf {S}{}^A\simeq \mathbf {S}{}^F)\) and fast convergence. \(\square \)
A classical tweak for fixed point iterations is to use relaxation. This enables to grant contraction property or to improve the convergence rate. In that case, the iteration writes:
For linear problems, it is well known that convergence is ensured for \(0<\omega < 2/\rho \left( {\mathbf {S}{}^G}^{1}\mathbf {S}{}^R\right) \) and the optimal value is \(\omega =2/(\mu _{\min }+\mu _{\max })\) where the \(\mu \)s’ stand for the minimal and maximal eigenvalues of \({\mathbf {S}{}^G}^{1}\mathbf {S}{}^R\). In [10] a technique to approximate the optimal relaxation was proposed. Anyhow dynamic tweaking of the coefficient using Aitken’s formula or, in the linear case, using Krylov solvers seem to be more pragmatic choices (see below). An alternative, when the area to be reanalyzed is known a priori, would be to choose the global model in order for the Auxiliary model to be slightly stiffer than the Fine one; this can be done by adding matter (filling holes), by using a homogenized model with Voigt bound (in case of heterogeneity in the Fine model) and by using a coarse mesh.
The results of the existence of sufficient and of optimal relaxations can be extended to the case of a monotone problem. Indeed in that case, the method can be interpreted as an operator splitting technique [35] on the condensed problem which inherits the useful properties of the original system (in particular monotonicity and coercivity). Reader may refer to [36] for detailed proof with weak assumptions.
In practice, it is convenient to have \(\omega \) adapted at each step. A good heuristic for the sequence \((\omega _n)\) is provided by Aitken’s \(\Delta ^2\). It was first tried in the global/local framework in [37]. The strategy is summedup in Algorithm 1.
QuasiNewton’s approaches for linear Global model
The system to solve associated with the fixed point iterations (26) writes:
which mechanically means that we seek the surface traction to impose inside the Global coarse model (ie the stress discontinuity) such that the Fine model with Dirichlet conditions issued from the Global solve is in balance with the Complement zone.
Applying a Newton iteration to system (29) leads to the sequence:
which was investigated in nonlinear relocalization techniques [14,15,16, 38], but which is in general not possible in the noninvasive framework. Anyhow, in the case of a linear Global model, it is possible to derive a quasiNewton approach.
If the Global problem is linear then the differential of the Global problem is constant and it is equal to the Schur complement \(\mathbf {D}\mathcal {S}{}^G=\mathbf {S}{}^G\). Regarding the nonlinear part, we have:
of course \(\mathbf {X}\) is not computable in a noninvasive manner, but a low rank approximation is possible using quasiNewton formulas. In particular, SR1 formula was tried with success in [1]. In practice, line search is not applied which makes the low rank update slightly lighter than usual, we note \(\varvec{\delta }^u_j=(\mathbf {u}_j\mathbf {u}_{j1})_{\Gamma }\) and \(\varvec{\delta }^p_j=(\mathbf {p}_j\mathbf {p}_{j1})\), the increment of the interface quantities:
with, for \(i>0\), \(\mathbf {R}_i = [\mathbf {r}_1 \ldots \mathbf {r}_{i}]\) and \(\varvec{\Delta }_i={\text {diag}}(\mathbf {r}_i^T \varvec{\delta }^u_i)\). Sherman–Morrison formula leads to:
It makes sense to first evaluate \(\mathbf {S}{}^{G^{1}}\mathbf {r}_i\) then apply corrections. For efficiency reasons, we also store the matrix \(\mathbf {W}_i:=\mathbf {S}{}^{G^{1}}\mathbf {R}_i\). The factorization of Matrix \(\left( \varvec{\Delta }_i+\mathbf {R}_i^T\mathbf {S}{}^{G^{1}}\mathbf {R}_i\right) \) is reused from one iteration to another, only one row and column must be computed. The method is recapitulated in Algorithm 2. Note that this algorithm is written in a way which makes no use of the linearity of the Global problem, so that it will be also tested in the full nonlinear case.
Conjugate gradient
Full linear case
This case occurs when all models are linear. Noninvasive global/local coupling can still be of interest in order to introduce complex local heterogeneities, stochastic behaviors or complex geometries in the Fine model.
For linear problems, it is rather classical to use Krylov accelerators on stationary iterations. In our case, the problem to solve (29) is governed by the operator \(\mathbf {S}{}^R{\mathbf {S}{}^G}^{1}\) which is symmetric in the \({\mathbf {S}{}^G}^{1}\) innerproduct. We then can derive a rightpreconditioned conjugate gradient. The algorithm being not so standard, it is given in Algorithm 3.
Beside the improved convergence compared to stationary iterations, using conjugate gradient allows an unconditional convergence (without necessity for the Auxiliary model to be sufficiently stiff).
Nonlinear case
Conjugate gradient can be extended to nonlinear cases using two ingredients:

A line search algorithm to optimize the length of the steps. For a given search direction \(\underline{\mathbf {p}}\), one tries to find the optimal length \(\alpha \) in term of the minimization of some norm of the residual. This can be done in a noninvasive manner by a sampling technique with several lengths \((\alpha _i)\) being tested in parallel. Classically these samples are used to interpolate the objective function and decide the final \(\alpha \). Because of the cost of the estimation of one configuration (one global solve followed by one local solve), we prefer to use directly the best sample already computed (except if the interpolated minimal let us expect a significantly better configuration).

A “conjugation” technique for the new search direction \(\underline{\mathbf {p}}_{j+1} = \mathbf {r}_{j+1} + \beta _j \underline{\mathbf {p}}_{j}\) given by a heuristic (using the notations of Algorithm 4) like:
$$\begin{aligned} \begin{aligned}&\text {FletcherReeves:}\, \beta _j=\frac{\mathbf {r}_{j+1}^T \mathbf {r}_{j+1}}{\mathbf {r}_{j}^T\mathbf {r}_{j} }&\text {PolacRibi}\grave{\hbox {e}}\text {re:} \ \beta _j=\frac{\mathbf {r}_{j+1}^T (\mathbf {r}_{j+1}\mathbf {r}_{j})}{\mathbf {r}_{j}^T\mathbf {r}_{j} } \\&\text {DaiYuan:}\ \beta _j=\frac{\mathbf {r}_{j+1}^T \mathbf {r}_{j+1}}{\underline{\mathbf {p}}_{j}^T(\mathbf {r}_{j+1}\mathbf {r}_{j}) }&\text {HestenesStiefel:}\ \beta _j=\frac{\mathbf {r}_{j+1}^T (\mathbf {r}_{j+1}\mathbf {r}_{j})}{\underline{\mathbf {p}}_{j}^T(\mathbf {r}_{j+1}\mathbf {r}_{j}) } \end{aligned} \end{aligned}$$(34)Moreover it is often chosen to avoid negative steps by using \(\beta _j \leftarrow \max (0,\beta _j)\). The reader may refer to [39] and associated bibliography for more details. In our examples, the Polac–Ribière formula appeared to be more stable.
Noninvasive implementation
The algorithms above are noninvasive in the sense that they can be implemented as a script driving commercial software with classical input and output instructions. The global/local process is represented in a flow chart on Fig. 2, where the different fields exchanged through the interface \(\Gamma \) are highlighted. In the one hand, the Global and Auxiliary models can be solved by a commercial software, the displacements and tractions are extracted as standard outputs by the script. In the other hand, Fine models can be solved by a research code or by commercial software. The script only needs to be able to transfer interface fields from one numbering to the other (and to compute and apply Matrix \(\mathbf {T}\) in case of nonmatching interface discretizations), and to apply basic algebraic operations in case of acceleration.
Numerical illustration
The method is illustrated on an academic 2D test case modeling a high pressure turbine blade of a plane engine (see Fig. 3) with an approximate size of 10 cm. The Reference model possesses local perforations with adapted mesh which are not present in the Global model. Note that, the Fine model is naturally more flexible than the Auxiliary model (because of the holes and of the refined mesh). The mechanical behavior is either isotropic linear elastic (Young’s modulus E and Poisson’s coefficient \(\nu \) are given in Table 1) or elastoviscoplastic of the form of [40] modeling a realistic IN100 material at hot temperature (\(\simeq 800\,^\circ \)C, parameters are reported in Table 1 using the notations of [40]).
The Fine model is granted an elastoviscoplastic behavior. We consider two configurations: in the first case all models are linear elastic, in the second case all models are elastoviscoplastic. In both cases, the structure is loaded by the body force f which models centrifugal forces due the rotation of the blade, and the traction g on the leading edge which represents the pressure brought by the air flux on the blade.
Linear models
We study the various acceleration strategies presented above in the case where all models are linear elastic. Figure 4a presents the evolution of the Euclidean norm of the residual on the interface. As expected conjugate gradient is faster than other acceleration techniques. Figure 4b presents the evolution of the error measured by the Mises stress on the most loaded element with respect to the Reference model (which should not be available in production cases). We observe an important practical difficulty: Abaqus’ truncation of Gauss point data makes it impossible to observe convergence beyond a relative precision of \(10^{6}\). This problem would appear much later on the residual which only involves nodal computations (which can be manipulated in double precision). Note that this problem was encountered in other studies based on Abaqus software [37] but not in studies which used Code_Aster for example [41].
Table 2 compares the duration of the computation. In that simple case all accelerations have close performance, but we observe that CG is faster than Aitken which is faster than SR1 which is 25% faster than Stationary iteration. The reported CPU time values concern the complete analysis reaching a global/local relative residual norm below \(10^{5}\). The main contribution to the CPU time comes from the finite element solves (proportional to the number of iterations). The remaining difference for a identical number of iterations is due to the complexity of the algorithm (typically SR1 and Aitken converged in the same number of iterations but SR1 involved more complex computations).
Nonlinear models
In that case, all models are granted the same nonlinear elastoviscoplastic behavior. As a consequence, plasticity may spread in the Complement zone. In that case the Fine and Auxiliary models only differ by their topology and their mesh.
Before comparing the acceleration techniques, we specifically study the choice of the parameters of the nonlinear conjugate gradient. Regarding the line search, the stationary iteration corresponds to \(\alpha =1\) (at least at the first iteration), and it behaves rather well (thanks to \(\mathbf {S}{}^A\simeq \mathbf {S}{}^F\)), so it is not absurd to take samples near 1. Of course the size of the sampled interval is also a question of experience for a given class of problems. In the studied case, prior experiments showed that, the optimal linesearch always belonged to the interval [.8, 1.4], we thus use either 4 sampling points \(\{.8,1.,1.2,1.4\}\), 9 sampling points \(\{.8,1.,1.1,1.15,1.2,1.25,1.3,1.35,1.4\}\) or 13 sampling points \(\{.8,1.,1.1,1.15,1.2, 1.22, 1.24, 1.25, 1.26, 1.28,1.3,1.35,1.4\}\).
Figure 5 presents the performance of conjugate gradient for various conjugation techniques and various samplings for the line search. We observe that, in that case, the Polac–Ribière conjugation gives best results, and testing 9 lengths seems to be significantly superior to only testing 4, whereas using a finer sampling discretization is not making improvements in this case. We recall that the line search is conducted in parallel so that oversampling does not take more wallclock time, it only “wastes” machine time and software licenses.
Figure 6 illustrates more precisely the effects of the sampling of the line search: the sampled values of the residual to be minimized are plotted for various lengths, for the initial and the third iterations. We observe that when using 9 or 13 samples, there is always a sample close to the interpolated optimal. One unfortunate effect of nonlinearity is that the finest sampling may lead conjugate gradient to follow a suboptimal solution path: in our case the 9sample lead almost always to a better residual than the 14sample. Another regrettable effect of nonlinearity is that the search direction may lead to an increase of the residual (see the 13sample case after the 8th iteration on Fig. 5b); in that case safeguards should be implemented and the algorithm should switch back to stationary iterations (without conjugation).
We compare, in Fig. 7, the conjugate gradient (in its best configuration, that is to say Polac–Ribière conjugation and 9 sampling points) to other acceleration techniques. We observe that nonlinear conjugate gradient also behaves better than the other techniques. Note that since the sampling of linesearch can be conducted in parallel, one iteration of CG is equivalent (in CPU time) to one iteration of the very cheap Aitken’s method.
Overlapping version
In previous sections, we had assumed that the interface was described as the boundary of elements for all models. In practice this hypothesis is not so restrictive because most often the zone of interest is detected after an initial computation on the coarse Global model, and it is constituted as a set of coarse elements satisfying a certain criterion. Even after remeshing, the boundary of the Fine description of the zone of interest matches a set of coarse faces (edges in 2D). Then a “simple” transfer matrix \(\mathbf {T}\) can be sufficient to communicate between models on the interface (25). In particular, the easy choice of \(\mathbf {T}\) being the interpolation matrix of the coarse kinematics in the fine kinematics can be implemented in most software. More evolved choices like mortar connections can also be employed in certain software [12].
We propose an alternative strategy which makes use of the possibility to have the models overlap. In that case, there is no restriction on the definition of the meshes. This idea can directly be connected to overlapping optimized Schwarz methods, yet we propose a mechanical interpretation of it.
Note that the use of the overlap can be advantageous in the situations where edge effects may affect the Fine model, even if meshes are conforming at the interfaces [34].
Handling of incompatible patches
The starting point is the observation that the method can be formulated as the search for \(\mathbf {p}\) which is the stress discontinuity on the Global model between the Complement zone and the Auxiliary description of the zone of interest. This discontinuity must be such that the Complement zone is in equilibrium with the Fine description of zone of interest loaded with Dirichlet conditions (29).
Since \(\mathbf {p}\) is a discontinuity, in order it to be well described in the coarse finite element model, it must be supported by the boundary of coarse elements. But there is no need for the support of \(\mathbf {p}\) to match the boundary of the zone of interest.
We thus propose to follow Fig. 8. The Fine subdomain \(\Omega {}^F\) is positioned where needed in the zone of interest, its mesh is independent from the coarse mesh. We note \(\Gamma {}^F=\partial \Omega {}^F\) the boundary of the Fine subdomain. The Auxiliary subdomain is the largest set of coarse elements fully contained in the zone of interest. We note \(\Gamma {}^A=\partial \Omega {}^A\) the boundary of the Auxiliary zone. The two interfaces \(\Gamma {}^F\) and \(\Gamma {}^A\) thus do not coincide. \(\Omega {}^C\) is defined as the Complement to \(\Omega {}^A\) in the Global model \(\Omega {}^C=\Omega {}^G\setminus \Omega {}^A\).
The Fine and Global models’ displacements are connected on \(\Gamma {}^F\); \(\mathbf {p}\) is applied to the Global model on \(\Gamma {}^A\) with the aim to reach balance between the nodal reaction from Complement model and the normal from the Fine model projected on \(\Gamma {}^A\). In the end, the coupled solution is \(u{}^F\) in \(\Omega {}^A\) and \(u{}^G\) in \(\Omega {}^G\setminus \Omega {}^F\). In the overlap, also called buffer zone \(\Omega {}^B=\Omega {}^F\setminus \Omega {}^A\), the Complement and the Fine model coexist.
Algorithm 5 gives the basic stationary iteration in the presence of overlap, all acceleration techniques could be considered. In order to distinguish between the interfaces, the Auxiliary and the Fine problems are written on separate lines even if they can be solved in parallel.
The main practical difficulty of this algorithm is the computation of the Fine reaction on \(\Gamma {}^A\) with \(\Omega {}^F\) not exactly represented on the coarse grid. This computation mixes the Fine stress \(\sigma _h{}^F\) and the coarse shape functions \(\phi {}^G_i\). Even if complex, this computation is feasible in certain software. Anyhow in the nonlinear case, \(\sigma _h{}^F\) is only known at Fine Gauss points and the integral can only be approximated.
Note that the conceptual difficulty caused by the coexistence of two models in the overlap is common to other coupling methods [42]. It seems wise not to introduce strong dissimilarities between the Fine and the Global models in the buffer zone, so that only the nonconformity of meshes lead to (hopefully only slightly) different solutions in it.
Usually enlarging the overlap leads to improved converge rate in Schwarz methods. This is not as clear in our method where the Fine and Global models may differ in the buffer zone. Since a good convergence rate is already ensured by the good approximation of the Fine model provided by the Auxiliary model, and because of the ambiguity of the model in the overlap, we see no interest in considering large overlaps (unless it also plays a role from a mechanical point of view, see the discussion on the decay of edge effects below).
Illustration of the coupling with overlap
For now, the version with overlap has only been implemented in the context of the coupling between a Global heterogeneous plate model and a Fine 3D model, for the study of bolted assemblies of thin structures [43], in Code_aster.
Despite the definition of smart transfer techniques between the 3D and the plate models [34], it appeared that 3D edge effects were impossible to avoid completely. Figure 9 shows the variation of the peeling stress in the Fine 3D model direction orthogonal to the interface for various lifting of the plate displacement (named “Lagrangian” and “warping” in the figure). We see that edge effects are important on the boundary of the Fine domain (\(\Gamma ^F\)) but they fade quickly so that the boundary of the Auxiliary domain (\(\Gamma ^A\)) can be positioned not too far inside the zone of interest. The buffer zone \(\Omega ^B\) between \(\Gamma ^F\) and \(\Gamma ^A\) was granted a width of two times the thickness, due to the exponential decay property of edge effects in plates. Note that in that case, the overlap was useful even with matching meshes (plate and 3D nodes were coincident on \(\Gamma {}^A\) and \(\Gamma {}^F\)).
As said earlier, in the case of nonmatching meshes, the difficulty for the coupling with overlap is the computation of the fine reaction on \(\Gamma ^A\) written \(\varvec{\lambda }{}^F_{j,i}\) in Algorithm 5. In [43], it was proposed to extract a band of Auxiliary elements connected to \(\Gamma ^A\) and project on it the Fine stress (defined at the Gauss point of the Fine mesh). This was implemented in Code_Aster using existing routines (PROJ_CHAMP() with keyword ECLA_PG).
Figure 10, from [43], presents the converged solution of the simple application of an isotropic plate in flexion where the Global model is a solid plate with unstructured mesh with triangular elements and the 3D Fine model bears a hole and is meshed with structured hexahedral elements. For now the code for the nonconforming coupling with overlap is just a proof of concept. During our limited number of experiments, the convergence rate was quite comparable with the conforming case. A more extensive numerical campaign should be conducted with varying position for the interfaces and different element sizes (so that situations where the transfers behave poorly may be encountered). Moreover, the method shall also be assessed in term of discretization error in the spirit of [41].
Conclusion
The global/local noninvasive coupling technique is a convenient way to enrich a global coarse model, handled by a commercial software, with local features, handled by the most adapted software. In this paper we proposed to interpret the method as an alternate optimized nonoverlapping Schwarz domain decomposition method. In this framework the coarse representation of the zone of interest is a clever way to build an approximation of the DirichlettoNeumann operator of the Fine model, which includes the effects of the imposed load. Belonging to the Schwarz family of domain decomposition method allows to benefit many theoretical results and practical shrewdness. We then derive a conjugate gradient solver in the linear and nonlinear cases, in that later case the line search is realized by a sampling which can be conducted in parallel in order not to penalize the wallclock time. Finally we show that an overlapping version can also be applied which enables to connect nonmatching meshes.
References
 1.
Gendre L, Allix O, Gosselet P, Comte F. Nonintrusive and exact global/local techniques for structural problems with local plasticity. Comput Mech. 2009;44(2):233–45.
 2.
Kelley F. Mesh requirements for the analysis of a stress concentration by the specified boundary displacement method. In: Proceedings of the second international computers in engineering conference. New York City: ASME; 1982. p. 39–42.
 3.
Ransom JB, McCleary SL, Aminpour MA, Knight NF Jr. Computational methods for global/local analysis. NASA STI/Recon Technical Report N, vol. 92. 1992. p. 33104.
 4.
Cormier NG, Smallwood BS, Sinclair GB, Meda G. Aggressive submodelling of stress concentrations. Int J Numer Methods Eng. 1999;46(6):889–909. 10.1002/(SICI)10970207(19991030)46:6<889::AIDNME699>3.0.CO;2F.
 5.
JaraAlmonte CC, Knight CE. The specified boundary stiffness/force SBSF method for finite element subregion analysis. Int J Numer Methods Eng. 1988;26(7):1567–78.
 6.
Whitcomb JD. Iterative global/local finite element analysis. Comput Struct. 1991;40(4):1027–31.
 7.
Whitcomb JD, Woo K. Application of iterative global/local finiteelement analysis. part 1: linear analysis. Commun Numer Methods Eng. 1993;9:745.
 8.
Hecht F, Lozinski A, Pironneau O. Numerical zoom and the Schwarz algorithm. In: Proceedings of the 18th conference on domain decomposition methods. 2009.
 9.
Duval M, Passieux JC, Salaün M, Guinard S. Nonintrusive coupling: recent advances and scalable nonlinear domain decomposition. Arch Comput Methods Eng. 2016;23(1):17–38. https://doi.org/10.1007/s118310149132x
 10.
Chevreuil M, Nouy A, Safatly E. A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties. Comput Methods Appl Mech Eng. 2013;255:255–74. https://doi.org/10.1016/j.cma.2012.12.003.
 11.
Guguin G, Allix O, Gosselet P, Guinard S. On the computation of plate assemblies using realistic 3D joint model: a nonintrusive approach. Adv Model Simul Eng Sci. 2016;3:16. https://doi.org/10.1186/s4032301600695.
 12.
Duval M, Passieux JC, Salaün M, Guinard S. Local/global nonintrusive parallel coupling for large scale mechanical analysis. In: 11th world congress on computational mechanics—5th European conference on computational mechanics. IACMECCOMAS. 2014.
 13.
Keyes DE. Aerodynamic applications of NewtonKrylovSchwarz solvers. In: Fourteenth international conference on numerical methods in fluid dynamics. Berlin: Springer; 1995. p. 1–20.
 14.
Cresta P, Allix O, Rey C, Guinard S. Nonlinear localization strategies for domain decomposition methods: application to postbuckling analyses. Comput Methods Appl Mech Eng. 2007;196(8):1436–46.
 15.
Hinojosa J, Allix O, Guidault PA, Cresta P. Domain decomposition methods with nonlinear localization for the buckling and postbuckling analyses of large structures. Adv Eng Softw. 2014;70:13–24.
 16.
Negrello C, Gosselet P, Rey C, Pebrel J. Substructured formulations of nonlinear structure problems—influence of the interface condition. Int J Numer Methods Eng. 2016;107(13):1083–105. https://doi.org/10.1002/nme.5195.
 17.
Bettinotti O, Allix O, Malherbe B. A coupling strategy for adaptive local refinement in space and time with a fixed global model in explicit dynamics. Comput Mech. 2014;53(4):561–74. https://doi.org/10.1007/s0046601309179
 18.
Bettinotti O, Allix O, Perego U, Oncea V, Malherbe B. A fast weakly intrusive multiscale method in explicit dynamics. Int J Num Methods Eng. 2014;100(8):577–95. https://doi.org/10.1002/nme.4750.
 19.
Bettinotti O, Allix O, Perego U, Oncea V, Malherbe B. Simulation of delamination under impact using a global local method in explicit dynamics. Finite Elem Anal Des. 2017;125:1–13. https://doi.org/10.1016/j.finel.2016.11.002.
 20.
Plews J, Duarte C, Eason T. An improved nonintrusive global–local approach for sharp thermal gradients in a standard FEA platform. Int J Numer Methods Eng. 2011;91(4):426–49.
 21.
Kim J, Duarte CA. A new generalized finite element method for twoscale simulations of propagating cohesive fractures in 3D. Int J Numer Methods Eng. 2015;103(13):1139–72. https://doi.org/10.1002/nme.4954.
 22.
Passieux JC, Réthoré J, Gravouil A, Baietto MC. Local/global nonintrusive crack propagation simulation using a multigrid XFEM solver. Comput Mech. 2013;52(6):1391–393. https://doi.org/10.1007/s0046601308823
 23.
Ciarlet PG. Linear and nonlinear functional analysis with applications. Philadelphia: SIAM; 2013.
 24.
Showalter RE. Hilbert space methods for partial differential equations. Monographs and studies in mathematics. London: Pitman; 1977.
 25.
Showalter RE. Monotone operators in Banach space and nonlinear partial differential equations. Mathematical surveys and monographs, vol. 49. Providence: American mathematical society; 1997.
 26.
Feistauer M, Zenisek A. Finite element solution of nonlinear elliptic problems. Numerische Mathematik. 1987;50:451–75.
 27.
Ladevèze P. Sur une famille d’algorithmes en mécanique des structures. Comptes Rendus Académie des Sciences  Mécanique, Paris. 1985;300(2):41–5.
 28.
Ladevèze P. Nonlinear computational structural mechanics—new approaches and nonincremental methods of calculation. Berlin: Springer; 1999.
 29.
Hauer D. The pDirichlettoNeumann operator with applications to elliptic and parabolic problems. J Differ Equ. 2015;259(8):3615–55.
 30.
Gander MJ, Halpern L. Méthodes de décomposition de domaines. Mathématiques pour l’ingénieur. Techniques de l’Ingénieur, SaintDenis, France. 2012.
 31.
Gendre L, Allix O, Gosselet P. A twoscale approximation of the Schur complement and its use for nonintrusive coupling. Int J Numer Methods Eng. 2011;87(9):889–905.
 32.
Badea L. On the Schwarz alternating method with more than two subdomains for nonlinear monotone problems. SIAM J Numer Anal. 1991;28(1):179–204. https://doi.org/10.1137/0728010.
 33.
Ladevèze P. Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation. Mechanical engineering series. New York: Springer; 1999. https://doi.org/10.1007/9781461214328.
 34.
Guguin G, Allix O, Gosselet P, Guinard S. Nonintrusive coupling of 3D and 2D laminated composite models based on finite element 3D recovery. Int J Numer Methods Eng. 2014;98(5):324–43.
 35.
Ryu EK, Boyd S. Primer on monotone operator methods. 2015.
 36.
Nouy A, Pled F. A multiscale method for semilinear elliptic equations with localized uncertainties and nonlinearities. 2017. arXiv:1704.05331v2.
 37.
Liu YJ, Sun Q, Fan XL. A nonintrusive global/local algorithm with nonmatching interface: derivation and numerical validation. Comput Methods Appl Mech Eng. 2014;277:81–103.
 38.
Pebrel J, Gosselet P, Rey C. Une approche par patchs pour les non linéarités localisées en calcul de structures. 18ème Congrès Français de Mécanique (Grenoble 2007). 2007.
 39.
Glowinski R. Variational methods for the numerical solution of nonlinear elliptic problems. Philadelphia: Society for Industrial and Applied Mathematics; 2015. https://doi.org/10.1137/1.9781611973785.
 40.
Longuet A, Burteau A, Comte F, CrouchezPilot A. Incremental lifing method applied to high temperature aeronautical component. CSMA 2013. 2013.
 41.
Duval M, Losinski A, Passieux JC, Salaün M. Residual error based adaptive mesh refinement with the nonintrusive patch algorithm. Comput Methods Appl Mech Eng. 2018;329:118–43.
 42.
Ben Dhia H. Multiscale mechanical problems: the Arlequin method. Comptes Rendus de l’Academie des Sciences Series IIB Mechanics Physics Astronomy. 1998;326(12):899–904.
 43.
Guguin G. Stratégie nonintrusive de couplage plaque/3d pour l’application aux plaques composites stratifiées. Ph.D., École Normale Supérieure de Cachan
 44.
Chaboche JL, Lemaitre J, Benallal A, Desmorat R. Mécanique des Matériaux Solides. Paris: Dunod; 2009.
Author's contributions
This work on the method initiated by OA is in the line of previous contributions by OA and PG. PG proposed the Schwarz’ interpretation. All authors contributed to the study of the acceleration methods. MB implemented the method and conducted the numerical experiments in the absence of overlap. Numerical experiments in the presence of overlap are extracted from the Ph.D. of GG. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
Not applicable.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Funding
This work was partially funded by the French National Research Agency as part of project ICARE (ANR12MONU000204).
Author information
Affiliations
Corresponding author
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Gosselet, P., Blanchard, M., Allix, O. et al. Noninvasive global–local coupling as a Schwarz domain decomposition method: acceleration and generalization. Adv. Model. and Simul. in Eng. Sci. 5, 4 (2018). https://doi.org/10.1186/s4032301800974
Received:
Accepted:
Published:
Keywords
 Noninvasive coupling
 Schwarz domain decomposition
 Nonconforming patch
 Krylov acceleration
Mathematics Subject Classification
 Primary 65N55
 65N85