# Substructured formulations of nonlinear structure problems - influence of the interface condition

Camille Negrello, Pierre Gosselet, Christian Rey, Julien Pebrel

► To cite this version:

Camille Negrello, Pierre Gosselet, Christian Rey, Julien Pebrel. Substructured formulations of non-linear structure problems - influence of the interface condition. International Journal for Numerical Methods in Engineering, 2016, 107 (13), pp.1083-1105. 10.1002/nme.5195 . hal-01332214

HAL Id: hal-01332214

https://hal.science/hal-01332214

Submitted on 15 Jun 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>
# Substructured formulations of nonlinear structure problems – Influence of the interface condition

Camille Negrello¹, Pierre Gosselet¹, Christian Rey¹,² and Julien Pebrel¹

(1) LMT Cachan, ENS Cachan/CNRS/Univ. Paris Saclay, 61 Avenue du Président Wilson, 94235 Cachan France.

(2) Safran Tech, rue des jeunes bois, Chateaufort, CS 80112, 78772 Magny-les-Hameaux France.

June 15, 2016

## Abstract

We investigate the use of non-overlapping domain decomposition (DD) methods for nonlinear structure problems. The classic techniques would combine a global Newton solver with a linear DD solver for the tangent systems. We propose a framework where we can swap Newton and DD, so that we solve independent nonlinear problems for each substructure and linear condensed interface problems. The objective is to decrease the number of communications between subdomains and to improve parallelism. Depending on the interface condition, we derive several formulations which are not equivalent, contrarily to the linear case. Primal, dual and mixed variants are described and assessed on a simple plasticity problem.

**Keywords:** Domain decomposition; nonlinear mechanics; Newton solver; Krylov solver; parallel processing

## 1 Introduction

In order to solve large nonlinear structure problems, an efficient strategy is to combine a Newton-based solver which leads to a sequence of linear systems, and a domain decomposition approach to solve the tangent systems. Indeed such a strategy combines well-known and robust methods for which many refinement are available: Newton can be tangent/constant/secant/modified/arc-length [24], and domain decomposition solvers [33, 13, 32, 10, 17] can be equipped with preconditioners and coarse problems which make them reliable and scalable [26, 27, 44, 19]. Moreover, the computation of stiffness matrices being done in parallel independently for each subdomain, and information being reusable from one system to another [37, 39, 18], the overall performance is very satisfying in general [11, 3].

But there are cases where such an approach is not as pertinent as expected. For instance, when dealing with strong localized nonlinearity, many (global) Newton iterations are required, whereas most of the structure undergoes a linear evolution. In this situation, it would be interesting to differentiate between the linear or nonlinear nature of each subdomain's behavior, in order to decrease the number of global iterations and communications. The possibility to conduct local nonlinear computations has been investigated for a long time in the Schwarz framework (with or without overlap between subdomains) [1, 8, 30, 4, 31, 22]. More recently the possibility to define nonlinear versions of the Schur complement methods (often called nonlinear relocalization techniques) was studied: BDD in [6, 2], FETI in [35], FETI-DP in [25], BDDC in [21, 25].

The aim of this paper is to give a formal framework to develop the nonlinear versions of the well-known linear solver FETI [13], BDD [33, 32] and FETI2LM[41]. They rely on the concept of nonlinear condensation which was exposed in conferences [36] and seminars [16].

In section 2 and 3, we present the nonlinear system in its monolithic and substructured form. In section 4, we introduce the concept of nonlinear condensation in its primal, dual and mixed versions. In section 5, we show that applying a Newton method to these condensed formulations leads, at each outer iteration, to the parallel solution of local nonlinear systems with, depending on the formulation, Dirichlet, Neumann or Robin boundary conditions, and to the linear interface systems which can be solved by classical BDD, FETI or FETI2LM. Iterative solvers being involved inside the outer loop, we interpret the method as an inexact Newton solver [7] which leads us, in section 6, to tune the convergence thresholds depending on the current residual in order to avoid oversolving. First assessments are given in section 7 on a simple problem, yet representative of the nonlinearity encountered in industrial problems; a discussion ends the paper.

## 2 Reference problem

We consider the classic problem of the evolution, under the small perturbation hypothesis, of a structure occupying the domain $\Omega$, submitted to body forces $f$, to traction forces $F$ on the part $\partial_f\Omega$ of its boundary and to given displacements
$u_g$ on the complementary part $\partial_u \Omega \neq \emptyset$. Note that the small perturbation hypothesis is crucial for the handling of rigid body motions in the dual approach, but it can be relaxed in other cases.

Let $\varepsilon(u)$ be the symmetric part of the gradient of displacement and $\sigma$ the Cauchy stress tensor. The problem to be solved can be written as:

at time $t \in [0, T]$, find $u \in \mathcal{U} / \forall v \in \mathcal{U}_0$

$$ \int_{\Omega} \sigma : \varepsilon(v) d\Omega = \int_{\Omega} v^T f d\Omega + \int_{\partial_f \Omega} v^T F dS $$

$$ \sigma = \sigma(\varepsilon(u(\tau)), \tau \in [0, t]) $$

where $\mathcal{U}$ is the space of kinematically admissible fields and $\mathcal{U}_0$ is the associated vector space:

$$ \mathcal{U} = \{u \in H^1(\Omega), u = u_g \text{ on } \partial_u \Omega\} $$

The notation $\sigma(\varepsilon(u(\tau)), \tau \in [0, t])$ means that the stress at one point depends on the whole history of the strain at this point (local nonlinear behavior). This history is most often materialized by internal variables like anelastic strain, hardening or damage.

The problem is discretized in space using the finite element method. The domain $\Omega$ is meshed; let $N$ be the matrix of shape functions, and $u$ be the nodal displacement unknowns such that $u = N u$. The problem is also supposed to be discretized in time, the discrete reference problem to be solved at Step $t_n$ can be written as:

$$ \text{Find } u(t_n) \text{ so that } f_{\text{int}}(u) + f_{\text{ext}} = 0 \quad (1) $$

with

$$ f_{\text{ext}} = \int_{\Omega} N^T f d\Omega + \int_{\partial_f \Omega} N^T F dS $$

$$ v^T f_{\text{int}} = - \int_{\Omega} \sigma_h : \varepsilon(Nv) d\Omega $$

$$ \sigma_h = \sigma_h(\varepsilon(Nu(t_j)), j \le n) $$

$\sigma_h$ is the discrete counterpart of $\sigma$, it depends of the whole discrete history. Integrals are classically computed numerically using Gauss quadrature, so that stress and internal variables are defined at Gauss points. Throughout the rest of the document, the time step will not be mentioned. In order to shorten expressions, Dirichlet boundary conditions $u_g$ are implicitly taken into account within $f_{\text{int}}$ and $f_{\text{ext}}$.

Note that in the linear case we have: $f_{\text{int}}(u) = -K u$ where $K$ is the stiffness matrix.

The classical solution strategy to Problem (1) is to use a Newton-Raphson algorithm to linearize the problem, and then solve a sequence of tangent systems.

## 3 Substructured formulation

We consider the conforming partition of $\Omega$ into $N$ non-overlapping subdomains $\Omega^{(s)}$, so that each element exactly belongs to one subdomain. Superscript (s) will refer to data attached to domain $\Omega^{(s)}$. In order to ease the treatment of groups of subdomains, we define the following block notations:

$$ x^{\Diamond} = \begin{pmatrix} \vdots \\ x^{(s)} \\ \vdots \end{pmatrix}, \quad x^{\Diamond} = (\dots x^{(s)} \dots), \quad x^{\Diamond} = \begin{pmatrix} \vdots & 0 \\ 0 & x^{(s)} \end{pmatrix} $$

Note that in the case of nonlinearly-dependent data, we use the same notation but the dependence is implicitly local:

$$ f_{\text{int}}^{\Diamond}(u^{\Diamond}) = \begin{pmatrix} \vdots \\ f_{\text{int}}^{(s)}(u^{(s)}) \\ \vdots \end{pmatrix} $$

The nodes on the interface between two subdomains play a specific role which we need to highlight. We note $\Gamma^{(i,j)}$ the set of nodes shared by subdomains $\Omega^{(i)}$ and $\Omega^{(j)}$, and $\Gamma^{(s)} = \cup_j \Gamma^{(s,j)}$ the interface nodes of subdomain $\Omega^{(s)}$. We use the subscript $b$ for nodes belonging to the interface $\Gamma^{(s)}$ and the subscript $i$ for internal nodes. The trace operator $t^{(s)}$ extracts interface nodal values (on $\Gamma^{(s)}$) from subdomain data (in $\Omega^{(s)}$):

$$ \forall s \in [1, N], \quad t^{(s)}u^{(s)} = u_b^{(s)} \quad \text{which also can be written as} \quad t^{\Diamond}u^{\Diamond} = u_b^{\Diamond} $$

$$ \text{assuming adapted ordering, } t^{(s)} = (o_{bi}^{(s)} \quad I_{bb}^{(s)}) $$

See Figure 1(a,b) for an example. We note $\mathbb{R}_b^{\Diamond} = \text{Range}(t^{\Diamond})$ the space of local interface vectors to which any vector $x_b^{\Diamond}$ belongs.
Let $\Gamma_A = \cup \Gamma^{(s)}$ denote the totality of interface nodes, and let $A^{(s)}$ be the canonical operator which injects nodes from $\Gamma^{(s)}$ to $\Gamma_A$. $A^{(s)}$ is a boolean full column-rank matrix. With one local interface node in $\Gamma^{(s)}$ associated with exactly one node in $\Gamma_A$, the global assembling operator $A^{\diamond}$ has no left-kernel. It is then a full row-rank matrix (sometimes called a primal assembling operator). See Figure 1(b,c) for an example. We note $\mathbb{R}^{\Gamma_A} = \text{Range}(A^{\diamond})$ the space of vectors defined on that “primal” interface.

Any matrix $B^{\diamond}$ satisfying $\text{Range}(B^{\diamond T}) = \text{Ker}(A^{\diamond})$ can be used as a dual assembling operator. In Figure 1(b,d) we give the most classical choice where $B^{\diamond}$ is a signed boolean matrix which describes each connection between interface nodes; we note $\Gamma_B$ the set of connections. Note that operator $B^{\diamond}$ needs not to be full row-rank (in the classical case it is not whenever one interface degree of freedom is shared by more than two subdomains). We note $\mathbb{R}^{\Gamma_B} = \text{Range}(B^{\diamond})$ the space of vectors defined on that “dual” interface.

In order to decouple equations between subdomains, we introduce the nodal reaction $\lambda_b^{(s)}$ imposed on subdomain $\Omega^{(s)}$ by its neighbors. The substructured reference problem can be written as:

$$ 
\begin{array}{l}
\text{Find } (u^{\diamond}, \lambda_b^{\diamond}) \text{ such that} \\
\begin{aligned}
& f_{\text{int}}^{\diamond}(u^{\diamond}) + f_{\text{ext}}^{\diamond} + t^{\diamond T} \lambda_b^{\diamond} = 0 \\
& A^{\diamond} \lambda_b^{\diamond} = 0 \\
& B^{\diamond} u_b^{\diamond} = 0
\end{aligned}
\end{array} \qquad (2)
$$

where the first equation expresses the equilibrium of each subdomain submitted to given forces and to unknown interface reactions $\lambda_b^{\diamond}$. The primal assembling operator $A^{\diamond}$ enables us to express the action-reaction principle while the dual assembling operator $B^{\diamond}$ enables us to express the continuity of the displacement field. In addition to Figure 1, note that in the case of two subdomains these operators can be written as $A^{\diamond} \lambda_b^{\diamond} = \lambda_b^{(1)} + \lambda_b^{(2)}$ and $B^{\diamond} u_b^{\diamond} = u_b^{(1)} - u_b^{(2)}$.

Assembling operators satisfy the following relationships:

* Assembling operators are orthogonal in the following sense:

$$ A^{\diamond} B^{\diamond T} = 0 \qquad (3) $$

* Assembling operators generate local interface nodal vectors:

$$ \text{Range}(A^{\diamond T}) \oplus \text{Range}(B^{\diamond T}) = \mathbb{R}_b^{\diamond} \qquad (4) $$

* Any local interface vector is uniquely defined as a combination of a balanced vector $B^{\diamond T} x_B$ and a continuous vector $A^{\diamond T} x_A$:

$$ 
\begin{array}{l}
\forall x_b^{\diamond}, \exists (x_B, x_A) \in \mathbb{R}^{\Gamma_B} \times \mathbb{R}^{\Gamma_A} / x_b^{\diamond} = B^{\diamond T} x_B + A^{\diamond T} x_A \\
\text{indeed} \quad \left\{ \begin{array}{l} x_A = (A^{\diamond} A^{\diamond T})^{-1} A^{\diamond} x_b^{\diamond} \\ x_B = (B^{\diamond} B^{\diamond T})^{+} B^{\diamond} x_b^{\diamond} \end{array} \right. \qquad (5)
\end{array}
$$

The use of the pseudo-inverse $(B^{\diamond} B^{\diamond T})^{+}$ is due to the potential presence of redundancies in the description of the connectivity between subdomains, though it is applied to a vector which belongs to $\text{Range}(B^{\diamond}) = \text{Range}(B^{\diamond} B^{\diamond T})$ so that $x_B$ is well defined, and $B^{\diamond T} x_B$ does not depend on the choice of the pseudo-inverse.

# 4 Nonlinear condensations

The substructured formulation (2) is strictly equivalent to the global formulation (1). We now propose various solution strategies which (under certain assumptions) all converge to the reference solution, though these methods, which also employ Newton-Raphson algorithm, are meant to generate a different sequence of linear systems.

## 4.1 Primal formulation

The primal formulation consists in rewriting system (2) in terms of one unknown interface displacement field $u_A$:

$$ 
\begin{array}{l}
\text{Find } u_A \in \mathbb{R}^{\Gamma_A} \text{ such that } A^{\diamond} \lambda_b^{\diamond} = 0 \\
\text{where } \lambda_b^{\diamond} := -[f_{\text{int}}^{\diamond}(u^{\diamond}) + f_{\text{ext}}^{\diamond}]_b \\
\text{and } u^{\diamond} \text{ solves } \left\{ \begin{array}{l} [f_{\text{int}}^{\diamond}(u^{\diamond}) + f_{\text{ext}}^{\diamond}]_i = 0 \\ t^{\diamond} u^{\diamond} = A^{\diamond T} u_A \end{array} \right. \qquad (6)
\end{array}
$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Local numberings, interface numberings, trace and assembly operators</figcaption></figure>

The last set of equations corresponds to the solution to independent mechanical problems for each subdomain with imposed displacement at the interface (a Dirichlet condition which implies that the continuity is automatically insured: $\mathbf{u}_b^\Diamond = \mathbf{A}^\Diamond T \mathbf{u}_A \Rightarrow \mathbf{B}^\Diamond \mathbf{u}_b^\Diamond = 0$). This displacement has to be found so that the associated reactions $\lambda_b^\Diamond$ are balanced on the interface.

If we assume that the last set of equations has a unique solution for any imposed interface displacement, then we can define an operator $S_{nl}^{(s)}$ so that:

$$ \lambda_b^{(s)} = S_{nl}^{(s)} (\mathbf{A}^{(s)T} \mathbf{u}_A; \mathbf{f}_{ext}^{(s)}) \quad (7) $$

This operator can be viewed as a nonlinear version of the Schur complement; it computes the reaction associated with a given displacement. In the linear case, an explicit expression can be given:

$$ \begin{gathered} S_t^{(s)}(\mathbf{u}_b^{(s)}; \mathbf{f}_{ext}^{(s)}) = S_t^{(s)}\mathbf{u}_b^{(s)} - \mathbf{b}_p^{(s)} \\ \text{with } S_t^{(s)} = \mathbf{K}_{bb}^{(s)} - \mathbf{K}_{bi}^{(s)} \mathbf{K}_{ii}^{(s)-1} \mathbf{K}_{ib}^{(s)} \quad (8) \\ \mathbf{b}_p^{(s)} = \mathbf{f}_{extb}^{(s)} - \mathbf{K}_{bi}^{(s)} \mathbf{K}_{ii}^{(s)-1} \mathbf{f}_{exti}^{(s)} \end{gathered} $$
The primal nonlinear condensed system reads:

$$ \text{Find } \mathbf{u}_A \in \mathbb{R}^{\Gamma_A} \text{ such that} \\ \mathbf{A}^{\Diamond} \mathbf{S}_{nl}^{\Diamond} (\mathbf{A}^{\Diamond T} \mathbf{u}_A; \mathbf{f}_{ext}^{\Diamond}) = 0 \qquad (9) $$

## 4.2 Dual formulation

The dual formulation consists in rewriting system (2) in terms of one unknown interface reaction field:

$$ \text{Find } \lambda_B \in \mathbb{R}^{\Gamma_B} \text{ such that } \mathbf{B}^{\Diamond} \mathbf{t}^{\Diamond} \mathbf{u}^{\Diamond} = 0 \\ \text{where } \mathbf{u}^{\Diamond} \text{ solves } \mathbf{f}_{int}^{\Diamond}(\mathbf{u}^{\Diamond}) + \mathbf{f}_{ext}^{\Diamond} + \mathbf{t}^{\Diamond T} \mathbf{B}^{\Diamond T} \lambda_B = 0 \qquad (10) $$

The last equation corresponds to the solution to independent mechanical problems for each subdomain with imposed traction at the interface (local Neumann problems). This traction has to be found so that displacements are continuous at the interface. Implicitly local reactions have been defined by $\lambda_b^{\Diamond} = \mathbf{B}^{\Diamond T} \lambda_B$, which guarantees the interface equilibrium $\mathbf{A}^{\Diamond} \lambda_b^{\Diamond} = 0$.

In order for problem (10) to be well-posed, one necessary condition is the self-equilibrium of the substructure. Let $\mathbf{R}^{(s)}$ be the basis of kinematically admissible rigid body motions of subdomain $\Omega^{(s)}$ (if they exist, these displacements are exactly represented in the finite element function space), then we have $\mathbf{R}^{(s)T} \mathbf{f}_{int}^{(s)} = 0$ (internal forces do not develop work within rigid body motions). This implies the admissibility condition for $\lambda_B$:

$$ \mathbf{R}^{\Diamond T} (\mathbf{f}_{ext}^{\Diamond} + \mathbf{t}^{\Diamond T} \mathbf{B}^{\Diamond T} \lambda_B) = 0 \qquad (11) $$

If we assume that (10) has a unique solution (up to a rigid body motion $\mathbf{R}^{\Diamond} \alpha^{\Diamond}$, where $\alpha^{\Diamond}$ is the unknown amplitude of rigid body motions) for any imposed interface traction, then we can define an operator $\mathbf{F}_{nl}^{(s)}$ such that:

$$ \mathbf{u}_b^{(s)} = \mathbf{F}_{nl}^{(s)} (\mathbf{B}^{(s)T} \lambda_B; \mathbf{f}_{ext}^{(s)}) + \mathbf{t}^{(s)} \mathbf{R}^{(s)} \alpha^{(s)} \qquad (12) $$

This operator can be viewed as a nonlinear version of the dual Schur complement (as employed in classical FETI methods [12] of which we borrow the notation $\mathbf{F}$); it computes the displacement associated with given reaction. In the linear case, an explicit expression can be given:

$$ \begin{gathered} \mathbf{F}_t^{(s)}(\lambda_b^{(s)}; \mathbf{f}_{ext}^{(s)}) = \mathbf{F}_t^{(s)}\lambda_b^{(s)} + \mathbf{b}_d^{(s)} \\ \text{with } \mathbf{F}_t^{(s)} = \mathbf{t}^{(s)}\mathbf{K}^{(s)} + \mathbf{t}^{(s)T} \qquad (13) \\ \mathbf{b}_d^{(s)} = \mathbf{t}^{(s)}\mathbf{K}^{(s)} + \mathbf{f}_{ext}^{(s)} \end{gathered} $$

We recall the following classical relationships [17]:

$$ \begin{gathered} \mathbf{R}^{(s)} = \text{Ker}(\mathbf{K}^{(s)}) \\ \mathbf{R}_b^{(s)} = \mathbf{t}^{(s)} \mathbf{R}^{(s)} = \text{Ker}(\mathbf{S}_t^{(s)}) \\ \mathbf{R}^{(s)T} \mathbf{f}_{ext}^{(s)} = \mathbf{R}_b^{(s)T} \mathbf{b}_p^{(s)} \end{gathered} \qquad (14) $$

Finally, the dual nonlinear condensed system reads:

$$ \text{Find } \lambda_B \in \mathbb{R}^{\Gamma_B}, \alpha^{\Diamond} \text{ such that} \\ \left\{ \begin{array}{l} \mathbf{B}^{\Diamond} (\mathbf{F}_{nl}^{\Diamond} (\mathbf{B}^{\Diamond T} \lambda_B; \mathbf{f}_{ext}^{\Diamond}) + \mathbf{t}^{\Diamond} \mathbf{R}^{\Diamond} \alpha^{\Diamond}) = 0 \\ \mathbf{R}^{\Diamond T} (\mathbf{f}_{ext}^{\Diamond} + \mathbf{t}^{\Diamond T} \mathbf{B}^{\Diamond T} \lambda_B) = 0 \end{array} \right. \qquad (15) $$

## 4.3 Mixed formulation

The mixed formulation consists in introducing a new interface variable for each subdomain $\mu_b^{(s)}$:

$$ \mu_b^{\Diamond} = \lambda_b^{\Diamond} + \mathbf{Q}_b^{\Diamond} \mathbf{u}_b^{\Diamond} \qquad (16) $$

The symmetric positive definite matrix $\mathbf{Q}_b^{\Diamond}$ is a parameter of the method. It can be interpreted as a stiffness (or an impedance) added to the interface (like in Robin-type boundary conditions).
The properties of the matrix $Q_b^\diamond$ imply a partitioning of the space of subdomain interface vectors:

$$ \mathrm{Ker}(\mathbf{B}^\diamond) \oplus \mathrm{Ker}(\mathbf{A}^\diamond Q_b^\diamond) = \mathbb{R}_b^\diamond \qquad (17) $$

This property enables us to reformulate both interface conditions (2) in one single boundary equation:

$$ \mathbf{A}^{\diamond T} (\mathbf{A}^{\diamond} \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T})^{-1} \mathbf{A}^{\diamond} \boldsymbol{\mu}_b^{\diamond} - \mathbf{u}_b^{\diamond} = 0 \qquad (18) $$

indeed when multiplying this equation on the left by $\mathbf{B}^\diamond$ or $\mathbf{A}^\diamond Q_b^\diamond$ we recover classical interface conditions.

The mixed formulation of Problem (2) can be written as:

$$ \text{Find } \boldsymbol{\mu}_b^{\diamond} \in \mathbb{R}_b^{\diamond} \text{ such that } \mathbf{A}^{\diamond T} (\mathbf{A}^{\diamond} \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T})^{-1} \mathbf{A}^{\diamond} \boldsymbol{\mu}_b^{\diamond} - \mathbf{u}_b^{\diamond} = 0 \qquad (19) $$

$$ \text{where } \mathbf{u}^{\diamond} \text{ solves } \mathbf{f}_{\text{int}}^{\diamond}(\mathbf{u}^{\diamond}) - \mathbf{t}^{\diamond T} \mathbf{Q}_b^{\diamond} \mathbf{t}^{\diamond} \mathbf{u}^{\diamond} + \mathbf{t}^{\diamond T} \boldsymbol{\mu}_b^{\diamond} + \mathbf{f}_{\text{ext}}^{\diamond} = 0 $$

The last equation corresponds to the solution of independent nonlinear problems for each subdomain under Robin boundary conditions. If we assume that it has a unique solution for any given $\boldsymbol{\mu}_b^{\diamond}$ (the addition of matrix $\mathbf{Q}_b^\diamond$ suppresses rigid body motions), then we can define a nonlinear interface operator:

$$ \mathbf{u}_b^{(s)} = \mathbf{M}_{nl}^{(s)}(\boldsymbol{\mu}_b^{(s)}; \mathbf{f}_{\text{ext}}^{(s)}, \mathbf{Q}_b^{(s)}) \qquad (20) $$

In the linear case, the operator can be written as:

$$ \begin{gathered} \mathbf{M}_t^{(s)}(\boldsymbol{\mu}_b^{(s)}; \mathbf{f}_{\text{ext}}^{(s)}, \mathbf{Q}_b^{(s)}) = \mathbf{M}_t^{(s)} \boldsymbol{\mu}_b^{(s)} + \mathbf{b}_m^{(s)} \\ \text{with } \mathbf{M}_t^{(s)} = \mathbf{t}^{(s)}(\mathbf{K}^{(s)} + \mathbf{t}^{(s)T} \mathbf{Q}_b^{(s)} \mathbf{t}^{(s)})^{-1} \mathbf{t}^{(s)T} \\ \mathbf{b}_m^{(s)} = \mathbf{t}^{(s)}(\mathbf{K}^{(s)} + \mathbf{t}^{(s)T} \mathbf{Q}_b^{(s)} \mathbf{t}^{(s)})^{-1} \mathbf{f}_{\text{ext}}^{(s)} \end{gathered} \qquad (21) $$

One can recognize the operators of the FETI2LM approach [41].

The mixed nonlinear condensed problem can be written as:

$$ \text{Find } \boldsymbol{\mu}_b^{\diamond} \in \mathbb{R}_b^{\diamond} \text{ such that } \mathbf{A}^{\diamond T} (\mathbf{A}^{\diamond} \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T})^{-1} \mathbf{A}^{\diamond} \boldsymbol{\mu}_b^{\diamond} - \mathbf{M}_{nl}^{\diamond} (\boldsymbol{\mu}_b^{\diamond}; \mathbf{f}_{\text{ext}}^{\diamond}, \mathbf{Q}_b^{\diamond}) = 0 \qquad (22) $$

Note that the method makes use of the operator $(\mathbf{A}^{\diamond} \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T})$ which has exactly the structure of an assembled condensed stiffness matrix. The factorization of such a matrix is expensive in the general case, but we can choose a specific fill-in of matrix $\mathbf{Q}_b^\diamond$ for the assembled matrix $(\mathbf{A}^{\diamond} \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T})$ to have a block diagonal structure, which makes its handling much cheaper.

Note that an equivalent formulation is possible where the boundary unknowns $\boldsymbol{\mu}_b^{\diamond}$ are replaced by two interface unknowns: one balanced force $\gamma_B$ and one continuous displacement fields $v_A$:

$$ \boldsymbol{\mu}_b^{\diamond} = \mathbf{B}^{\diamond T} \gamma_B + \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T} v_A \qquad (23) $$

This is due to (17) which gives the following decomposition, similar to (5):

$$ \forall \mathbf{x}_b^{\diamond} \in \mathbb{R}_b^{\diamond}, \exists (\mathbf{x}_B, \mathbf{x}_A) \in \mathbb{R}^{\Gamma_B} \times \mathbb{R}^{\Gamma_A} / \mathbf{x}_b^{\diamond} = \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T} \mathbf{x}_A + \mathbf{B}^{\diamond T} \mathbf{x}_B \\ \text{indeed} \quad \left\{ \begin{array}{l} \mathbf{x}_A = (\mathbf{A}^{\diamond} \mathbf{Q}_b^{\diamond} \mathbf{A}^{\diamond T})^{-1} \mathbf{A}^{\diamond} \mathbf{x}_b^{\diamond} \\ \mathbf{x}_B = (\mathbf{B}^{\diamond} \mathbf{Q}_b^{\diamond}^{-1} \mathbf{B}^{\diamond T})^{+} \mathbf{B}^{\diamond} \mathbf{Q}_b^{\diamond}^{-1} \mathbf{x}_b^{\diamond} \end{array} \right. \qquad (24) $$

## 4.4 Note on the existence of nonlinear Schur complements

In previous sections we assumed the existence of nonlinear primal/dual/mixed Schur complement. For a well-posed global mechanical problem, local existence of such subdomain operator is most probable (see for instance [5] for a list of nonlinear frameworks where solutions exist). Yet the formulation, the material and the shape of the subdomain may strongly limit the domain of existence and uniqueness of the Dirichlet/Neumann/Robin problems: typically the large displacements hypothesis may prevent us from using too large compressive loads because of the possibility of buckling, as well as damage may prevent us from using too large tensile loads (see Section 7 for an illustration of these difficulties).
The global existence of the operators can be proved in the case of coercive continuous monotone operators; see for instance [42, 43] for an analysis at the level of the variational formulation and [14] for the analysis of the finite element approximation. Mechanically this framework is associated with positive hardening behaviors and certain contact laws, in small strains [29, 30]. In that context, Dirichlet, Neumann (assuming balanced load) and Robin problems are well-posed: the solution exists, is unique and varies continuously with respect to the load. Thanks to the continuity of the trace, the Schur operators can be defined. Moreover, the Schur complements inherit properties from the global problem (typically monotonicity, coercivity and continuity; see for instance [20] and associated bibliography), which makes the global condensed problems (9,15,22) well-posed.

# 5 Solution strategy

At this point, assuming the existence of local “Schur-type” nonlinear operators on subdomains, we have obtained one global nonlinear interface problem with an additive structure since it can be written as the assembly of subdomain contributions:

$$L(x) := \sum_{s} L^{(s)}(x) = 0 \qquad (25)$$

We propose to apply a Newton-Raphson procedure to that system. Iteration $k$ consists in solving:

$$\frac{\partial L}{\partial x}(x_k) \dot{x}_k = -L(x_k) \quad \text{and} \quad x_{k+1} = x_k + \dot{x}_k \\ \text{i.e.} \quad \left( \sum_s \frac{\partial L^{(s)}}{\partial x}(x_k) \right) \dot{x}_k = - \sum_s L^{(s)}(x_k) \quad \text{and} \quad x_{k+1} = x_k + \dot{x}_k \qquad (26)$$

We observe that the right-hand side, which corresponds to the evaluation of the residual, is obtained as the assembly of the result of independent solutions to nonlinear problems for each subdomain; whereas the left-hand side corresponds to the construction of the tangent operator as the assembly of subdomains' tangent operators.

The key point of this solution procedure is that the tangent condensed operators can actually be computed: indeed the nonlinear condensed operators are defined as an assembly of the trace (which are linear operations) of the solution to subdomains' problems, the linearization of which leads to the definition of their tangent stiffness matrices. In other words, the tangent of the nonlinear condensed operator (primal, dual or mixed) is the condensation of tangent subdomain operators.

We now explain in more detail the computation of the tangent operator for the primal approach. Let $S_t^{(s)}$ be the tangent operator of the nonlinear Schur complement $S_{nl}^{(s)}$. It can be defined as the unique linear operator satisfying the following relation for any $\dot{u}_b^{(s)}$:

$$\dot{\lambda}_b^{(s)} := S_{nl}^{(s)}(\dot{u}_{b_k}^{(s)} + \dot{u}_b^{(s)}; f_{ext}^{(s)}) - S_{nl}^{(s)}(u_{b_k}^{(s)}; f_{ext}^{(s)}) = S_{t_k}^{(s)}\dot{u}_b^{(s)} + o(\dot{u}_b^{(s)})$$

The associated subdomain problems can be written as:

$$f_{int}^{(s)}(u_k^{(s)}) + f_{ext}^{(s)} + t^{(s)T} \lambda_{b_k}^{(s)} = 0 \quad \text{with} \quad t^{(s)}u_k^{(s)} = u_{b_k}^{(s)}$$

The differentiation around $u_k^{(s)}$ involves the tangent stiffness matrix $K_{t_k}^{(s)}$:

$$-K_{t_k}^{(s)} \dot{u}_k^{(s)} + t^{(s)T} \lambda_b^{(s)} + o(\dot{u}_b^{(s)}) = 0 \quad \text{with} \quad t^{(s)} \dot{u}_k^{(s)} = \dot{u}_b^{(s)}$$

If we assume the well-posedness of the tangent Dirichlet problem $(\| (K_{t_k ii}^{(s)})^{-1} \| < \infty)$, then we can condense the previous equation without amplifying the negligible terms and then we obtain that

$$\dot{\lambda}_b^{(s)} = S_{t_k}^{(s)} \dot{u}_b^{(s)} + o(\dot{u}_b^{(s)}) \\ = (K_{t_k bb}^{(s)} - K_{t_k bi}^{(s)} (K_{t_k ii}^{(s)})^{-1} K_{t_k ib}^{(s)}) \dot{u}_b^{(s)} + o(\dot{u}_b^{(s)} - K_{t_k bi}^{(s)} (K_{t_k ii}^{(s)})^{-1} \dot{u}_i^{(s)})$$

Hence the identification between the tangent Schur complement and the Schur complement of the tangent stiffness matrix.

In the following, we derive the important steps of the different strategies.
## 5.1 Primal formulation

In the primal case, applying the Newton algorithm to equation (9) leads to:

$$ (A^{\diamond} S_{t_k}^{\diamond} A^{\diamond T}) \dot{u}_{A_k} = -A^{\diamond} S_{nl}^{\diamond} (A^{\diamond T} u_{A_k}; f_{ext}^{\diamond}) \quad (27) $$

where $S_{t_k}^{(s)} = (\frac{\partial S_{nl}^{(s)}}{\partial u_k^{(s)}}(u_{b_k}^{(s)}; f_{ext}^{(s)}))$ is the subdomain tangent primal Schur complement as given in (8).

First, let us focus on the right-hand side. It corresponds to the evaluation of a nonlinear Dirichlet problem for each subdomain with imposed interface displacements $u_{A_k}$:

$$ u_k^{\diamond} \text{ such that } \left\{ \begin{array}{l} (f_{int}^{\diamond}(u_k^{\diamond}) + f_{ext}^{\diamond})_i = 0 \\ t^{\diamond} u_k^{\diamond} = A^{\diamond T} u_{A_k} \end{array} \right. , \quad \lambda_{b_k}^{\diamond} := -(f_{int}^{\diamond}(u_k^{\diamond}) + f_{ext}^{\diamond})_b \quad (28) $$

From the parallel solution of these systems, we obtain the internal displacements $u_k^{\diamond}$ and the reactions $\lambda_{b_k}^{\diamond}$ whose lack of balance $A^{\diamond} \lambda_{b_k}^{\diamond} = A^{\diamond} S_{nl}^{\diamond} (A^{\diamond T} u_{A_k}; f_{ext}^{\diamond})$ is the right-hand side of the tangent system.

The left-hand side of (27) is the assembly of the subdomain's tangent primal Schur complements, that is to say a classical primal domain decomposition formulation of the tangent problem.

Thus, after an arbitrary initialization, the straightforward solution consists in repeating the following steps: (i) solving independent nonlinear Dirichlet systems for each subdomain, (ii) computing the lack of balance at the interface, and (iii) solving the global tangent interface problem by the application of the BDD algorithm [33, 32].

## 5.2 Dual formulation

In the dual case, applying Newton algorithm to equation (15) leads to:

$$ \left\{ \begin{array}{l} (B^{\diamond} F_{t_k}^{\diamond} B^{\diamond T}) \dot{\lambda}_{B_k} + B^{\diamond} R_b^{\diamond} \dot{\alpha}_k^{\diamond} = -B^{\diamond} (F_{nl}^{\diamond} (B^{\diamond T} \lambda_{B_k}; f_{ext}^{\diamond}) + R_b^{\diamond} \alpha_k^{\diamond}) \\ R_b^{\diamond T} B^{\diamond T} \dot{\lambda}_{B_k} = 0 \end{array} \right. \quad (29) $$

where $F_{t_k}^{(s)} = (\frac{\partial F_{nl}^{(s)}}{\partial \lambda_b^{(s)}}(\lambda_{b_k}^{(s)}; f_{ext}^{(s)}))$ is the subdomains' tangent dual Schur complement in configuration $k$ as given in (13).

As usually done in the FETI method, we introduce a projector $P_B$ on $Ker(R_b^{\diamond T} B^{\diamond T})$, so that we can seek $P_B \dot{\lambda}_{B_k}$ instead of $\dot{\lambda}_{B_k}$, solution to:

$$ P_B^T (B^{\diamond} F_{t_k}^{\diamond} B^{\diamond T}) P_B \dot{\lambda}_{B_k} = -P_B^T B^{\diamond} F_{nl}^{\diamond} (B^{\diamond T} \lambda_{B_k}; f_{ext}^{\diamond}) \quad (30) $$

and the contribution of the rigid body motions is sought after $\lambda_B$ is determined.

The right-hand side of (29) corresponds to the evaluation of nonlinear Neumann problems for each subdomain with imposed interface reaction $\lambda_{B_k}$:

$$ f_{int}^{\diamond}(u_k^{\diamond}) + f_{ext}^{\diamond} + t^{\diamond T} B^{\diamond T} \lambda_{B_k} = 0 \quad (31) $$

From the parallel solution of these systems, we obtain the displacements $u_k^{\diamond}$ whose interface gap $P_B^T B^{\diamond} t^{\diamond} u_k^{\diamond} = P_B^T B^{\diamond} F_{nl}^{\diamond} (B^{\diamond T} \lambda_{B_k}; f_{ext}^{\diamond})$ is the right-hand side of the tangent system.

The tangent matrix is the assembly of local tangent dual Schur complements computed at the value $u_k^{\diamond}$ obtained when computing the right-hand side. Then the first term of (29) is the assembly of the subdomain's tangent dual Schur complements, that is to say a classical dual domain decomposition formulation of the tangent problem.

In order for system (31) to be well posed, the rigid-body admissibility condition (11) has to be satisfied. This is realized by initializing $\lambda_B$ by an admissible reaction $\lambda_{B_0}$ which verifies:

$$ R^{\diamond T} (f_{ext}^{\diamond} + t^{\diamond T} B^{\diamond T} \lambda_{B_0}) = 0 \quad (32) $$

Indeed, since increments $\dot{\lambda}_{B_k}$ satisfy a zero-admissibility condition (29), and because the admissibility is a linear condition thanks to the small perturbation hypothesis, any coming $\lambda_{B_k}$ will satisfy the admissibility condition. Finding an initial value $\lambda_{B_0}$ corresponds to the initialization of the FETI method.

Thus, the straightforward solution consists in solving a FETI coarse grid problem to initialize $\lambda_{B_0}$, then repeating the following steps: (i) solving independent nonlinear Neumann systems for each subdomain, (ii) computing the displacement gap at the interface, and (iii) solving the global tangent interface problem by the application of FETI algorithm [12]. Note that the initialization can be improved by a full solution of the initial FETI system.
## 5.3 Mixed formulation

In the mixed case applying Newton algorithm to equation (22) leads to:

$$ (A^{\diamond T} (A^{\diamond} Q_b^{\diamond} A^{\diamond T})^{-1} A^{\diamond} - M_{t_k}^{\diamond}) \dot{\mu}_{b_k}^{\diamond} = M_{nl}^{\diamond} (\mu_{b_k}^{\diamond}; f_{ext}^{\diamond}, Q_b^{\diamond}) - A^{\diamond T} (A^{\diamond} Q_b^{\diamond} A^{\diamond T})^{-1} A^{\diamond} \mu_{b_k}^{\diamond} \quad (33) $$

where $M_{t_k}^{(s)} = (\frac{\partial M_{nl}^{(s)}}{\partial \mu_b^{(s)}}(\mu_{b_k}^{(s)}; f_{ext}^{(s)}, Q_b^{(s)}))$ is the subdomains' tangent mixed Schur complement in configuration k as given in (21).

The first term of the right-hand side of (33) corresponds to the evaluation of nonlinear Robin problems for each subdomain with interface impedance $Q_b^{\diamond}$ and imposed interface reaction $\mu_{b_k}^{\diamond}$:

$$ f_{int}^{\diamond}(\mathbf{u}_k^{\diamond}) - t^{\diamond T} Q_b^{\diamond} t^{\diamond} \mathbf{u}_k^{\diamond} + t^{\diamond T} \mu_{b_k}^{\diamond} + f_{ext}^{\diamond} = 0 \quad (34) $$

from the parallel solution of these systems, we obtain the displacement $\mathbf{u}_k^{\diamond}$ from which we deduce the mixed residual $b_{m_k}^{\diamond} = \mathbf{u}_k^{\diamond} - A^{\diamond T} (A^{\diamond} Q_b^{\diamond} A^{\diamond T})^{-1} A^{\diamond} \mu_{b_k}^{\diamond}$ which is the right-hand side of the tangent system.

The tangent operator is obtained by mixed condensation of the tangent stiffness matrix computed at the value $\mathbf{u}_k^{\diamond}$ obtained when evaluating the right-hand side. Then the first term of (33) is a classical mixed domain decomposition formulation of the tangent problem.

Thus after an arbitrary initialization, the straightforward solution consists in repeating the following steps: (i) solving independent nonlinear Robin systems for each subdomain, (ii) computing the mixed residual, and (iii) solving the global tangent interface problem by the application of FETI-2LM algorithm [41].

## 5.4 Swapping the linear solvers

In any case, the solution procedure alternates independent nonlinear solutions for each subdomain with different boundary conditions, and global tangent solutions similar to a domain decomposition method applied to linear problems. As we will show later, the choice of the boundary condition strongly impacts the nonlinear subdomains' solving. As the tangent operator is always associated with the tangent stiffness matrix of subdomains, it is possible to choose any formulation for the linear problem.

Indeed, starting from (8,13,21) we have the following classical relationships [17]:

$$ \begin{aligned} F_t^{\diamond} &= S_t^{\diamond +} \\ M_t^{\diamond} &= (S_t^{\diamond} + Q_b^{\diamond})^{-1} \\ b_d^{\diamond} &= F_t^{\diamond} b_p^{\diamond} \\ b_m^{\diamond} &= M_t^{\diamond} b_p^{\diamond} \end{aligned} \quad (35) $$

This implies that even if a formulation results in a specific residual after the subdomains' nonlinear solutions, any linear solver can be employed. Typically, the mixed formulation has interesting capabilities for nonlinear problems, but it is not very classical for linear system for which primal (BDD) or dual (FETI) formulations are more standard and embed powerful preconditioners. After computing $M_{nl}^{\diamond}(\mu_{b_k}^{\diamond}; f_{ext}^{\diamond}, Q_b^{\diamond})$, one can compute $b_{m_k}^{\diamond} = M_{nl}^{\diamond}(\mu_{b_k}^{\diamond}; f_{ext}^{\diamond}, Q_b^{\diamond}) - A^{\diamond T} (A^{\diamond} Q_b^{\diamond} A^{\diamond T})^{-1} A^{\diamond} \mu_{b_k}^{\diamond}$, deduce the equivalent primal right-hand side $b_{p_k}^{\diamond} = (S_{t_k}^{\diamond} + Q_b^{\diamond}) b_{m_k}^{\diamond}$ and use the associated interface lack of balance $A^{\diamond} b_{p_k}^{\diamond}$ as input to the tangent BDD solver. This property can directly be observed by introducing (23) in the mixed tangent system (33), since after premultiplying by $(Q_b^{\diamond} + S_{t_k}^{\diamond})$ we obtain:

$$ \begin{aligned} S_{t_k}^{\diamond} A^{\diamond T} \dot{v}_{A_k} - B^{\diamond T} \dot{\gamma}_{B_k} &= (Q_b^{\diamond} + S_{t_k}^{\diamond}) (M_{nl}^{\diamond}(\mu_{b_k}^{\diamond}; f_{ext}^{\diamond}, Q_b^{\diamond}) - A^{\diamond T} v_{A_k}) \\ &= (Q_b^{\diamond} + S_{t_k}^{\diamond}) b_{m_k}^{\diamond} = b_{p_k}^{\diamond} \end{aligned} \quad (36) $$

which leads to the system set in terms of $\dot{v}_{A_k}$:

$$ (A^{\diamond} S_{t_k}^{\diamond} A^{\diamond T}) \dot{v}_{A_k} = A^{\diamond} b_{p_k}^{\diamond} \quad (37) $$

The following quantities can be deduced (even from an inexact solution):

$$ \begin{aligned} \dot{\mu}_{b_k}^{\diamond} &= S_{t_k}^{\diamond} A^{\diamond T} \dot{v}_{A_k} - b_{p_k}^{\diamond} \\ \dot{u}_k^{\diamond} &= (K_{t_k}^{\diamond} + t^{\diamond T} Q_b^{\diamond} t^{\diamond})^{-1} t^{\diamond T} (b_{p_k}^{\diamond} + \dot{\mu}_{b_k}^{\diamond}) \\ \dot{u}_{b_k}^{\diamond} &= t^{\diamond T} \dot{u}_k^{\diamond} \\ \dot{\lambda}_{b_k}^{\diamond} &= S_{t_k}^{\diamond} \dot{u}_{b_k}^{\diamond} - b_{p_k}^{\diamond} = \dot{\mu}_{b_k}^{\diamond} - Q_b^{\diamond} \dot{u}_{b_k}^{\diamond} \end{aligned} \quad (38) $$
Note that $\dot{\gamma}_{B_k}$ plays no role in the algorithm but it could also be post-processed:

$$ \dot{\gamma}_{B_k} = B^{\diamond T +} \left( \dot{\mu}_{b_k}^{\diamond} - Q_b^{\diamond} A^{\diamond T} \dot{v}_{A_k}^{\diamond} \right) \qquad (39) $$

where $B^{\diamond T +}$ is a scaled assembling operator commonly employed in the preconditioning step of the FETI method [40, 28]; because of the space splitting properties (24), $\dot{\gamma}_{B_k}$ does not depend on the choice of the scaling. Of course, if a dual formulation had been preferred then $\dot{\gamma}_{B_k}$ would have been computed by a FETI-like system (left-multiply (36) by $B^{\diamond F_t^{\diamond}}$).

## 5.5 Typical algorithm

Algorithm 1 sums up the main steps of the method with the mixed nonlinear local problems and primal global solver. For simplicity reasons, only one increment was considered.

As can be seen in this algorithm, besides classical global convergence criterion $\varepsilon_{NG}$ (as in a standard Newton approach), two precision thresholds are used: the local nonlinear thresholds $\varepsilon_{NL}^{\diamond}$ (associated with the Newton processes carried out independently on subdomains) and the global linear threshold of the domain decomposition (Krylov) solver $\varepsilon_K$ (here BDD). The convergence criteria are discussed in more detail in the following section.

Other parameters can be the initializations of the various iterative solvers and, for the mixed approach, the choice of the impedance matrices $Q_b^{\diamond}$.

<caption>**Algorithm 1:** Mixed nonlinear approach with BDD tangent solver</caption>

**Define:**

$$ r_{nl}^{m\diamond}(u^{\diamond}, \mu_b^{\diamond}) = f_{int}^{\diamond}(u^{\diamond}) - t^{\diamond T} Q_b^{\diamond} t^{\diamond} u^{\diamond} + t^{\diamond T} \mu_b^{\diamond} + f_{ext}^{\diamond} $$

**Initialization:**

$$ (u_0^{\diamond}, \lambda_{b_0}^{\diamond}) \text{ such that } B^{\diamond t^{\diamond}} u_0^{\diamond} = 0 \text{ and } A^{\diamond} \lambda_{b_0}^{\diamond} = 0 $$

$$ \text{Set } k = 0 $$

$$ \text{Define } \mu_{b_k}^{\diamond} = \lambda_{b_k}^{\diamond} + Q_b^{\diamond} t^{\diamond} u_k^{\diamond} $$

$$ \text{while } \|r_{nl}^{m\diamond}(u_k^{\diamond}, \mu_{b_k}^{\diamond})\| + \|B^{\diamond t^{\diamond}} u^{\diamond}\|_B > \varepsilon_{NG} \text{ do} $$

**Local nonlinear step:**

$$ \text{Set } u_{k,0}^{\diamond} = u_k^{\diamond} \text{ and } j = 0 $$

$$ \text{while } \|r_{nl}^{m\diamond}(u_{k,j}^{\diamond}, \mu_{b_k}^{\diamond})\| > \varepsilon_{NL}^{\diamond} \text{ do} $$

$$ \left| \begin{array}{l} u_{k,j+1}^{\diamond} = u_{k,j}^{\diamond} - (K_{t_{k,j}}^{\diamond} + t^{\diamond T} Q_b^{\diamond} t^{\diamond})^{-1} r_{nl}^{m\diamond}(u_{k,j}^{\diamond}, \mu_{b_k}^{\diamond}) \\ \text{Set } j = j + 1 \end{array} \right. $$

**end**

**Linear right-hand side:**

$$ b_{m_k}^{\diamond} = A^{\diamond T} (A^{\diamond} Q_b^{\diamond} A^{\diamond T})^{-1} A^{\diamond} \mu_{b_k}^{\diamond} - t^{\diamond} u_{k,j}^{\diamond} $$

$$ b_{p_k}^{\diamond} = (S_{t_{k,j}}^{\diamond} + Q_b^{\diamond}) b_{m_k}^{\diamond} $$

**Global linear step:**

$$ \text{Set } \dot{v}_{A_k}^{0} = 0 \text{ and } i = 0 $$

$$ \text{while } \|b_{p_k}^{\diamond} - (A^{\diamond} S_{t_{k,j}}^{\diamond} A^{\diamond T})^{\circ} \dot{v}_{A_k}^{i}\| > \varepsilon_K \text{ do} $$

| Make BDD iterations (index $i$)

**end**

$$ \text{Set } u_{k+1}^{\diamond} = u_k^{\diamond} + \dot{u}_k^{i\diamond} \text{ and } \lambda_{b_{k+1}}^{\diamond} = \lambda_{b_k}^{\diamond} + \dot{\lambda}_{b_k}^{i\diamond} \text{ using (38)} $$

$$ \text{Set } k = k + 1 $$

**end**

## 6 Error analysis

The aim of this section is to further analyze the effect of the different thresholds and the way they are connected. When possible, we will make reference to the global residual $r$ which is defined as:

$$ f_{int}(u) + f_{ext} = r \qquad (40) $$

The objective of the solver is that $\|r\| < \varepsilon_{NG} \|f_{ext}\|$ (global nonlinear criterion).

We first show how the convergence criteria for the Local Newton solvers $\varepsilon_{NL}^{\diamond}$ and for the Krylov interface solver $\varepsilon_K$ are related for the different approaches. Then we interpret the whole process as an inexact Newton solver [7] and
we consider adapting at each step the criteria $\varepsilon_K$ and $\varepsilon_{NL}^{\Diamond}$ to the current error status $\|\mathbf{r}\|$. $\varepsilon_K$ and $\varepsilon_{NL}^{\Diamond}$ are interpreted as forcing terms which need to be close to the objective $\varepsilon_{NG}$ only when global convergence is almost reached whereas they should be relaxed at the beginning of the process in order to avoid oversolving (see Section 6.4).

# 6.1 Primal approach

## 6.1.1 Local Newton solvers

The inner (local) Newton loops are associated with Dirichlet problems. Assuming these conditions are exactly taken into account, the convergence is controlled by the internal node residue:

$$f_{int_i}^{\Diamond}(\mathbf{u}^{\Diamond}) + f_{ext_i}^{\Diamond} = \mathbf{r}_i^{\Diamond} \text{ with } \|\mathbf{r}_i^{\Diamond}\| < \varepsilon_{NL}^{\Diamond} \|f_{ext_i}^{\Diamond}\| \quad (41)$$

From the $\mathbf{u}^{\Diamond}$ computed above, the reactions are defined by the following relation:

$$\lambda_b^{\Diamond} := -(f_{int_b}^{\Diamond}(\mathbf{u}^{\Diamond}) + f_{ext_b}^{\Diamond}) \quad (42)$$

## 6.1.2 Link to the global error

In the primal case, the displacement search space is the same as in the non-substructured case (displacements are continuous across the interface). This simplifies the error analysis. Indeed, we directly have:

$$\|f_{int}(\mathbf{u}) + f_{ext}\| \le \|f_{int_i}^{\Diamond}(\mathbf{u}^{\Diamond}) + f_{ext_i}^{\Diamond}\| + \|A^{\Diamond}(f_{int_b}^{\Diamond}(\mathbf{u}^{\Diamond}) + f_{ext_b}^{\Diamond})\| \\ \le \|\varepsilon_{NL}^{\Diamond}\|\|f_{ext_i}^{\Diamond}\| + \|A^{\Diamond}\lambda_b^{\Diamond}\| \quad (43)$$

This means that global convergence occurs when both the local Newton criteria and the resulting lack of balance at the interface (the initial residual of the tangent interface solver) are small.

## 6.1.3 Tangent interface solver

From the point of view of the nonlinear Schur complement, the computation of the reactions suffers from the error made on internal displacements, if $S_{nl}^{\Diamond}(A^{\Diamond T}\mathbf{u}_{A_k}; f_{ext}^{\Diamond})$ denotes the exact reaction (corresponding to $\mathbf{r}_i^{\Diamond} = 0$), we write:

$$\lambda_b^{\Diamond} = S_{nl}^{\Diamond}(A^{\Diamond T}\mathbf{u}_{A_k}; f_{ext}^{\Diamond}) + \mathbf{r}_b^{\Diamond} \quad (44)$$

In the linear case, the internal residual would propagate on the boundary according to the relation $\mathbf{r}_b^{\Diamond} = -K_{bi}^{\Diamond}K_{ii}^{\Diamond-1}\mathbf{r}_i^{\Diamond}$, so that we can assume that the error on the reaction is of the same order of the inner solver error $\|\mathbf{r}_b^{\Diamond}\| \approx \|\mathbf{r}_i^{\Diamond}\|$.

We then can study the error made on the tangent solution. The computed vector $\mathbf{u}_{A_k}^{\circ}$ satisfies the equation:

$$(A^{\Diamond} S_{lk}^{\Diamond} A^{\Diamond T}) \mathbf{u}_{A_k}^{\circ} = -A^{\Diamond} S_{nl}^{\Diamond} (A^{\Diamond T} \mathbf{u}_{A_k}; f_{ext}^{\Diamond}) - A^{\Diamond} \mathbf{r}_b^{\Diamond} + \mathbf{r}_A^{\Diamond} \quad (45)$$

where $\|\mathbf{r}_A^{\Diamond}\| < \varepsilon_K \|\mathbf{r}_{A_0}^{\Diamond}\|$ controls the error of the linear (Krylov) solver. Since $\|A^{\Diamond}\| = O(1)$, we see that the error of the local Newton solver and the global interface iterative solver should be of the same order of magnitude.

# 6.2 Dual approach

## 6.2.1 Local Newton solvers

The inner (local) Newton loops are associated with Neumann problems, and the convergence is controlled by the following residue:

$$f_{int}^{\Diamond}(\mathbf{u}^{\Diamond}) + f_{ext}^{\Diamond} + t^{\Diamond T} B^{\Diamond T} \lambda_B = \mathbf{r}^{\Diamond} \text{ with } \|\mathbf{r}^{\Diamond}\| < \varepsilon_{NL}^{\Diamond} \|f_{ext}^{\Diamond}\| \quad (46)$$

The output from these computations is the interface displacement. Unfortunately there is no direct control on the associated error $\delta_b^{\Diamond}$. A good estimation is provided by the last correction brought to $\mathbf{u}_b^{\Diamond}$ by the inner Newton loop.

## 6.2.2 Link to the global error

The difficulty with the dual approach is that the (broken) search space is larger than the original space (where displacements are continuous). In the linear case, a costless processing enables us to obtain a continuous displacement [34]. Such a strategy is not realistic for nonlinear problems.

One possibility is to define an error extended to the broken space ($e_{bs}$):

$$e_{bs}^2 = \|f_{int}^{\Diamond}(\mathbf{u}^{\Diamond}) + f_{ext}^{\Diamond} + t^{\Diamond T} B^{\Diamond T} \lambda_B\|^2 + \|B^{\Diamond} t^{\Diamond} \mathbf{u}^{\Diamond}\|^2_B \quad (47)$$

where $\|\cdot\|_B$ is a well chosen norm for the interface discontinuity. We recover an expression where the inner criterion is cumulated with the initial residual of the tangent interface solver.
### 6.2.3 Tangent interface solver

The tangent interface solution can be written as:

$$P_B^T (B^{\diamond} F_{t_k}^{\diamond} B^{\diamond T}) P_B \lambda_{B_k} = -P_B^T (B^{\diamond} F_{nl}^{\diamond} (B^{\diamond T} \lambda_{B_k}; f_{ext}) + B^{\diamond} \delta_b^{\diamond}) + r_B^L \quad (48)$$

where $\|r_B^L\| < \varepsilon_K \|r_{B_0}^L\|$ controls the error of the linear solver. Since $\|B^{\diamond}\| = O(1)$, we see that the forward error of the local Newton solver and the backward error of the global interface iterative solver should be of the same order of magnitude.

## 6.3 Mixed approach

### 6.3.1 Local Newton solvers

The inner (local) Newton loops are associated with Robin problems. The convergence is controlled by the following residue:

$$f_{int}^{\diamond}(u^{\diamond}) - t^{\diamond T} Q_b^{\diamond} t^{\diamond} u^{\diamond} + t^{\diamond T} \mu_b^{\diamond} + f_{ext}^{\diamond} = r^{\diamond} \quad \text{with} \quad \|r^{\diamond}\| < \varepsilon_{NL} \|f_{ext}^{\diamond}\| \quad (49)$$

As in the dual approach, the output from these computations is the interface displacement, of which the error $\delta_b^{\diamond}$ can be estimated by the last correction brought to $u_b^{\diamond}$ by the inner Newton loop.

### 6.3.2 Link to the global error

As in the dual case, the search space is a broken space. This implies the use of an extended norm to evaluate the global error: not only the displacements are not continuous at the interface but neither are the reactions balanced. With reaction $\lambda_b^{\diamond}$ defined by (42), we see that the trace of the local residue can be written as:

$$r_b^{\diamond} = Q_b^{\diamond} (A^{\diamond T} v_A - t^{\diamond} u^{\diamond}) + (B^{\diamond T} \gamma_B - \lambda_b^{\diamond}) \quad (50)$$

This relation links the two interface errors, which means that the dual criterion (47) is sufficient to monitor the convergence: if both $\|r_b^{\diamond}\|$ and $\|B^{\diamond}t^{\diamond}u^{\diamond}\|_B$ are small, then so must $\|A^{\diamond}\lambda_b^{\diamond}\|$ be.

### 6.3.3 Tangent interface solver

The tangent interface solution can be written as (in primal form):

$$(A^{\diamond} S_{t_k}^{\diamond} A^{\diamond T}) \dot{v}_{A_k} = A^{\diamond} (Q_b^{\diamond} + S_{t_k}^{\diamond}) (u_{b_k}^{\diamond} + \delta_b^{\diamond} - A^{\diamond T} v_{A_k}) + r_A^L \quad (51)$$

where $\|r_A^L\| < \varepsilon_K \|r_{A_0}^L\|$ controls the error of the linear solver. From the definition of $\delta_b^{\diamond}$, we can expect $\|A^{\diamond}(Q_b^{\diamond} + S_{t_k}^{\diamond})\delta_b^{\diamond}\| \simeq \|A^{\diamond}t^{\diamond}r^{\diamond}\|$ which means that the convergence criterion of the Krylov solver must be of the order of magnitude of the local Newton solvers' criteria.

## 6.4 Tuning of criteria

The previous subsections established connections between the convergence criteria of the local Newton solvers and of the global interface Krylov solver. We can now analyze the algorithm in the framework of inexact Newton methods [7]. From previous analysis, and using the notations of (26), we solve (25) by the sequence:

$$\frac{\partial L}{\partial x}(x_k) \dot{x}_k = -L(x_k) + r_k, \quad \text{with} \quad \|r_k\| \le \varepsilon_k \|L(x_k)\| \quad (52) \\ x_{k+1} = x_k + \dot{x}_k$$

$\varepsilon_k$ is called the forcing term, and its value can be updated at each Newton iteration in order to optimize the convergence rate. More precisely, [9] shows that:

1. If $0 \le \varepsilon_k \le \varepsilon_{max} < 1$ $\forall k$, then the series converges linearly, with an asymptotic rate $a \le \varepsilon_{max}$.

2. If $\lim_{k \to +\infty} \varepsilon_k = 0$, convergence is superlinear.

3. If $\varepsilon_k = O(\|L(x_k)\|)$, convergence is quadratic.
These convergence properties lead to more or less simple expressions for the forcing terms (see [9]):

$$ \varepsilon_k = \frac{1}{2^{k+1}} \qquad (53) $$

$$ \varepsilon_k = \gamma \left( \frac{\| \mathbf{L}(\mathbf{x}_k) \|}{\| \mathbf{L}(\mathbf{x}_{k-1}) \|} \right)^\alpha, \quad \gamma \in [0, 1[, \alpha \in ]0, 2] \qquad (54) $$

$$ \varepsilon_k = \frac{\| \mathbf{L}(\mathbf{x}_k) \| - \| \mathbf{L}(\mathbf{x}_{k-1}) + \mathbf{L}'(\mathbf{x}_{k-1}) \mathring{\mathbf{x}}_{k-1} \|}{\| \mathbf{L}(\mathbf{x}_{k-1}) \|} \qquad (55) $$

Expression (53) is easy to implement but does not relate the linear criterion with $\|\mathbf{L}(\mathbf{x}_k)\|$ at each iteration: the convergence of the series $(\mathbf{x}_k)_k$ will theoretically not be better than superlinear. Expressions (54) and (55) are more complicated to implement, but can a priori allow faster convergence rates. For example taking $\alpha = 2$ in expression (54) should lead to quadratic convergence.

# 7 Assessments

The aim of these assessments is to test the various configurations of the approach to an academic problem with a realistic nonlinear behavior. We are aware that the methods we propose are sensitive to many parameters, among others the type of nonlinearities (with potential instabilities), the partitioning of the domain, the localization of the linearities (whether the nonlinearity spreads amongst subdomains or not), the variants of the solvers, the loading increments, the convergence criteria.

The following study is thus far from pretending to be exhaustive, it simply aims at giving major trends of the methods and at proving that they can be computationally interesting. At the end of the section, we discuss the applicability of the methods to other cases and the expected performance.

## 7.1 Perfect plasticity problem

We consider a rectangular plate in traction with a central hole (length $L_x = 95$ mm, height $L_y = 31.5$ mm, radius $r = 5$ mm) as depicted in Figure 2 (after taking into account the symmetries). The material is a perfectly plastic steel (Young modulus $E = 210$ GPa, Poisson's ratio $\nu = 0.3$, elasticity limit $\sigma_0 = 420$ MPa), and we assume the plane stress conditions. The mesh is constituted by triangle elements with three nodes, with characteristic length $h = 0.01 L_y$. The eight subdomains are delimited by concentric circles. They are numbered starting from the center.

<figure><img src="image_3.png" /><figcaption>Figure 2: Localized plasticity problem / numbering of subdomains</figcaption></figure>

Small perforations were added near the central hole in order to trigger a large level of localized plasticity. Traction loading is incrementally applied on the perforated plate by imposing a displacement $u_D$ on its right side. Increments are taken as follows, with $u_e = 25$ µm:

$$ u_D = [1 \ 2 \ 3 \ 4 \ 5 \ 5.75 \ 6.5] u_e $$

The first step, corresponding to an increment of loading $u_D = u_e$, remains elastic. Steps $u_D = 2u_e$ to $u_D = 5u_e$ activate plasticity in only the central subdomain, reaching a maximum of $p = 0.051$ of cumulated plasticity. Plasticity then spreads in the second and third subdomains until $\max(p) = 0.15$.

This sequence of loadings was identified using a classical (outer Newton) approach. It could be reused as is with primal and mixed approaches. Solving a local Neumann problem with perfectly plastic subdomains is known to be a tougher problem than a Dirichlet or a Robin problem. Indeed, it was necessary to adopt a finer loading sequence for the dual approach; the following steps needed to be inserted:

$$ u_D = [4.5 \ 5.25 \ 5.5 \ 6 \ 6.125 \ 6.25 \ 6.375] u_e $$

It was verified that, for the primal and mixed approaches, this refined sequence leads to results of comparable quality as the coarse sequence.
In the following, we compare the classical approach (outer Newton, inner Schur-Krylov), the primal and dual nonlinear approaches, and the mixed nonlinear approach with parameter $Q_b$ equal to the stiffness of the neighbors, or $Q_b$ equal to the optimal choice (Schur complement of the remainder of the structure). Of course the latter choice is not computationally realistic, and efficient approximations will be the subject of further studies.

In order to assess the performance, we measure:

* The number of outer Newton loops, since nonlinear approaches shall enable us to follow a faster path to convergence,

* The maximum number of computations of the tangent amongst subdomains. For the classical approach, it is simply equal to the number of outer iterations. The inner Newton loops makes that number larger in non-linear approaches.

* The number of cumulated Krylov (conjugate gradient) iterations which is proportional to the number of exchanges between subdomains.

The aim of the nonlinear approaches is to lower the number of outer Newton and Krylov iterations, at the price of some extra local Newton iterations. Of course the number of computations of the tangent varies from one subdomain to another, depending on the spread and intensity of the nonlinearity, which can cause poor load balancing. This is a real issue which will be the subject of future studies. Nevertheless we believe that the reduction of the number of outer iterations and of the exchanges (Krylov iterations) already makes the proposed methods interesting.

### 7.1.1 Classical technique

We briefly give the performance for the classical Newton method with a domain decomposition solver for the tangent systems. Note that BDD is used but FETI gives a quite similar number of Krylov iterations. The outer Newton's convergence criterion is set to $\epsilon_{NG} = 10^{-5}$. For the results given in Table 1, the inner Krylov solver criteria for the condensed system is fixed to $\epsilon_K = 10^{-6}$.

<table><caption>Table 1: Performance of the classical approach with BDD solver for the tangent systems, with fixed criteria.</caption><thead><tr><th>Increment of loading</th><th>u<sub>e</sub></th><th>2u<sub>e</sub></th><th>5u<sub>e</sub></th><th>5.75u<sub>e</sub></th><th>6.5u<sub>e</sub></th></tr></thead><tbody><tr><th>Spread of nonlinearity</th><td>Elastic</td><td>1SD plastifies</td><td></td><td>Several SD plastify</td><td></td></tr><tr><th>Global iterations</th><td>1</td><td>2</td><td>13</td><td>17</td><td>23</td></tr><tr><th>Cumulated Krylov iterations</th><td>19</td><td>38</td><td>255</td><td>336</td><td>460</td></tr></tbody></table>

### 7.1.2 Nonlinear localization with fixed criteria

Before tuning the different thresholds, as developed in section 6.4, a comparison is made between the classical outer Newton and the nonlinear approaches with fixed stopping criteria. In order to achieve a meaningful comparison, thresholds are fixed to typical values:

$$ \epsilon_{NG} = 10^{-5} \\ \epsilon_{NL} = \epsilon_K = 10^{-6} $$

Relative results of primal, dual and mixed approaches normalized by the classical approach are given in Table 2 for chosen loading increments.

Table 2 clearly indicates a decrease in the number of cumulated global iterations when using the primal and mixed approaches: 20% less iterations for the primal and $K_{bb}$-mixed approaches, 30% less iterations with the optimal Robin condition. For the dual method, the number of global iterations increases because of the steps inserted additionally to insure convergence.

The gain in Krylov iterations is quite similar to the gain in global Newton iterations, because each solution of a linear system requires an almost constant number of Krylov iterations (note that before requiring extra increments, the dual approach needs one Krylov iteration less than the primal approach, hence the 0.95 ≈ 18/19 value in the table). These figures emphasize the need to adapt the precision of the Krylov solver as a function of the current nonlinear residual.

Regarding local Newton iterations, we see that, except for the dual method which is inefficient in that case, no more than 3 times more computations of the tangent are required for the nonlinear approaches. In particular, the optimal Robin only requires 2.4 times more computations of the tangent.
<table><caption>Table 2: Ratios of cumulated iterations between nonlinearly localized methods and classic method, with fixed stopping criteria</caption><thead><tr><th></th><th>Increment of loading</th><th>u<sub>e</sub></th><th>2u<sub>e</sub></th><th>5u<sub>e</sub></th><th>5.75u<sub>e</sub></th><th>6.5u<sub>e</sub></th></tr><tr><th></th><th>Spread of nonlinearity</th><th>Elastic</th><th colspan="2">SD 1</th><th>SD 1-2</th><th>SD 1-3</th></tr></thead><tbody><tr><th rowspan="4">Global Newton</th><th>Primal/Classic</th><td>1</td><td>1</td><td>0.77</td><td>0.82</td><td>0.83</td></tr><tr><th>Dual/Classic</th><td>1</td><td>1</td><td>1.46</td><td>2.12</td><td>3.35</td></tr><tr><th>Mixed/Classic, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>1</td><td>0.77</td><td>0.82</td><td>0.78</td></tr><tr><th>Mixed/Classic, Q<sub>b</sub> opti</th><td>1</td><td>1</td><td>0.62</td><td>0.65</td><td>0.70</td></tr><tr><th rowspan="4">Krylov</th><th>Primal/Classic</th><td>1</td><td>1</td><td>0.76</td><td>0.82</td><td>0.82</td></tr><tr><th>Dual/Classic</th><td>0.95</td><td>0.95</td><td>1.36</td><td>1.96</td><td>3.08</td></tr><tr><th>Mixed/Classic, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>1</td><td>0.76</td><td>0.82</td><td>0.78</td></tr><tr><th>Mixed/Classic, Q<sub>b</sub> opti</th><td>1</td><td>1</td><td>0.61</td><td>0.64</td><td>0.69</td></tr><tr><th rowspan="4">Local Newton</th><th>Primal/Classic</th><td>1</td><td>2</td><td>2.62</td><td>2.88</td><td>3.04</td></tr><tr><th>Dual/Classic</th><td>1</td><td>2</td><td>5.69</td><td>9.24</td><td>16.74</td></tr><tr><th>Mixed/Classic, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>2</td><td>2.69</td><td>2.94</td><td>3</td></tr><tr><th>Mixed/Classic, Q<sub>b</sub> opti</th><td>1</td><td>2</td><td>1.85</td><td>2.12</td><td>2.43</td></tr></tbody></table>

As a conclusion for this test case, the primal and mixed versions of nonlinearly localized methods lead to interesting performance. Moreover they seem to be more effective when plasticity remains localized in only one subdomain. Indeed, global and Krylov numbers of iterations present a strong decrease when comparing them to those of the classic method (ratios from 0.61 to 0.76 at step $u_D = 5u_e$), which means a decrease in communications between subdomains. The cost of this improvement stays limited, with less than 3 times more factorizations at step $5u_e$ (and even less than 2 times more for optimal $Q_b$). This is interesting because, as previously said, local iterations are independent and do not require communications between processors, which results in a gain of time for large-scale solutions.

The dual approach does not provide significant results, since additional increments had to be used to achieve convergence. Thus, it will not be considered in the following part, where better performance is achieved with an adaptation of stopping criteria.

### 7.1.3 Adaptation of stopping criteria

In order to observe the influence of the tuning of $\varepsilon_K$ and $\varepsilon_{NL}$ with respect to the current convergence state, the two expressions of section 6.4 are tested with a set of coefficients chosen for their performance (see Table 3). As discussed in section 6, we choose $\varepsilon_K = \varepsilon_{NL}$. The same sequence of load increments as that of the previous section is applied to the perforated plate.

<table><caption>Table 3: Chosen parameters for $\varepsilon_K = \varepsilon_{NL}$</caption><tr><td>Choice 1</td><td>Eq (54)</td><td>γ = 0.7, α = 1.5, ε<sub>K0</sub> = 10<sup>-6</sup></td></tr><tr><td>Choice 2</td><td>Eq (55)</td><td>ε<sub>K0</sub> = 10<sup>-4</sup></td></tr></table>

The relative performance for choice 1 is summed up in Table 4, whereas the relative performance for choice 2 is given in Table 5. Depending on the quantity, we compare the variants with adapted criteria (INexact in the tables) to the fixed criteria approaches, or the nonlinear approaches with the classical approach, both with adapted criteria. Note that choice 2 was too aggressive for classic and primal nonlinear localized methods to converge, so that only the mixed approaches were studied.

The aim of adaption is to avoid useless inner iterations without modifying the global iterations. This absence of perturbation corresponds to ratios of 1 in the “Global Newton” rows. This happens almost every time.

Regarding the number of Krylov iterations, the figures in Table 4 are not as good, but still of the same order as the ones in Table 2. Table 5 shows that mixed approaches can be really efficient with a good adaption of criteria. This means that adaption is not in contradiction with nonlinear localization and shall be used to reduce the number of exchanges.
<table><caption>Table 4: Performance of adapted criteria with choice 1.</caption><thead><tr><th></th><th>Increment of loading</th><th>u<sub>e</sub></th><th>2u<sub>e</sub></th><th>5u<sub>e</sub></th><th>5.75u<sub>e</sub></th><th>6.5u<sub>e</sub></th></tr><tr><th></th><th>Spread of nonlinearity</th><th>Elastic</th><th colspan="2">SD 1</th><th>SD 1-2</th><th>SD 1-3</th></tr></thead><tbody><tr><th rowspan="4">Global Newton</th><th>Classic IN/Classic</th><td>1</td><td>1</td><td>1</td><td>1.06</td><td>1</td></tr><tr><th>Primal IN/Primal</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>Mixed IN/Mixed, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>1</td><td>1</td><td>1</td><td>1.06</td></tr><tr><th>Mixed IN/Mixed, Q<sub>b</sub> opti</th><td>1</td><td>1</td><td>1</td><td>1.09</td><td>1.06</td></tr><tr><th rowspan="4">Krylov</th><th>Classic IN/Classic</th><td>1</td><td>1</td><td>0.77</td><td>0.77</td><td>0.69</td></tr><tr><th>Primal IN/Classic IN</th><td>1</td><td>1</td><td>0.86</td><td>0.85</td><td>0.89</td></tr><tr><th>Mixed IN/Classic IN, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>1</td><td>0.86</td><td>0.85</td><td>0.91</td></tr><tr><th>Mixed IN/Classic IN, Q<sub>b</sub> opti</th><td>1</td><td>1</td><td>0.73</td><td>0.73</td><td>0.79</td></tr><tr><th rowspan="4">Local Newton</th><th>Classic IN/Classic</th><td>1</td><td>1</td><td>1</td><td>1.06</td><td>1</td></tr><tr><th>Primal IN/Classic IN</th><td>1</td><td>2</td><td>2.54</td><td>2.56</td><td>2.78</td></tr><tr><th>Mixed IN/Classic IN, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>2</td><td>2.31</td><td>2.33</td><td>2.57</td></tr><tr><th>Mixed IN/Classic IN, Q<sub>b</sub> opti</th><td>1</td><td>2</td><td>1.85</td><td>2.06</td><td>2.30</td></tr></tbody></table>

<table><caption>Table 5: Performance of adapted criteria with choice 2.</caption><thead><tr><th></th><th>Increment of loading</th><th>u<sub>e</sub></th><th>2u<sub>e</sub></th><th>5u<sub>e</sub></th><th>5.75u<sub>e</sub></th><th>6.5u<sub>e</sub></th></tr><tr><th></th><th>Spread of nonlinearity</th><th>Elastic</th><th colspan="2">SD 1</th><th>SD 1-2</th><th>SD 1-3</th></tr></thead><tbody><tr><th rowspan="2">Global Newton</th><th>Mixed IN/Mixed, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>1</td><td>1</td><td>1</td><td>1.06</td></tr><tr><th>Mixed IN/Mixed, Q<sub>b</sub> opti</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th rowspan="2">Krylov</th><th>Mixed IN/Classic IN, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>0.79</td><td>0.71</td><td>0.67</td><td>0.66</td></tr><tr><th>Mixed IN/Classic IN, Q<sub>b</sub> opti</th><td>1</td><td>0.79</td><td>0.58</td><td>0.56</td><td>0.57</td></tr><tr><th rowspan="2">Local Newton</th><th>Mixed IN/Classic IN, Q<sub>b</sub> = K<sub>bb</sub></th><td>1</td><td>2</td><td>2.23</td><td>2.28</td><td>2.43</td></tr><tr><th>Mixed IN/Classic IN, Q<sub>b</sub> opti</th><td>1</td><td>2</td><td>1.85</td><td>1.83</td><td>2.04</td></tr></tbody></table>

Concerning local Newton iterations, which correspond to the extra cost of nonlinear localization techniques, we see that adaption makes it possible to reduce their number, with slightly smaller figures in Tables 4 than in Table 2, and more significant improvement in Table 5 for mixed approaches.

## 7.2 Discussion

Assessments proved that the nonlinear localization methods were able to reduce the number of Krylov iterations (and then the number of communication) at the price of some extra local computations (factorization), without modifying the number of global (outer Newton) iterations. Unfortunately our `octave-mpi` implementation does not allow reliable time measurements and large simulations. Moreover, the time performance of the methods would be highly dependent on the hardware and in particular on the performance of the network compared to the performance of the cores.

We chose the previous experiment in order to prove that performance is highly dependent on the chosen boundary condition. In particular Neumann boundary conditions for perfectly plastic domains lead to very stiff problems, because in that case plasticity is estimated from the above (plasticity is overestimated and Newton iterations to reduce it are conducted on very unrealistic configurations). During our experimentations, we were able to design configurations with positive hardening where the dual approach was equivalent or even better than the primal.

More generally, it is easy to build test cases where one variant behaves poorly, or even fails. Figure 3 presents two such cases. In the first case, the decomposition makes subdomains very slender and prone to buckling (assuming large transformations). In that case, the method may converge (slowly) to a suboptimal solution [21]. In the second case, the decomposition amplifies the initial deformation, and assuming damageable behavior, subdomains may fail even if
the global problem was well-posed.

On the contrary, any experiment where the chosen boundary condition is barely influenced by the nonlinearity will lead to better parallel performance than a linear problem: an extreme case is an isostatic elastoplastic lattice solved with a dual approach.

<figure><img src="image_4.png" /><figcaption>Figure 3: Two critical cases</figcaption></figure>

We can thus draw the following conclusions:

* The method may not converge or follow an inadequate path if local instabilities are possible. Based on the literature of Schwarz methods for nonlinear problems [1, 31], we expect the method to converge unconditionally (assuming small enough increments) in the case of coercive maximal monotone operators (typically positive hardening in small strains) where the existence and interesting properties of nonlinear Schur complements can be established [42].

* The shape of the subdomains shall play an important role in the potential development of nonlinear effects.

* The initial boundary condition is critical, zero may be a very bad guess which triggers unduly complex local nonlinear computations. It is worth considering starting the method by a global elastic prediction.

* Well chosen Robin conditions should help follow the good path (the one of the monolithic solution). In the example above, we could compute the tangent Schur complement of the remainder of the structure which is (first order) optimal, but it is of pure academic interest, since the associated numerical cost would not be acceptable for actual industrial simulations; approximations of that operator is a classical question [23, 15] which we will address in future works. The approximation studied in the example, constituted by the stiffness of the neighbors, is not flexible enough, and it is not expected to give good results when the remainder of the structure is very compliant.

# 8 Conclusion

In this article, we have proposed nonlinear versions of FETI, BDD and FETI2LM algorithms where nonlinear computations are carried out independently on the subdomains. The starting point is a “nonlinear condensation” of the problems which leads, after linearization, to the parallel evaluation of the residual by nonlinear problems with boundary conditions characteristic of the method and to the solution to a classical DD system using any classical linear solver amongst FETI, BDD or FETI2NL.

The method modifies the nonlinear system to be solved and introduces inner inexact solvers, which implies a careful tuning of the stopping criteria. We proposed to exploit a classical inexact Newton formula.

A first assessment proved the potential of the method with an interesting trade-off between Krylov iterations (which corresponds to network communications) and Local Newton iterations (associated with local CPU computations).

However the method suffers from several drawbacks. First, the performance depends on the chosen boundary condition, which is a strong difference with the linear case. Second, the shape of the subdomains plays an even greater role than in the linear case, since it can trigger unphysical nonlinear effects. Third, in the case of localized nonlinearity, the method makes it possible to concentrate computations on the nonlinear subdomains but linear subdomains are then inactive, leading to poor load balancing (but potential power saving when idling processors enter standby mode). The question of load balancing will be addressed in future work when a more capable implementation is available. Typically we believe that it will be possible to propose an efficient strategy by encompassing the nonlinear problem in a mesh/substructuring adaption loop when a guaranteed precision is required by the user [38].

Much work is thus still needed before the method is reliable. We believe that well-chosen Robin conditions are an important ingredient to make the method stable and efficient.
References

[1] L. Badea. On the Schwarz alternating method with more than two subdomains for nonlinear monotone problems. *SIAM J. Numer. Anal.*, 28(1):179–204, 1991.

[2] Felipe Bordeu, Pierre-Alain Boucard, and Pierre Gosselet. Balancing domain decomposition with nonlinear relocalization: parallel implementation for laminates. In *Proceedings of the first international conference on parallel, distributed and grid computing for engineering*, pages CCP: 90, paper : 4, Pécs, 2009.

[3] D. Brands, A. Klawonn, O. Rheinbach, and J. Schröder. Modelling and convergence in arterial wall simulations using a parallel FETI solution strategy. *Comput. Meth. Appl. Mech. Engrg.*, 11(5):569–583, October 2008.

[4] Xiao-Chuan Cai and David E. Keyes. Nonlinearly preconditioned inexact Newton algorithms. *SIAM J. Sci. Comput.*, 24(1):183–200, 2002.

[5] P.G. Ciarlet. *Linear and nonlinear functional analysis with applications*. SIAM, 2013.

[6] P. Cresta, O. Allix, C. Rey, and S. Guinard. Nonlinear localization strategies for domain decomposition methods: Application to post-buckling analyses. *Computer Methods in Applied Mechanics and Engineering*, 196:1436–1446, 2007.

[7] Ron S Dembo, Stanley C Eisenstat, and Trond Steihaug. Inexact Newton methods. *SIAM Journal on Numerical analysis*, 19(2):400–408, 1982.

[8] Maksymilian Dryja and Wolfgang Hackbusch. On the nonlinear domain decomposition method. *BIT*, 37(2):296–311, 1997.

[9] Stanley C Eisenstat and Homer F Walker. Choosing the forcing terms in an inexact Newton method. *SIAM Journal on Scientific Computing*, 17(1):16–32, 1996.

[10] C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, and D. Rixen. FETI-DP: a dual-primal unified FETI method - part i: a faster alternative to the two-level FETI method. *International Journal for Numerical Methods in Engineering*, 50(7):1523–1544, 2001.

[11] C. Farhat, K. Pierson, and M. Lesoine. The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems. *Computer Methods in Applied Mechanics and Engineering*, 184(2-4):333–374, 2000.

[12] C. Farhat and F. X. Roux. The dual Schur complement method with well-posed local Neumann problems. *Contemporary Mathematics*, 157:193, 1994.

[13] C. Farhat and F. X. Roux. Implicit parallel processing in structural mechanics. *Computational Mechanics Advances*, 2(1):1–124, 1994. North-Holland.

[14] M. Feistauer and A. Zenisek. Finite element solution of nonlinear elliptic problems. *Numerische Mathematik*, 50:451–475, 1987.

[15] Lionel Gendre, Olivier Allix, and Pierre Gosselet. A two-scale approximation of the Schur complement and its use for non-intrusive coupling. *International Journal for Numerical Methods in Engineering*, 87(9):889–905, 2011.

[16] Pierre Gosselet. Domain decomposition methods for nonlinear problems in structural mechanics. Seminars in mathematical faculty of Essen (Germany), lab. N. Oresme of Caen, lab. MSSMat of école centrale de Paris, lab. GeM of école centrale Nantes (France); slides available on hal.archives-ouvertes.fr, 2011.

[17] Pierre Gosselet and Christian Rey. Non-overlapping domain decomposition methods in structural mechanics. *Archives of computational methods in engineering*, 13(4):515–572, 2007.

[18] Pierre Gosselet, Christian Rey, and Julien Pebrel. Total and selective reuse of Krylov subspaces for the solution to a sequence of nonlinear structural problems. *International Journal for Numerical Methods in Engineering*, 94(1):60–83, 2013.

[19] Pierre Gosselet, Daniel Rixen, François-Xavier Roux, and Nicole Spillane. Simultaneous-FETI and Block-FETI: robust domain decomposition with multiple search directions. *International Journal for Numerical Methods in Engineering*, 104(10):905–927, 2015.

[20] D. Hauer. The p-Dirichlet-to-Neumann operator with applications to elliptic and parabolic problems. *Journal of differential equations*, 259(8):3615–3655, 2015.
[21] J. Hinojosa, O. Allix, P.A. Guidault, and P. Cresta. Domain decomposition methods with nonlinear localization for the buckling and post-buckling analyses of large structures. *Advances in Engineering Software*, 70:13–24, 2014.

[22] Feng-Nan Hwang and Xiao-Chuan Cai. A class of parallel two-level nonlinear Schwarz preconditioned inexact Newton algorithms. *Computer Methods in Applied Mechanics and Engineering*, 196(8):1603–1611, 2007.

[23] C. Japhet, F. Nataf, and F. Rogier. The optimized order 2 method. application to convection-diffusion problems. *Future Generation Computer Systems*, 18(1):17–30, 2001.

[24] C. Kelley. *Solving Nonlinear Equations with Newton's Method*. Society for Industrial and Applied Mathematics, 2003.

[25] Axel Klawonn, Martin Lanser, and Oliver Rheinbach. Nonlinear FETI-DP and BDDC methods. *SIAM J. Sci. Comput.*, 36(2):A737–A765, 2014.

[26] Axel Klawonn and Oliver Rheinbach. Robust FETI-DP methods for heterogeneous three dimensional elasticity problems. *Computer Methods in Applied Mechanics and Engineering* and *Computer Methods in Applied Mechanics and Engineering*, 196(8):1400–1414, 2007.

[27] Axel Klawonn, Oliver Rheinbach, and Olof B. Widlund. An analysis of a FETI-DP algorithm on irregular subdomains in the plane. *SIAM J. Numer. Anal.*, 46(5):2484–2504, 2008.

[28] Axel Klawonn and Olof Widlund. FETI and Neumann-Neumann iterative substructuring methods: Connections and new results. *Communications on Pure and Applied Mathematics*, 54(1):57–90, 2001.

[29] P. Ladevèze. Sur une famille d'algorithmes en mécanique des structures. *Comptes Rendus Académie des Sciences - Mécanique*, Paris, 300(2):41–45, 1985.

[30] P. Ladevèze. *Nonlinear Computational Structural Mechanics – New Approaches and Non-Incremental Methods of Calculation*. Springer Verlag, 1999.

[31] P. Ladevèze, D. Néron, and P. Gosselet. On a mixed and multiscale domain decomposition method. *Computer Methods in Applied Mechanics and Engineering*, 196(8):1526–1540, 2007.

[32] Patrick Le Tallec. Domain decomposition methods in computational mechanics. *Comput. Mech. Adv.*, 1(2):121–220, 1994.

[33] Jan Mandel. Balancing domain decomposition. *Communications in Numerical Methods in Engineering*, 9(3):233, 1993.

[34] Augustin Parret-Fréaud, Christian Rey, Pierre Gosselet, and Frédéric Feyel. Fast estimation of discretization error for FE problems solved by domain decomposition. *Computer Methods in Applied Mechanics and Engineering*, 199(49-52):3315–3323, 2010.

[35] J. Pebrel, Christian Rey, and P. Gosselet. A nonlinear dual-domain decomposition method: Application to structural problems with damage. *International Journal for Multiscale Computational Engineering*, 6(3):251, 2008.

[36] Julien Pebrel, Pierre Gosselet, and Christian Rey. Etude du choix des conditions d'interface pour des stratégies non linéaire de décomposition de domaine. In *Actes du neuvième colloque national en calcul des structures*, volume 2, pages 393–398, Giens (Var), 2009.

[37] Christian Rey and Franck Risler. A Rayleigh-Ritz preconditioner for the iterative solution to large scale nonlinear problems. *Numerical Algorithms*, 17(3/4):279–311, 1998.

[38] Valentine Rey, Christian Rey, and Pierre Gosselet. A strict error bound with separated contributions of the discretization and of the iterative solver in non-overlapping domain decomposition methods. *Computer Methods in Applied Mechanics and Engineering*, 270(1):293–303, 2014.

[39] F. Risler and C. Rey. Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale non-linear problems. *Numerical algorithms*, 23:1, 2000.

[40] Daniel J. Rixen and Charbel Farhat. A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems. *International Journal for Numerical Methods in Engineering*, 44(4):489–516, 1999.

[41] François-Xavier Roux. A feti-2lm method for non-matching grids. In Michel Bercovier, Martin J. Gander, Ralf Kornhuber, and Olof Widlund, editors, *Domain Decomposition Methods in Science and Engineering XVIII*, volume 70 of *Lecture Notes in Computational Science and Engineering*, pages 121–128. Springer Berlin Heidelberg, 2009.
[42] R.E. Showalter. *Hilbert space methods for partial differential equations*. Monographs and studies in mathematics. Pitman, London, 1977.

[43] R.E. Showalter. *Monotone operators in Banach space and nonlinear partial differential equations*, volume 49 of *Mathematical surveys and monographs*. American mathematical society, 1997.

[44] N. Spillane, V. Dolean, P. Hauret, F. Nataf, C. Pechstein, and R. Scheichl. Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps. *Numer. Math.*, 126(4):741–770, 2014.