Dealing with nonlinear phenomena has become one of the predominant issues for mechanical engineers, in the objective of virtual testing. Whether they are geometrical or related to the material behavior, nonlinearities can be treated by a combination of Newton and linear solvers. Newton algorithms can be modified, secant, quasi-Newton [1, 2, 3, 4]
, depending mostly on the complexity of tangent operators computation. If the meshed structure has a large number of degrees of freedom, linear solvers are chosen to be iterative and parallel, belonging to the class of Domain Decomposition Methods for instance[5, 6, 7, 8, 9].
This article focuses on the nonlinear substructuring and condensation method, which has been investigated in previous studies [10, 11, 12, 13]. The substructured formulation involves a choice of interface transmission conditions, which can be either primal, dual or mixed, referring either to interface displacements, nodal interface reactions, or a linear combination of the two previous types – i.e. Robin interface conditions. In this context, the mixed formulation has shown good efficiency [14, 15, 16, 13], mostly due to a sound choice of the parameter introduced in the linear combination of interface conditions. Being homogeneous to a stiffness, and often refered to as an interface impedance, this parameter can indeed be optimized, depending on the mechanical problem . However, the computational cost of the optimal value involves in general storage and manipulation of global matrices, and is consequently not affordable in the framework of parallel computations.
The interface impedance, in DDM methods for structural mechanics, should model, from the point of view of one substructure, its interactions with the complement of the whole structure. In order to achieve good convergence rates without degrading computational speed, interface impedance can generally be approximated either by short scale or long scale formulations, depending on the predominant phenomena which must be accounted for. In the mechanical context, for instance, a common short scale approximation can be built by assembling interface stiffness of the neighbors [14, 15, 16, 13].
However, filtering long range interactions gives quite a coarse approximation of interface impedance, and does not give an accurate representation of the environment of each substructure. A good evaluation of the remainder of the structure should indeed couple these two strategies. Starting from this consideration, we propose here a new construction process of the interface impedance, based on a “spring in series” modeling of the structure, which couples the long and short range interactions with the structure. The heuristic we develop is strongly influenced by the availability of the various terms involved in our approximation.
The first section of this paper introduces the reference (mechanical or thermal) problem and the notations used in the following. A succinct presentation of the nonlinear substructuring and condensation method then recalls the principles of its mixed formulation: how the interface nonlinear condensed problem is built from nonlinear local equilibriums, and the basics of the whole solving process, involving a global Newton algorithm combined with two internal solvers (parallel local Newton algorithms and a multi-scale linear preconditioned Krylov solver for tangent interface system). At Section 4, the question of finding a relevant Robin parameter for mixed interface conditions is developed, mainly based on the observation that for each substructure, the optimal interface impedance is the nonlinear discretized Dirichlet-to-Neumann operator of its complementary part. The new heuristic is then introduced, starting from the model of two springs in series, and a possible nonlinear multi-scale interpretation is given. The details of the two-scale approximation can be found at subsection 4.5, its efficiency is evaluated at last section on several academic numerical examples.
2 Reference problem, notations
2.1 Global nonlinear problem
We consider here a nonlinear partial differential equation on a domain, representative of a structural mechanical or thermal problem, with Dirichlet conditions on the part of its boundary, and Neumann conditions on the complementary part . After discretization with the Finite Element method, the problem to be solved reads:
Vector takes into account boundary conditions (Dirichlet or Neumann) and dead loads, operator refers to the discretization of homogeneous partial differential equation.
In linear elasticity, under the small perturbations hypothesis, one has:
with the stiffness matrix of the structure.
Classical DDM notations will be used – see figure 6: global domain is partitioned into subdomains . For each subdomain, a trace operator restricts local quantities defined on to boundary quantities defined on :
Quantities defined on internal nodes (belonging to ) are written with subscript : .
Global primal (resp. dual) interface are noted (resp ). Primal assembly operators are defined as canonical prolongation operators from to : is a full-ranked boolean matrix of size - where is the size of global primal interface and the number of boundary degrees of freedom belonging to subdomain .
Diamond notations are used in the following: for a domain substructured in subdomains , concatenated local variables are superscripted , or , depending on the alignment.
Any matrix satisfying can be assigned to dual assembly operator – see figure 6 for the most classical choice.
3 Nonlinear substructuring and condensation: mixed formulation
This section recalls the principle of nonlinear substructuring and condensation, which is explained in details in .
3.1 Formulation of the condensed problem
Nonlinear problem (1) is decomposed into nonlinear subproblems:
where is the unknown local interface nodal reaction, introduced to represent interactions of the subdomain with neighboring subdomains.
Transmission conditions hold:
The mixed formulation consists in introducing a new interface unknown:
where the matrix is a parameter of the method. It has to be symmetric positive definite, and can be interpreted as a stiffness added to the interface, per subdomain: is called interface impedance.
Local equilibriums can then be reformulated as:
We assume the existence, at least locally, of a nonlinear mixed analogue of the Schur complement (ie. a discrete Robin-to-Dirichlet operator):
The tangent operator to can be explicitly computed in function of the tangent stiffness :
Moreover, in the linear case, the Robin-to-Dirichlet operator written is affine, with the constant term associated with external forces:
For the upcoming discussion, we will make use of the nonlinear primal Schur complement (Dirichlet-to-Neumann, noted ) which is such that . The tangent primal Schur complement can be computed from the tangent stiffness matrix:
and we have . Note that the tangent dual Schur complement (Neumann-to-Dirichlet) can be written as . In the linear case, the primal Schur complement is an affine operator with the constant term due to the external load:
Thanks to the complementarity between balanced and continuous quantities, and to the symmetry positive definiteness of , any boundary displacement (defined independently on neighboring subdomains) can be split in a unique way into a continuous field belonging to and a balanced field belonging to . Thus, the transmission conditions can be written in terms of and , and gathered in a single equation:
Finally, interface condensed problem reads:
3.2 Solving strategy
3.2.1 Newton-Krylov algorithm
Local solutions of nonlinear equilibriums (2) are computed by applying local Newton algorithms.
The interface mixed residual is assembled.
The interface tangent problem is solved by a DDM solver.
Newton global algorithm can be written, with previous notations:
Tangent problem then reads:
3.2.2 Alternative formulation
Tangent problem (5) could be treated by a FETI-2LM solver . An equivalent formulation of problem (4) is also possible, where the boundary interface unknown is replaced by a couple of interface unknowns , being a nodal reaction and an interface displacement. Couple is made unique by imposing the three following conditions:
With this formulation, tangent problem is expressed by:
3.2.3 Typical algorithm
Algorithm LABEL:alg:robin-bdd sums up the main steps of the method with the mixed nonlinear local problems and primal tangent solver. For simplicity reasons, only one load increment was considered.
As can be seen in this algorithm, several convergence thresholds are needed:
Global convergence criterion : since our approach is mixed, the criterion not only controls the quality of the subdomains balance (as in a standard Newton approach) but also the continuity of the interface displacement which is measured by an appropriate norm written .
Local nonlinear thresholds , which are associated with the Newton processes carried out independently on subdomains.
The global linear threshold of the domain decomposition (Krylov) solver (here BDD).
The other parameters of the method are the initializations of the various iterative solvers and the choice of the impedance matrices .
4 New heuristic for the interface impedance
The parameter is involved all along the solving, and a special care should be paid to its computation.
In order to frame the ideas, let us consider the Robin-Robin algorithm with stationary iteration, in the nonlinear case, with nonlinear impedance. Starting from the initial guess , we have the iterations of Algorithm LABEL:alg:robin-robin.
Assembled quantities and are defined such that the interface conditions can be written as:
and we assume the nonlinear local operators to be such that the equivalence between (8) and the following equation is ensured:
Considering a given subdomain , and writing its complement, we can condense the whole problem on its interface; the boundary displacement must then be the solution to:
Comparing this equation with line (1) of algorithm LABEL:alg:robin-robin, one can see that, starting from a zero initial guess , the method converges in only one iteration with : the ideal impedance is the Dirichlet-to-Neumann operator of the complement.
In order to further discuss the problem, we now consider the linear case, and we recall that we have: . In that case, the optimal impedance is thus an affine operator whose linear part (which we will write in agreement with the development of previous section) accounts for the stiffness of the complement domain, whereas the constant part accounts for the external load on the complement part. Note that another point of view is to use a strictly linear impedance together with a good (non-zero) initialization for which should account for the external load on the complement domain.
The construction of a good constant part for the impedance is usually realized, in the linear case, by the introduction of a well-chosen coarse problem; this is discussed in subsection 4.4. In the nonlinear case, building a coarse problem which would connect all subdomains during their inner Newton loop seems complex; more, it would break the independent computations. It looks simpler to rely on a good initialization in order to propagate the right-hand side: this can be done at low cost by easing accuracy constraints in the first inner Newton loop (adapting ), and then using the multiscale solver of the global linear step. Note that in  a coarse problem is built for nonlinear versions of FETIDP and BDDC but, again, it mainly serves to find a good initialization before independent parallel nonlinear solves.
We now focus on the construction of , i.e. the linear part of the impedance. In the linear case, one can show [20, 21] that for a slab-wise decomposition of the structure (or a tree-like decomposition, i.e. whose connectivity graph has no cycles), the setting is optimal, in the sense that the convergence is reached in a maximum number of iterations equal to the number of subdomains (iterations are only needed to propagate the right-hand side). If the convenient coarse grid is added, convergence can be extremely fast. For an arbitrary decomposition, the optimality of such a setting can theoretically be lost, because of the unclear propagation rate of the right-hand side [22, 21]. However, the pertinence of this value still seems to be ongoing, especially being given the difficulty to define a more relevant setting for a matrix operator .
Starting from these considerations, let us further analyze the terms of the following expression of the interface impedance for a given subdomain :
The first term, , accounts for very local interactions. It is sparse, and exactly has the fill-in of matrix . The second term, , accounts for long range interactions, it depends on the whole structure (geometry and material), and couples all degrees of freedom together via in-depth interactions. It is thus a full matrix; this property can be seen as the consequence of the pseudo-differentiability of the underlying Steklov-Poincaré operator of which the Schur complement is the discretization. It is important to note the minus sign: the short range part is very stiff and the global effects mitigate it.
Obviously, formula (10) is intractable in a distributed environment. However, different strategies have been investigated to compute approximations at low cost – see next subsection for a quick review.
In the nonlinear context, the use of a linear impedance is of course non-optimal. Moreover, the best linear impedance probably resembles the Schur complement of the remainder of the subdomain in the final configuration, which is of course unknown a priori. Our aim is then to try to find a heuristic which gives an easy-to-compute approximation of the formula (10) to be applied to the initial tangent stiffness.
4.2 Quick review
The question of finding a good approximation of the Schur complement of a domain is at the core of mixed domain decomposition methods like optimized Schwarz methods  or the Latin method . Studies have proved that they needed to reproduce short-range effects (like local heterogeneity) but also structural effects (like the anisotropy induced by the slenderness of plate structures 26]. A possibility in order to better model short-range interaction between interface nodes is to use Ventcell conditions instead of simple Robin conditions ; this enables to recover the same sparsity for the impedance as for the stiffness of the subdomain. An extreme strategy is to use (scalar) Robin conditions on the Riesz’ image of the normal flux leading to a fully populated impedance matrix . A more reasonable strategy is to use a strip approximation of the Schur complement , which can also be computed by adding elements to the subdomains , in the spirit of restricted additive Schwarz methods .
From an algebraic point of view, short range approximation (or even ) is sometimes used for FETI’s preconditioner , where it is called lumped approximation. Let be the set of the neighbors of subdomain , we have
Being an assembly among a few subdomains of sparse block-diagonal matrices, this term is quite cheap to compute, and does not require any extra-computations, since local tangent stiffnesses are calculated anyway at each iteration of the solving process. The efficiency of the simple approximation (11) has been studied, in the context of nonlinear substructuring and condensation, in some research works [10, 16, 13], and has given good results when tested on rather homogeneous structures of standard shape.
In the domain decomposition framework for linear problems, long range interactions are taken into account thanks to the coarse grid problems [32, 5, 33], which enables the method to comply with Saint-Venant’s principle. These are closely related to projection techniques inspired by homogenization [34, 35, 24, 36, 37] in order to get low rank approximations. Let be an orthonormal basis of a well chosen subspace of displacements, the approximation can be written as:
Saint-Venant’s principle imposes to contain at least the rigid body motions of , for computational efficiency it can be complemented by affine deformation modes or by displacements defined independently by interfaces.
However, if short range approximations do not provide enough information to give a good representation of the faraway structure influence on a substructure , neither do long range approximation give a good estimation of the near-field response to a sollicitation. Besides, in the context of small displacements, a lack of precision on the close structure is more problematic than the filtering of long range interactions: predominant mechanical reactions usually come from nearby elements of the mesh.
The best strategy for would combine both short and long range formulations, however this version has not been much investigated yet. In particular it is not that easy to ensure the positivity of the impedance if the two approximations are computed independently. In  an expensive scale separation was introduced in the context of non intrusive global/local computations where was somehow available (which is not the case in our distributed framework). We propose here a new expression for parameter , in the context of nonlinear substructuring and condensation with mixed interface conditions, which combines short and long scale formulations, at low computational cost.
4.3 Spring in series model
Our heuristic for the impedance relies on the simple observation that finding a two-scale approximation of the flexibility of may be more patent than for the stiffness. It is inspired by the simple model of two springs assembled in series: one spring models the stiffness of the neighboring subdomains whereas the second models the stiffness of the faraway subdomains (see figure 7). The resulting equivalent flexibility is the sum of the two flexibilities.
In practice, in order to recover the structure of (10), while remaining tractable, we propose the local flexibility to be the inverse of a sparse matrix, and the long-range flexibility to be low-rank. The latter condition is also motivated by [39, 40], where it is shown that low-rank approximants of fully populated inverse operators, arising from FE discretization of elliptic problems, can be derived from the hierarchical-matrices theory. Typically we have:
where can refer for instance to expressions (11) or (12), is a small-sized square matrix, and an interface vectors basis of size . Writing the local contribution of basis , expression (14) can be inversed using the Sherman-Morrisson formula:
This stiffness is a sparse matrix corrected by a low-rank term; then, when solving the (generalized) Robin problems, the Sherman-Morrisson formula can be used again:
The short-range term enables to regularize the problem without impairing the sparsity of the stiffness matrix.
4.4 A multi-scale interpretation
Starting from Algorithm LABEL:alg:robin-robin, a macroscopic condition, inspired from the Latin method, can be imposed on the nodal reactions: the nodal reactions should satisfy a weak form of the interface balance, defined by a macroscopic basis :
After algebraic calculations, writing local equilibriums with this new condition leads to:
where is a projector on the low-dimension subspace .
Not only the coarse space associated to the macroscopic constraint (17) results in the propagation of the right-hand side on the whole structure ( is not sparse) but also in the modification of the impedance by the symmetric negative low rank term .
In our nonlinear context, considering the basic setting (11) or (12), the modification proposed in (14) can be seen as the introduction of a multi-scale computation inside the mixed nonlinear substructuring and condensation method. As said earlier, the propagation of the right-hand side is ensured by a well-built initialization, which can be realized by adapting the inner Newton criterion at each global iteration.
4.5 Two-scale approximation of the flexibility
4.5.1 General idea
From previous analysis, we try to derive an approximation of the (linear) optimal flexibility (10) which takes the additive form of (14). Being given a substructure , we write the assembly of local tangent Schur complements on the remainder .
Using the quotient and the inverse formulas for the Schur complement, we have:
We here assume a substructuring ensuring the inversibility of and , i.e. Dirichlet conditions are not concentrated on only one subdomain, and the complementary part of each subdomain is connected. In practice, this is almost always the case; if not, a simple subdivision can overcome the problem.
Classical preconditioners of BDD-algorithm can then be used as approximations of the inverse of . We hence introduce the concatenation of the scaled local traces of rigid body motions () of subdomains belonging to , with scaled assembly operators taking into account the absence of matter inside subdomain . Considering the classical definition of scaled assembly operators , modified operators can be defined as:
Let be the -orthogonal projector on :
The BDD-theory states that in the first term, can be conveniently approximated by a scaled sum of local inversesThe GENEO theory  states that, if needed, computable extra modes shall be inserted in in order to maintain the quality of the approximation.. After developing and factorizing, we have a first approximation of the flexibility:
4.5.2 Long range interactions term
The second term of expression (20), written , is a matrix of low rank , where is the number of neighbors rigid body motions. It could be used as is, however its computation involves the inversion of quantity , an interface matrix of rank , where is the number of local rigid body modes of the whole remainder . In the context of large structures with a high number of subdomains, can increase drastically; saving the computation and factorization of such a matrix could then become quite interesting. Moreover, during the computation of the structure coarse problem, a close quantity is already assembled and factorized: the matrix – with and . Compared to , the addition of the local term linked to in somewhat balances the classical scaling on its boundary (taking into account non-existant matter inside ), we thus propose:
4.5.3 Short range interactions term
The first term of expression (20), written , can also be simplified. First, for numerical efficiency, a diagonal lumping technique is used to approximate the local Schur complements (as explained in section 4.2). Then, in order to preserve sparsity, the projectors are removed. Assuming stiffness scaling is used we then directly recover the inverse of the superlumped stiffness of the neighbors:
4.5.4 Scaling issue
A way to avoid building the modified scaled assembly operators is to notice that for , the following relation holds between modified and classical scaling operators :
and we observe that the local diagonal matrix can be extracted without cost from :
With evident notations, for a scaling based on the material stiffness, the diagonal coefficient of associated with degree of freedom is equal to:
Final expression. To conclude, we propose the following two-scale impedance:
4.6 Attempt to enrich the short-range approximation
The short range part of the impedance, corresponding to the sparse approximation of by , seems very crude. In particular, we most probably underestimate the flexibility of the neighbors by using a diagonal operator.
We believe it is worth mentioning the tentative improvement which consisted in adding another low rank term:
where is a diagonal matrix and an orthonormal basis, approximations of the eigen-elements of associated with the higher part of the spectrum. They could be obtained at a moderate cost by post-processing the tangent BDD iterations in the spirit of 
(but considering the classical eigenvalues instead of the generalized ones).
This low rank term could be concatenated with the one associated with rigid body motions , and thus did not modify the usability of the approximation. We observed that it led to a stiffness which was closer to our reference (measured with the Frobenius norm). But in practice when using it as the impedance in our numerical experiments, the reduction achieved in iterations numbers was not worth the additional cost of the enrichment term – this is why we do not present it in detail. This “improvement” may be more useful on other classes of nonlinear problems for which it would be important not to overestimate the stiffness of the remainder of the structure.
5.1 Two test cases
The efficiency of the expression (22) is evaluated on two numerical test cases. First test case is a bi-material beam with bending load, represented on figure 8. Material and geometrical parameters are given in table 1: one of the two materials is chosen to be elastoplastic with linear hardening, the other one is chosen to remain elastic. Load is applied with imposed displacement on the edge defined by .
5.2 Elastic analysis
The ultimate goal of this paper is to assess the performance of the new impedance (22) in the nonlinear multi-scale distributed context. Before we reach that point, a preliminary mono-scale elastic study is performed in order to verify that the heuristic developed in previous sections is actually able to capture both short and long range interactions within the structure.
Sollicitations are here keeped low enough to remain in the elastic domain of every materials: bi-material beam and multiperforated beam are both submitted to a bending load of intensity . More, decomposition is for now only performed along -axis (multiple points will be involved in next section, where the nonlinear multi-scale context is considered). One of the interest of the elastic linear case with slab-wise decomposition relies on the ability to express the optimal interface impedance: (see 4.1). Even if the computational cost of this parameter would be, in a real situation, absolutely not affordable in the context of parallel resolutions, it was calculated here for the purpose of our analysis. A comparison with an optimal reference can thus be made for the two following expressions:
Being given the alternative formulation we chose for the mixed nonlinear substructuring and condensation method (see section 3.2.2), an elastic resolution would be strictly equivalent to a primal BDD resolution. Therefore, no comparison of different interface impedances is possible with Algorithm LABEL:alg:robin-bdd. A mono-scale FETI-2LM solver  was hence implemented, corresponding to the first formulation of the mixed interface problem with the unknown (5). This algorithm enables to solve linear problems with Robin interface transmission conditions.
Note that an optimal coarse problem could be added in order to recover an efficient multi-scale solver [45, 46, 47]. However, this augmentation strategy would make it impossible to discern the efficiency of the long range interactions term of our two-scale impedance. Again, our aim is not to compete with augmented Krylov solvers for linear problems but to find an alternative way, compatible with nonlinear problems, to introduce long-range effects. The mono-scale formulation is thus preserved in order to evaluate the ability of (22) to introduce in local equilibriums information related to the interactions with the far structure, in a linear context where the optimal parameter is known.
Results are given on table 4 for the two previously introduced test cases.
As expected, for both test cases, the optimal interface impedance