A domain decomposition method employing displacement-only partitioned equations for quasi-static structural analysis

Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 Available online 5 August 2024 0045-7825/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma A domain decomposition method employing displacement-only partitioned equations for quasi-static structural analysis Seung-Hoon Kanga, K.C. Parkb, José A. Gonzálezc, SangJoon Shina,d,∗ aDepartment of Aerospace Engineering, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, 08826, Republic of Korea bAerospace Mechanics Research Center, Smead Aerospace Engineering Sciences, University of Colorado Boulder, Boulder, 80309-429, CO, USA cEscuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, Sevilla, 41092, Spain dInstitute of Advanced Aerospace Technology, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, 08826, Republic of Korea ARTICLE INFO Keywords: Parallel computation Displacement-only partitioned formulation Static structural analysis Domain decomposition method ABSTRACT The present study reports a family of iterative domain decomposition method for the static structural analysis, labeled as AFETI-DO and its variants (AFETI-DOc and AFETI-DOci), all of which employ a recently developed Displacement-Only (DO) partitioned formulation (Park et al., 2023). The DO partitioned equation (𝑑(𝐟−𝐊𝐝)= 0) consists of the applied force 𝐟, the block-diagonal stiffness matrix (𝐊) for each partition, the coupling projection operator (𝑑), and the unknown partitioned displacements (𝐝). Numerical experiments with the AFETI-DO and its two variants demonstrate that, while their performance is comparable to existing iterative solution algorithms for homogeneous problems, for highly heterogeneous two and threedimensional problems the proposed methods yield an accurate convergence threshold compared with the FETI-DP method (Farhat et al., 2001). For fourth-order elasticity problems (i.e., thin plates and shells), the present methods offer a comparable accuracy and computational effort required with the FETI-DP method. In addition, while not extensive, the numerical performance of the proposed AFETI-DO method and its variants appear to exhibit scalability and improved accuracy, provided the preconditioners are paired with the commensurate heterogeneities. 1. Introduction Parallel computation is nowadays one of the common strategies for the solution of large structural finite element (FE) models. Varieties of parallel algebraic solvers [1,2] have been implemented into the matrix representation of FE analysis, in conjunction with both shared and distributed memory architectures. The multigrid method is one of the most efficient parallel solvers. In the case of purely algebraic solvers, the algebraic multigrid [3] will stand out, which was recently enhanced by machine learning techniques [4] to increase its performance further. Another equally popular approach is the use of iterative algorithm based on domain decomposition or partitioning methods [5]. In this, modification of the governing partial differential equation itself with additional constraints has been attempted. The partitioning method separates the globally assembled equation into the ones for the local subdomains and interface, and conventionally solves for an interfacial unknown. This solution process is suitable for distributed memory architectures while providing significant scalability for large-scale analysis. In particular, this paper deals with a non-overlapping partitioning method. A host of methods of this type have been proposed in the literature, and one of the most remarkable progresses in the structural FE analysis was shown by the FE tearing and interconnecting (FETI) method [6]. The original one-level FETI (FETI-1) method enforces the interface continuity via dual Lagrange ∗Corresponding author at: Institute of Advanced Aerospace Technology, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, 08826, Republic of Korea. E-mail address: [email protected] (S. Shin). https://doi.org/10.1016/j.cma.2024.117273 Received 29 February 2024; Received in revised form 28 July 2024; Accepted 28 July 2024 Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 2 S.-H. Kang et al. multipliers, which is the fundamental variable to be solved during the iteration. To satisfy the self-equilibrium state of the floating subdomains, the partitioned solution corresponding to rigid body modes has been used to construct the so-called coarse space. Farhat et al. [7] determined that such coarse space plays a crucial role in the global imbalance propagation during the iteration. The balancing domain decomposition (BDD) method is another representative partitioning method proposed by Mandel [8]. It satisfies the interfacial equilibrium using the interfacial displacement, which is a primal variable. The theoretical and numerical similarities between FETI and BDD methods have been extensively investigated [9,10]. To date, numerous variants have been proposed based on these two methods. The algebraically partitioned FETI (AFETI) method [11,12] employed the localized Lagrange multipliers at the subdomain interface, which provided a unique rank-sufficient constraint definition without duplicities. Its similarity to the FETI-1 method was reported with respect to the numerical and theoretical aspects [13–15]. The two-level FETI (FETI-2) method [16,17] additionally imposed the subdomain corner continuity, to improve convergence for the fourth-order plate and shell elements. Beyond that, Farhat et al. [18,19] proposed the dual-primal FETI (FETI-DP) method, in which a coarse space was defined by the primal displacement constraint at the subdomain corners. It did not introduce rigid body modes, and universal usage for static and dynamic cases was accomplished with a fast convergence speed. For three-dimensional elasticity problems, scalability has been obtained by augmentation of the coarse space with auxiliary modes [18]. Klawonn et al. [20,21] proposed the adaptive coarse space for highly heterogeneous case. Furthermore, the recursive multi-level approach of the FETI-DP method [22] has exhibited scalability for static analyses up to 21 billion degrees of freedom (DOFs). On the other hand, the BDD method with the corner constraints was proposed by Dohrmann [23]. Furthermore, its variants emerged with various approaches such as the inexact preconditioner [24,25], multi-space, multi-level [26], adaptive [27], and mixed approaches [28]. Recently, a novel displacement-only (DO) partitioned formulation [29] (labeled as the DP formulation in that paper) has been introduced. In this formulation, the unknown vector comprises partitioned displacements for the structural problem, each corresponding to a subdomain. In essence, the DO (displacement-only) partitioned equations are analogous to the assembled finiteelement equations in that the only unknowns of the DO partitioned equations are partitioned displacements. This is illustrated below: Assembled equations: 𝐌𝑔 𝐝𝑔+𝐊𝑔𝐝𝑔=𝐟𝑔,(𝐌𝑔,𝐊𝑔)=𝐋⊤ 𝑔(𝐌,𝐊)𝐋𝑔 DO partitioned equations: 𝐌 𝐝+𝑑𝐊𝐝 =𝑑𝐟 Kinematic relations: 𝐝=𝐋𝑔𝐝𝑔 Herein, subscript 𝑔refers to assembled (or global) system, 𝐋𝑔is the assembly Boolean operator, (𝐌,𝐊)are partitioned block diagonal mass and stiffness matrices, (𝐝,𝐝𝑔)are partitioned and assembled displacements, (𝐟,𝐟𝑔)are partitioned and assembled force vectors, and 𝑑is the coupling projector operator. Notations  □and □⊤denote the second-order temporal differentiation and transpose of □, respectively. This is in contrast to the use of saddle-point equations whose unknowns consist of dual variables, i.e., subdomain displacements and the interface Lagrange multipliers, in most domain decomposition methods. The DO partitioned equations adopt the unassembled block-diagonal stiffness matrix in which each subdomain maintains its intact characteristics independently. Thus, in addition to offering implementation simplicity and computational efficiency, such properties of the DO formulation may provide robust accuracy for the heterogeneous problems. By exploiting the aforementioned desirable properties of the DO partitioned formulation, the objective of the present paper is to develop iterative static solution algorithms for the static structural analysis, labeled as AFETI-DO method (and its two variants AFETI-DOc and AFETI-DOci). The rest of the paper is organized as follows. Section 2offers a brief review of the DO partitioned formulation that is applicable to structural dynamical systems. Section 3adapts the DO partitioned formulation to be suitable for the solution of static structural problems by modifying the coupling projection operator and illustrates that the static DO partitioned equations yield the desired correct static displacements. AFETI-DO method is presented in Section 4, followed by the development of several preconditioners that play a key role in accelerating iterative convergence in Section 5. The numerical performance of the AFETI-DO method and its variants is reported in Section 6, focusing on simple and yet representative two-dimensional and three-dimensional heterogeneous elasticity and plate problems. The example problems illustrate that the present AFETI-DO method offers perhaps the simplest implementation among domain decomposition methods and yields robust accuracy for heterogeneous problems whose accurate solutions have challenged some of the existing domain decomposition methods that utilize saddle-point formulations. 2. Variational framework and the DO partitioned formulation Before introducing the static DO partitioned formulation, this section will briefly explain the variational framework based on localized Lagrange multipliers [30,31] and the second-order dynamic DO partitioned formulation [29]. Fig. 1 illustrates the situation. The assembled FE discretization is partitioned into multiple subdomains, and each subdomain interfacial node is connected to a frame node. The subdomain interface problem consists of the interface nodes and its associated Lagrange multipliers. Because such approach localizes the Lagrange multipliers to each subdomain, the interfacial connectivity will be uniquely defined even at the subdomain corners. As shown in Eq. (1), the total virtual work of the partitioned system 𝛿𝑊𝑡becomes the sum of the virtual work of subdomains 𝛿𝑊 and that of the interfacial constraints 𝛿𝑊𝑐. 𝛿𝑊𝑡=𝛿𝑊 +𝛿𝑊𝑐(1) Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 3 S.-H. Kang et al. Fig. 1. Partitioning of an assembled FE discretization into 4 subdomains using localized Lagrange multipliers. For a linear elastic undamped system, 𝛿𝑊 will be expressed as follows: 𝛿𝑊 (𝐝)= 𝑁𝑠 ∑ 𝑠=1 𝛿𝑊 𝑠(𝐝𝑠)= 𝑁𝑠 ∑ 𝑠=1 𝛿𝐝𝑠⊤{𝐌𝑠 𝐝𝑠+𝐊𝑠𝐝𝑠−𝐟𝑠} = 𝛿𝐝⊤{𝐌 𝐝+𝐊𝐝 −𝐟},(2) 𝐝=⎧ ⎪ ⎨ ⎪ ⎩ 𝐝1 ⋮ 𝐝𝑁𝑠 ⎫ ⎪ ⎬ ⎪ ⎭ ,𝐟=⎧ ⎪ ⎨ ⎪ ⎩ 𝐟1 ⋮ 𝐟𝑁𝑠 ⎫ ⎪ ⎬ ⎪ ⎭ ,𝐌=⎡⎢⎢⎣ 𝐌1⋯𝟎 ⋮ ⋱ ⋮ 𝟎⋯𝐌𝑁𝑠⎤⎥⎥⎦ ,𝐊=⎡⎢⎢⎣ 𝐊1⋯𝟎 ⋮ ⋱ ⋮ 𝟎⋯𝐊𝑁𝑠⎤⎥⎥⎦ .(3) Herein, 𝛿𝑊 𝑠is the virtual work of the 𝑠th subdomain, and 𝑁𝑠is the total number of the subdomains. Vectors 𝐝𝑠and 𝐟𝑠are the displacement and external force of the 𝑠th subdomain, respectively. Matrices 𝐌𝑠and 𝐊𝑠are the mass and stiffness of the 𝑠th subdomain. Hereafter, 𝐌will be assumed to be symmetric positive definite (SPD). Notation 𝛿□denotes the first-order variation of □. Meanwhile, the displacement compatibility relationship at the subdomain interface will be written as below: ⎧ ⎪ ⎨ ⎪ ⎩ 𝐁1⊤𝐝1−𝐋1 𝑓𝐝𝑓 ⋮ 𝐁𝑁𝑠⊤𝐝𝑁𝑠−𝐋𝑁𝑠 𝑓𝐝𝑓 ⎫ ⎪ ⎬ ⎪ ⎭ =𝐁⊤𝐝−𝐋𝑓𝐝𝑓=𝟎,(4) 𝐁=⎡⎢⎢⎣ 𝐁1⋯𝟎 ⋮ ⋱ ⋮ 𝟎⋯𝐁𝑁𝑠⎤⎥⎥⎦ ,𝐋𝑓=⎡⎢⎢⎢⎣ 𝐋1 𝑓 ⋮ 𝐋𝑁𝑠 𝑓 ⎤⎥⎥⎥⎦ .(5) Herein, 𝐝𝑓is the frame displacement vector. And, 𝐁𝑠and 𝐋𝑠 𝑓are the interfacial and frame Boolean matrices, in which 𝐁𝑠⊤𝐝𝑠and 𝐋𝑠 𝑓𝐝𝑓become the interfacial displacement vector and its associated frame displacement vector for the 𝑠th subdomain. Matrix 𝐁𝑠maps the interfacial DOFs to the subdomains ones in the 𝑠th subdomain, while 𝐋𝑠 𝑓extracts the interfacial DOFs of the 𝑠th subdomain from the frame DOFs. Note that in the present variational framework based on localized Lagrange multipliers [30,31], both matrices 𝐁𝑠 and 𝐋𝑠 𝑓are unsigned, boundary comprised of only 1 and 0, due to the intermediate frame between the interfaces. In contrast, using a conventional Lagrange multiplier [6,16,18] will require the signed interfacial Boolean matrix that consists of 1, −1, and 0. Detailed comparison between localized and conventional Lagrange multipliers is found in Ref. [13]. Then, the Lagrange multiplier vector 𝝀will be introduced for the interfacial energy functional. Herein, 𝝀⊤=[𝝀1⊤⋯𝝀𝑁𝑠⊤], in which 𝝀𝑠is the Lagrange multiplier vector connected to the 𝑠th subdomain. The virtual work introduced by the interfacial constraint will be written as in Eq. (6). 𝛿𝑊𝑐(𝐝,𝐝𝑓,𝝀)=𝛿(𝝀⊤{𝐁⊤𝐝−𝐋𝑓𝐝𝑓})(6) After substituting Eqs. (2) and (6) into Eq. (1), the partitioned equation is derived from the stationary-point condition of the total virtual work as ⎡⎢⎢⎢⎣ 𝐊+𝐌𝑑2 𝑑𝑡2𝐁 𝟎 𝐁⊤𝟎−𝐋𝑓 𝟎−𝐋⊤ 𝑓𝟎 ⎤⎥⎥⎥⎦ ⎧ ⎪ ⎨ ⎪ ⎩ 𝐝 𝝀 𝐝𝑓 ⎫ ⎪ ⎬ ⎪ ⎭ =⎧ ⎪ ⎨ ⎪ ⎩ 𝐟 𝟎 𝟎⎫ ⎪ ⎬ ⎪ ⎭ ,(7) where 𝑑2∕𝑑𝑡2is the second-order temporal differentiation operator. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 4 S.-H. Kang et al. Introducing the nodal discrete accelerations vector, the DO partitioned dynamic equation [29] is expressed as 𝐌 𝐝=𝑑{𝐟−𝐊𝐝},(8) where 𝑑=𝐈−𝐁𝑓𝐌𝑏𝐁⊤𝐌−1,𝑓=𝐈−𝐌𝑏𝐋𝑓𝐌−1 𝑓𝐋⊤ 𝑓,𝐌𝑏=[𝐁⊤𝐌−1𝐁]−1 ,𝐌𝑓=𝐋⊤ 𝑓𝐌𝑏𝐋𝑓.(9) Herein, matrix 𝑑is a projector enforcing the interface equilibrium. The interfacial and frame mass matrices, 𝐌𝑏and 𝐌𝑓, are SPD because the partitioned mass matrix 𝐌is previously assumed to be SPD. In addition, 𝐈is the identity matrix. As shown in Eq. (8), the DO formulation is described by only the partitioned displacement vector 𝐝, while the dual and primal interfacial variables 𝝀and 𝐝𝑓are eliminated. Vector 𝐌 𝐝is the partitioned inertial force, and 𝐟−𝐊𝐝 is the uncoupled partitioned sum of the external and internal forces. By using the idempotent projector 𝑑, the term 𝑑{𝐟−𝐊𝐝}will become the partitioned sum of the external, internal, and interfacial forces, which will satisfy the self-equilibrium condition in conjunction with the partitioned inertia 𝐌 𝐝. A more detailed explanation and derivation on the DO formulation is provided in Ref. [29]. Remark 1. DO formulation treats the interfacial displacement differently from the primal method (e.g., BDD [8]). The primal method solves for the assembled interfacial displacement, which corresponds to the frame displacement 𝐝𝑓(see, Fig. 1 and Eq. (4)). Meanwhile, vector 𝐝in DO formulation (Eq. (8)) consists of all the partitioned displacement. Thus, its interfacial components, 𝐁1⊤𝐝1 to 𝐁𝑁𝑠⊤𝐝𝑁𝑠, are defined independently along the subdomain interface. 3. Static DO partitioned formulation In this section, the DO formulation is particularized to statics. In the static DO partitioned formulation, the partitioned acceleration (  𝐝) will vanish from Eq. (8), as follows: 𝑑𝐊𝐝 =𝑑𝐟,(10) with the mass matrix remaining in the definition of the projector in Eq. (9). A classical adaptation of solving the preceding quasi-static equilibrium equation (Eq. (10)) is obtained as the following dynamic relaxation form:  𝐌 𝐝+ 𝐂 𝐝+𝑑𝐊𝐝 =𝑑𝐟( 𝐝, 𝐝)→𝟎 ⟹𝑑{𝐟−𝐊𝐝}→𝟎,(11) where  𝐌and  𝐂are the pseudo-mass and pseudo-damping matrices, respectively; and the projector 𝑑is constructed using the pseudo-mass matrix  𝐌. When developing the proposed partitioned solution algorithm for statics, named as AFETI-DO method, we will adopt the pseudomass concept from classical relaxation techniques (see, e.g., Ref. [32]), while employing the preconditioned conjugate gradient (PCG) technique in lieu of relaxation pseudo-time integration with artificial damping term (  𝐂 𝐝). Before the AFETI-DO method is described in the next section, three prerequisite aspects will be addressed. First, the construction of the projector 𝑑employing a pseudo-mass matrix will be discussed. Second, the interfacial compatibility equation for Eq. (10) will be introduced. Finally, it will be shown that a direct solution of Eq. (10) yields the same solution obtained by the assembled FE equations so that an iterative solution of Eq. (10) also yields the correct solution as well. 3.1. Projector 𝑑for the static DO partitioned formulation Although the partitioned inertial force is not considered in Eq. (10), the partitioned mass matrix 𝐌will still remain in the projector 𝑑. As shown in Eq. (9),𝑑consists of the mass and Boolean coupling matrices (𝐌,𝐋𝑓, and 𝐁). Of several pseudo-mass matrices  𝐌that were considered, the following diagonal stiffness will be selected: 𝐌← 𝐌=𝐃,𝐃= diag (𝐊).(12) Such selection will provide the stiffness regularization effect for the iterative solution algorithm, as found in Ref. [29]. In addition, introducing a diagonal matrix will significantly simplify the constitution of projector 𝑑. This projector, in which 𝐌←𝐃is substituted into 𝑑in Eq. (9), will be from now on denoted as  𝑑. For a clearer explanation, the static DO partitioned formulation in Eq. (10) will be re-expressed in the form as follows:  𝑑𝐊𝐝 = 𝑑𝐟,(13) meaning that 𝑓,𝐌𝑏, and 𝐌𝑓in Eq. (9) will be re-written as follows: 𝑓← 𝑓=𝐈−𝐃𝑏𝐋𝑓𝐃−1 𝑓𝐋⊤ 𝑓,𝐌𝑏←𝐃𝑏=[𝐁⊤𝐃−1𝐁]−1 =𝐁⊤𝐃𝐁,𝐌𝑓←𝐃𝑓=𝐋⊤ 𝑓𝐃𝑏𝐋𝑓.(14) Here, 𝐃𝑏is a diagonal matrix in which the entries are the diagonal stiffness of the interfacial DOFs. 𝐃𝑓is also a diagonal matrix, in which each entry is the sum of the interfacial diagonal stiffness sharing the same frame DOF. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 5 S.-H. Kang et al. Because the diagonal matrix 𝐃−1 has no coupling entry between the interfacial and internal DOFs, the property 𝐁⊤𝐃−1 =𝐃−1 𝑏𝐁⊤ is fulfilled. Using this property, the projector  𝑑can be expressed as 𝑑← 𝑑=𝐈−𝐁 𝑓𝐃𝑏𝐁⊤𝐃−1 =𝐈−𝐁 𝑓𝐃𝑏𝐃−1 𝑏𝐁⊤=𝐈−𝐁 𝑓𝐁⊤,(15) where 𝐈− 𝑑is enclosed by the matrices 𝐁and 𝐁⊤. Thus, the projector  𝑑will operate only on the interfacial DOFs through the frame projector  𝑓. Example 1. As the unassembled DO equation, Eqs. (10) and (13), may be relatively brand new, an example of 4 bar elements partitioned into 4 subdomains will be introduced as shown in Fig. 2. In the figure, 𝑒and 𝑠are the elemental and subdomain indices. 𝐸(𝑒),𝐴(𝑒), and 𝐿(𝑒)denote the elastic modulus, cross-sectional area, and length of 𝑒th element, respectively. The assembled displacement 𝐝𝑔, partitioned displacement 𝐝, frame displacement 𝐝𝑓, and the associated Lagrange multiplier 𝝀will be expressed as 𝐝𝑔=⎧ ⎪ ⎨ ⎪ ⎩ 𝑑𝑔,1 𝑑𝑔,2 𝑑𝑔,3 ⎫ ⎪ ⎬ ⎪ ⎭ ,𝐝=𝐋𝑔𝐝𝑔= ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝑑1 1 𝑑2 1 𝑑3 1 𝑑3 2 𝑑4 1 𝑑4 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ ,𝐝𝑓={𝑑𝑓},𝝀=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜆1 1 𝜆2 1 𝜆3 1 𝜆4 1 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ ,(16) where 𝑑𝑔,𝑗 is the displacement of 𝑗th DOF of the assembled FE. 𝑑𝑠 𝑗is partitioned displacement of the 𝑗th DOF in terms of the 𝑠th subdomain. 𝑑𝑓is the frame displacement that connects the partitioned displacements 𝑑1 1,𝑑2 1,𝑑3 1, and 𝑑4 1.𝜆𝑠 𝑗is the Lagrange multiplier connected to the 𝑗th DOF of 𝑠th subdomain. In Eq. (16),𝐋𝑔is an unsigned assembly Boolean matrix that relates the assembled displacement to the partitioned displacement. Once, 𝐋𝑔is constructed by a standard FE routine, the partitioned interface Boolean matrix 𝐁and the frame Boolean matrix 𝐋𝑓will be determined in a usual way (see, Appendix A in Ref. [33] that computes 𝐁and 𝐋𝑓). For the example shown in Fig. 2, the Boolean matrices will be determined as in Eq. (17), relating the displacement vectors to be 𝐝=𝐋⊤ 𝑔𝐝𝑔and 𝐁⊤𝐝=𝐋𝑓𝐝𝑓. 𝐋𝑔=⎡⎢⎢⎢⎢⎢⎢⎣ 100 100 100 010 100 001 ⎤⎥⎥⎥⎥⎥⎥⎦ ,𝐁=⎡⎢⎢⎢⎢⎢⎢⎣ 1000 0100 0010 0000 0001 0000 ⎤⎥⎥⎥⎥⎥⎥⎦ ,𝐋𝑓=⎡⎢⎢⎢⎢⎣ 1 1 1 1 ⎤⎥⎥⎥⎥⎦ .(17) Given the partitioned stiffness matrix 𝐊and the assembly Boolean matrix 𝐋𝑔are computed via a suitable FE software, 𝐌,𝐌𝑏, and 𝐌𝑓will be obtained as follows: 𝐊= blockdiag (𝐤fix,4𝐤fix,3𝐤free,2𝐤free),𝐤fix =[1],𝐤free =[1 −1 −1 1 ], 𝐌←𝐃= diag (𝐊)= blockdiag (1,4,3,3,2,2), 𝐌𝑏←𝐃𝑏=𝐁⊤𝐃𝐁 = blockdiag (𝐷1 1, 𝐷2 1, 𝐷3 1, 𝐷4 1)= blockdiag (1,4,3,2), 𝐌𝑓←𝐃𝑓=𝐋⊤ 𝑓𝐃𝑏𝐋𝑓=[𝐷𝑓]=[10]. (18) where blockdiag is meant to form block-diagonal matrix from the entries, scalars, or submatrices. Then, the frame projector  𝑓will be obtained as below: 𝐈− 𝑓=𝐃𝑏𝐋𝑓𝐃−1 𝑓𝐋⊤ 𝑓=1 𝐷𝑓 ⎡⎢⎢⎢⎢⎢⎣ 𝐷1 1𝐷1 1𝐷1 1𝐷1 1 𝐷2 1𝐷2 1𝐷2 1𝐷2 1 𝐷3 1𝐷3 1𝐷3 1𝐷3 1 𝐷4 1𝐷4 1𝐷4 1𝐷4 1 ⎤⎥⎥⎥⎥⎥⎦ =1 10 ⎡⎢⎢⎢⎢⎣ 1111 4444 3333 2222 ⎤⎥⎥⎥⎥⎦ .(19) As shown, by extracting the diagonal interfacial terms from the partitioned stiffness matrix, the entries of  𝑑will be obtained without conducting the complicated matrix multiplication in the original projector 𝑑. Therefore,  𝑑does not have to be explicitly constructed, which will demonstrate the implementation simplicity of the present iterative algorithm. In addition, the stiffness regularization effect will be expected from the projector  𝑓. Fig. 3 illustrates the computation involved in the matrix–vector multiplication with the projectors  𝑑and  ⊤ 𝑑, based on the example of 4 bar elements in Fig. 2 and its related values in Eqs. (16) to (19). In it, 𝐱and 𝐲are the input and projected output vectors, entries of which are notated as same as those of vector 𝐝in Eq. (16). In advance to the projection, the preparation stage will extract the interfacial diagonal stiffness for each subdomain. Subsequently, for each frame DOF, the frame diagonal stiffness 𝐷𝑓will be obtained by summing the associated interfacial diagonal stiffness. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 6 S.-H. Kang et al. Fig. 2. Example of 4 bar elements partitioned into 4 subdomains. The DOF numbering is localized for each subdomain in the partitioned system. Fig. 3. Procedure of the matrix–vector multiplication with projectors  𝑑and  ⊤ 𝑑for the example of 4 bar elements partitioned into 4 subdomains. Circled alphabets indicate the sequence of operations. The matrix–vector multiplication by the projector  𝑑begins by summing the associated interfacial components of input vector 𝐱for each frame DOF. Then, each interfacial DOF obtains the associated frame summation 𝑥𝑓and computes the output through the diagonal stiffness scaling, i.e., 𝑦𝑠 𝑗=(𝐷𝑠 𝑗∕𝐷𝑓)𝑥𝑓. For internal DOF, the output is the same as the input. Multiplication with the projector  ⊤ 𝑑is similar to that with  𝑑, except that the scaling sequence is inverted. The scaling, 𝑥𝑠 𝑗=(𝐷𝑠 𝑗∕𝐷𝑓)𝑥𝑠 𝑗, is performed first for each interfacial DOF. Subsequently, each frame DOF sums the associated interfacial DOFs 𝑥𝑠 𝑗, which will become the output. Same as the projector  𝑑, no operation is required for the internal DOFs. 3.2. Interfacial compatibility for the static DO partitioned formulation The static DO partitioned formulation in Eq. (13) cannot solely produce a correct displacement solution. The displacement will have a deficiency regarding the interfacial compatibility, which is the second row in Eq. (7). Instead, constraints are imposed on the acceleration. Such condition will not be valid for the static analysis, which is not equipped with the time integration scheme. Therefore, the interfacial compatibility will be additionally introduced to the static DO partitioned formulation in Eq. (13). Referring to the kinematic relation between partitioned and assembled displacements, 𝐝=𝐋𝑔𝐝𝑔, the interfacial compatibility becomes 𝐝∈(𝐋𝑔), in which the (□)indicates the range space of □. Recalling Eqs. (16) and (17), assembling Boolean operator 𝐋𝑔is a full rank matrix and contains column vector with multiple nonzero entries. Such column vector maps single assembled displacement to multiple partitioned displacements while retaining their values to be the same among the subdomains. In this Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 7 S.-H. Kang et al. Table 1 Direct solution of the displacement vector for the example of 4 bar elements partitioned into 4 subdomains. The function lsqr of MATLAB R2023b is used for the partitioned system. Assembled FE DO formulation Eq. (13) Eqs. (13) and (20) 𝐝𝑔=𝐊−1 𝑔𝐟𝑔𝐝=[ 𝑑𝐊]+ 𝑑𝐟 𝐝 = ⊤ 𝑑[ 𝑑𝐊 ⊤ 𝑑]+ 𝑑𝐟 𝐝assembled 𝑔=⎧ ⎪ ⎨ ⎪ ⎩ 1 2 2⎫ ⎪ ⎬ ⎪ ⎭ 𝐝= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0.29 1.18 −0.50 0.50 −0.50 0.50 ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ ,𝐝partitioned 𝑔=⎧ ⎪ ⎨ ⎪ ⎩ 0.12 0.50 0.50⎫ ⎪ ⎬ ⎪ ⎭ 𝐝= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 1 1 2 1 2 ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ ,𝐝partitioned 𝑔=⎧ ⎪ ⎨ ⎪ ⎩ 1 2 2⎫ ⎪ ⎬ ⎪ ⎭ paper, the interfacial compatibility will be enforced by utilizing the useful property of projector  𝑑,(𝐋𝑔)=( ⊤ 𝑑), derivation of which is included in Appendix A.1. Then, the partitioned displacement vector will be assumed as 𝐝= ⊤ 𝑑 𝐝,(20) which symmetrizes the matrix of the static DO partitioned formulation in Eq. (13). Note that notation  𝐝, which appears here and hereafter, does not contain physical interpretation. It is employed as an unknown vector for the iterative algorithm, and becomes a partitioned displacement vector 𝐝after it is multiplied by projectors. 3.3. Solutions by static DO and assembled FE equations are identical Example 2. The example of 4 bar elements in Fig. 2 will now be revisited to validate whether solutions by the static DO and the assembled FE equations are identical. To this end, we note that the assembled equations will be given by 𝐊𝑔𝐝𝑔=𝐟𝑔,𝐊𝑔=𝐋⊤ 𝑔𝐊𝑔𝐋𝑔=⎡⎢⎢⎣ 10 −3 −2 −3 3 0 −2 0 2 ⎤⎥⎥⎦ ,𝐟𝑔=⎧ ⎪ ⎨ ⎪ ⎩ 0 3 2⎫ ⎪ ⎬ ⎪ ⎭ ,𝐝assembled 𝑔=𝐊−1 𝑔𝐟𝑔=⎧ ⎪ ⎨ ⎪ ⎩ 1 2 2⎫ ⎪ ⎬ ⎪ ⎭ .(21) The static DO partitioned equations will be obtained as  𝑑𝐊𝐝 = 𝑑𝐟,𝐝= ⊤ 𝑑 𝐝, 𝑑=𝐈−𝐁 𝑓𝐁⊤,  𝑓=𝐈−𝐃𝑏𝐋𝑓𝐃−1 𝑓𝐋⊤ 𝑓=1 10 ⎡⎢⎢⎢⎢⎣ 9 −1 −1 −1 −4 6 −4 −4 −3 −3 7 −3 −2 −2 −2 8 ⎤⎥⎥⎥⎥⎦ ,𝐟=𝐋𝑔[𝐋⊤ 𝑔𝐋𝑔]−1 𝐟𝑔= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0 0 0 3 0 2 ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ ,𝐝= ⊤ 𝑑[ 𝑑𝐊 ⊤ 𝑑]+ 𝑑𝐟= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 1 1 2 1 2 ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ .(22) where □+indicates Moore–Penrose pseudo-inverse of □. Finally, assembled solution from the partitioned solution 𝐝will be as follows: 𝐝partitioned 𝑔=[𝐋⊤ 𝑔𝐋𝑔]−1 𝐋⊤ 𝑔=⎧ ⎪ ⎨ ⎪ ⎩ 1 2 2⎫ ⎪ ⎬ ⎪ ⎭ .(23) The results are also summarized in Table 1. As shown, Eq. (20) (𝐝= ⊤ 𝑑 𝐝) provides a correct solution which is the same as those of assembled FE. In addition, Appendix A.2 provides a demonstration of the way the projector ⊤ 𝑑, a general form of  ⊤ 𝑑, enforces the vector 𝐝to satisfy the interface compatibility condition. Given such condition, the projector 𝑑will operate on the partitioned system correctly, while conserving the energy to be the same as that in the assembled one. The static DO formulation retains an identical solution to the assembled FE for the case of an iterative algorithm. Fig. 4 shows that the PCG algorithm based on DO formulation, AFETI-DO, provides the displacement solution that converges to that of the assembled FE. It also describes the efficacy of diagonal stiffness 𝐃employed in the projector  𝑑. Due to the stiffness regularization effect, 𝐌←𝐃provides superior convergence performance to 𝐌←𝐌physical, where 𝐌physical is a physical structural mass with unit density. As shown in Fig. 4, the effect of stiffness regularization is retained for the increased internal DOFs. The proof-of-concept study by Park et al. [29] also supported the correctness of the iterative DO solution and the effect of 𝐃 via the example of a simple but heterogeneous plate bending element. Detailed implementation of the AFETI-DO algorithm will be explained in the next sections. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 8 S.-H. Kang et al. Fig. 4. Convergence of the iterative DO solution for the bar element example. Discretization in Fig. 2 is used as a baseline mesh. 𝐊+is selected as the preconditioner. Remark 2. The rigid body modes are not introduced in the present static DO partitioned formulation, which is the combination of Eqs. (13) and (20). This is in contrast to the traditional partitioning methods such as FETI-1 and one-level AFETI (AFETI-1) methods [6,11,12], in which the rigid body modes and their amplitudes will be naturally introduced and constitute the coarse space. Nevertheless, the rigid body modes still provide useful properties, which will be elaborated in the next sections. 4. Present AFETI-DO algorithm Based on the aforementioned characteristics, this section introduces an iterative algorithm for the present static DO partitioned formulation. Also, several preconditioning matrices will be presented. We have labeled the ensuing iterative algorithm as AFETI-DO. The acronym AFETI is inherited from the Algebraically partitioned Finite Element Tearing and Interconnecting method [11,12] that employs the method of localized Lagrange multipliers [13,30,31]. The label DO signifies the Displacement-Only characteristics as alluded to in Introduction. 4.1. Preconditioned conjugate gradient (PCG) algorithm Considering that the traditional partitioning methods [6,8,11,16,19,23] usually employ the PCG algorithm to obtain the solution, this paper will also employ PCG algorithm for the present static AFETI-DO method. In addition to Eqs. (13) and (20), the partitioned residual vector 𝐫will be introduced for the iterative solution algorithm, as follows: 𝐫= 𝑑{𝐟−𝐊𝐝},𝐝∈( ⊤ 𝑑).(24) Solving the above equation will be equivalent to obtaining the solution  𝐝in Eq. (25). 𝐫= 𝑑𝐟− 𝑑𝐊 ⊤ 𝑑 𝐝(25) Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 9 S.-H. Kang et al. Because the partitioned stiffness matrix 𝐊is symmetric positive semi-definite (SPSD), the matrix  𝑑𝐊 ⊤ 𝑑will also be SPSD. The vector 𝐟is determined to be 𝐟=𝐋⊤ 𝑔 +𝐟𝑔=𝐋𝑔[𝐋⊤ 𝑔𝐋𝑔]−1 𝐟𝑔in the partitioning method, which means that the external nodal force in the assembled system will be equally divided to the subdomains if it is applied to the interface. Because 𝐟∈(𝐋𝑔), the vector 𝐟is not included in the null space of 𝐊, which consists of the partitioned rigid body modes. Therefore, the vector  𝑑𝐟is included in the range space of  𝑑𝐊 ⊤ 𝑑, and then the PCG algorithm will be applicable [34]. Using the idempotent characteristics of the projector  𝑑, the PCG algorithm for the present static AFETI-DO method follows in Algorithm 1. Algorithm 1 AFETI-DO algorithm. 1: 𝐝0∈(𝐋𝑔),𝐟=𝐋+ 𝑔𝐟𝑔 2: 𝐫0= 𝑑{𝐟−𝐊𝐝0} 3: 𝐳0=𝐀+𝐫0 4: 𝐬0= ⊤ 𝑑𝐳0 5: 𝐩0=𝐬0 6: for 𝑘= 0,⋯, 𝑚 do 7: 𝐪𝑘=𝐊𝐩𝑘⊳Fully parallelizable stage 8: 𝛼𝑘=𝐫⊤ 𝑘𝐬𝑘 𝐩⊤ 𝑘𝐪𝑘 9: 𝐝𝑘+1 =𝐝𝑘+𝛼𝑘𝐩𝑘 10: 𝐫𝑘+1 =𝐫𝑘−𝛼𝑘 𝑑𝐪𝑘⊳Inter-processor stage 11: if the solution converges then 12: exit the loop 13: end if 14: 𝐳𝑘+1 =𝐀+𝐫𝑘+1 ⊳Fully parallelizable stage 15: 𝐬𝑘+1 = ⊤ 𝑑𝐳𝑘+1 ⊳Inter-processor stage 16: 𝛽𝑘=𝐫⊤ 𝑘+1𝐬𝑘+1 𝐫⊤ 𝑘𝐬𝑘 17: 𝐩𝑘+1 =𝐬𝑘+1 +𝛽𝑘𝐩𝑘 18: end for In Algorithm 1, the subscript 𝑘indicates the iteration count. The preconditioner is indicated as 𝐀+, in which a partitioned block-diagonal matrix will be employed as follows: 𝐀+=⎡⎢⎢⎣ 𝐀1+⋯𝟎 ⋮ ⋱ ⋮ 𝟎⋯𝐀𝑁𝑠+⎤⎥⎥⎦ .(26) As shown in Algorithm 1, the displacement vector 𝐝𝑘is enforced to be 𝐝𝑘∈( ⊤ 𝑑)=(𝐋𝑔)in each iteration. And, each iteration includes four matrix–vector multiplications. Two of them are multiplications with the block-diagonal partitioned matrices 𝐊and 𝐀+, which will be fully parallelizable for the case of the parallel computation. Remaining multiplications will occur with the projection matrices  𝑑and  ⊤ 𝑑. As inferred from Eq. (19), such multiplication requires the information of the adjacent subdomains for each subdomain, and becomes the inter-processor stage in the parallel computation. In the aspect of the distributed memory architecture, each processor stores the stiffness and preconditioning matrices, 𝐊𝑠and 𝐀𝑠+, of the assigned subdomains. In addition, it stores the diagonal stiffness information of the adjacent subdomains for the multiplication with the projection matrices  𝑑and  ⊤ 𝑑, which requires negligible storage when compared against the storage for 𝐊𝑠and 𝐀𝑠+. The initial solution vector 𝐝0is required to satisfy the interfacial compatibility condition 𝐝0∈(𝐋𝑔)in the AFETI-DO algorithm. Such condition will also be satisfied by multiplying an arbitrary vector  𝐝0by the projector  ⊤ 𝑑, i.e., 𝐝0= ⊤ 𝑑 𝐝0, using the property ( ⊤ 𝑑)=(𝐋𝑔). In the present paper, the initial solution vector for AFETI-DO algorithm is selected to be 𝐝0=𝟎among several candidates investigated. To determine whether the iterative algorithm terminates, it will be convenient to employ PCG residual vector 𝐫𝑘and observe its relative residual norm as 𝜖PCG 𝑘=‖𝐫𝑘‖2 ‖𝐫0‖2 .(27) Alternatively, the relative residual norm of the global force equilibrium might be employed, defined as 𝜖𝑔 𝑘=‖𝐟𝑔−𝐊𝑔𝐝𝑔,𝑘‖2 ‖𝐟𝑔‖2 ,(28) where 𝐝𝑔,𝑘 =𝐋+ 𝑔𝐝𝑘=[𝐋⊤ 𝑔𝐋𝑔]−1 𝐋⊤ 𝑔𝐝𝑘is the assembled displacement vector at 𝑘th iteration. Comparing Eqs. (27) and (28), the relative global force residual norm 𝜖𝑔 𝑘is regarded as a more rigorous error measure. As discussed in the study on FETI-1 method [6], the convergence of the relative PCG residual under a certain tolerance does not ensure Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 16 S.-H. Kang et al. Fig. 7. Plane stress example with element-by-element heterogeneity: Convergence history of global force residual for the case of 𝐻∕ℎ= 8 and 𝑁𝑠= 82. Fig. 8. Plane stress example with heterogeneous channel: Square plate for the case of 𝑁𝑠= 32. The gray- and black-colored elements are assigned with elastic modulus of 𝜒−0.5𝐸and 𝜒0.5𝐸, respectively. The corner point for FETI-DP is indicated by a circle. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 17 S.-H. Kang et al. Fig. 9. Plane stress example with heterogeneous channel: Convergence history for the case of 𝜒= 106,𝐻∕ℎ= 8, and 𝑁𝑠= 82. converges. In addition, Schur complement form of the traditional method contributes to the improved condition number and convergence rate. In contrast, PCG residual of the present DO formulation enforces the internal and interfacial equilibrium simultaneously during the iteration. Although such an expanded form of the residual vector is not favorable in the aspect of the condition number, it significantly enhances the solution accuracy. In Fig. 9, the global force residual norm of the deluxe FETI-DP method terminates to decrease at an early iteration, while its PCG residual norm gradually converges to a significantly lower value. However, the global force residual norm of the present methods decreases to the order from 10−6to 10−7, showing the most accurate solution despite the relatively large PCG residual norm. Meanwhile, low deformational modes might be employed for AFETI-DOc method to enhance the convergence rate. A similar approach is also found in Ref. [37]. Then, the initial vector 𝐝0and the coarse space projector  𝑐will be replaced to be 𝐝0← 𝐔𝑑[ 𝐔⊤ 𝑑𝐊 𝐔𝑑]−1  𝐔⊤ 𝑑𝐟, 𝑐←𝐈− 𝑑𝐊 𝐔𝑑[ 𝐔⊤ 𝑑𝐊 𝐔𝑑]−1 𝐔⊤,(55) where  𝐔𝑑= ⊤ 𝑑𝐑,𝐔=⎡⎢⎢⎣ [𝐑1|Φ1]⋯𝟎 ⋮ ⋱ ⋮ 𝟎⋯[𝐑𝑁𝑠|Φ𝑁𝑠]⎤⎥⎥⎦ ,(56) which modifies the Lines 1, 4, and 15 of the Algorithm 2in Appendix B. In Eq. (56),Φ𝑠is the matrix of low deformational modes of the 𝑠th subdomain, which is obtained from the following generalized eigenvalue analysis: Φ𝑠=[Φ𝑠 1,…,Φ𝑠 𝑁𝑠 𝛷],[𝐊𝑠−𝛽𝑠 𝑝𝐃𝑠]Φ𝑠 𝑝=𝟎, 𝛽𝑠 𝑝>0,(57) where 𝑠= 1,…, 𝑁𝑠and 𝑝= 1,…, 𝑁𝑠 𝛷.𝑁𝑠 𝛷is the number of low deformational modes for 𝑠th subdomain. Φ𝑠 𝑝and 𝛽𝑠 𝑝are 𝑝th deformational mode vector and corresponding eigenvalue in 𝑠th subdomain. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 18 S.-H. Kang et al. Table 7 Plane stress example with heterogeneous channel: Performance comparison among the partitioning methods for the case of 𝑁𝑠= 32,𝑁𝑟= 18, and 𝑁𝑐= 20. N/A indicates that the method does not reach the convergence criterion. 𝜒𝐻 ℎIteration count for 𝜖𝑔 𝑘<10−6 (condition number) AFETI-1 FETI-DP Deluxe FETI-DP AFETI-DO AFETI-DOc AFETI-DOci 10014 14 (1.6e1) 15 (4.9e0) 16 (7.3e0) 85 (2.3e3) 25 (6.7e1) 26 (7.1e1) 28 16 (2.0e1) 18 (7.9e0) 17 (1.0e1) 125 (4.7e3) 35 (1.4e2) 37 (1.5e2) 10214 48 (1.8e2) 50 (1.1e2) 25 (1.2e1) 396 (4.6e4) 62 (5.7e2) 61 (5.7e2) 28 N/A 57 (2.2e2) 27 (1.5e1) 572 (8.3e4) 81 (8.6e2) 79 (8.7e2) 10414 N/A 192 (1.1e4) 39 (3.7e1) 1377 (4.2e6) 275 (5.0e4) 229 (5.0e4) 28 N/A N/A N/A 2036 (7.6e6) 393 (7.4e4) 318 (7.4e4) 10514 N/A N/A N/A 2434 (4.2e7) 489 (5.0e5) 377 (5.0e5) 28 N/A N/A N/A 3599 (7.6e7) 1131 (9.7e5) 568 (7.4e5) 10614 N/A N/A N/A N/A 1224 (5.7e6) 606 (5.0e6) 28 N/A N/A N/A N/A N/A N/A Table 8 Plane stress example with heterogeneous channel: Performance of AFETI-DOc method with 6 deformation modes (𝑁𝑠 𝛷= 6) per subdomain for the case of 𝑁𝑠= 32. The number of coarse space DOFs is 𝑁𝑐= 20 for deluxe FETI-DP. Those for AFETI-DOc without and with deformation modes are 𝑁𝑟=𝑁𝑢= 18 and 𝑁𝑢= 72, respectively. N/A indicates that the method does not reach the convergence criterion. 𝜒𝐻 ℎIteration count for 𝜖𝑔 𝑘<10−6 (condition number) Deluxe FETI-DP AFETI-DOc (𝑁𝑠 𝛷= 0) AFETI-DOc (𝑁𝑠 𝛷= 6) 10014 16 (7.3e0) 25 (6.7e1) 14 (4.1e0) 28 17 (1.0e1) 35 (1.4e2) 20 (8.3e0) 10214 25 (1.2e1) 62 (5.7e2) 24 (1.9e1) 28 27 (1.5e1) 81 (8.6e2) 33 (3.0e1) 10414 39 (3.7e1) 275 (5.0e4) 40 (1.0e2) 28 N/A 393 (7.4e4) 53 (1.3e2) 10514 N/A 489 (5.0e5) 47 (1.1e2) 28 N/A 1131 (9.7e5) 76 (9.6e2) 10614 N/A 1224 (5.7e6) 105 (7.7e4) 28 N/A N/A 95 (3.6e3) The results are summarized in Table 8 and Fig. 10. Herein, 𝑁𝑢is the number of coarse space DOFs for such augmented AFETI-DOc method, which is a sum of the number of rigid body modes (𝑁𝑟) and deformational modes. As shown, adding 6 deformational modes per subdomain significantly improves the convergence behavior of AFETI-DOc, while still showing superior accuracy to FETI-DP. However, increasing the coarse space DOFs should be extensively examined and refined. For example, investigation is required for the trade-off between convergence improvement and computational cost, especially for large-scale analysis. In addition, other types of coarse mode preconditioning, e.g., balancing preconditioner [39], should be examined to explore the best strategy that is robust for the variety of numerical examples. Such study is one of the future works for the present approach. 6.1.3. Three-dimensional solid with element-by-element heterogeneity In Fig. 11, the example is extended to a three-dimensional situation and a larger number of DOFs. A cube of side length 𝑙= 1 m is discretized by eight-node hexahedral (HEX8) solid elements with 3 translational DOFs per node. One side (𝑋= 0) of the cube is clamped, and the other side (𝑋=𝑙) is under a uniformly distributed tensile load. Poisson’s ratio is 𝜈= 0.3, and the element-by-element material heterogeneity is considered. Elastic modulus is 𝐸= 100 MPa for the homogeneous example (𝜒= 1). For FETI-DP method, complemented Definition MD1 in Ref. [19] is selected as follows: corner points are defined to be located at the non-colinear points belonging to at least the three or more subdomains. Because FETI-DP method is not scalable for three-dimensional configuration, coarse space is augmented via the translational rigid body modes at the subdomain edges [18]. The results are summarized in Table 9 and Fig. 12. For the homogeneous material (𝜒= 1), the situation is similar to that in the previous two-dimensional examples. The present AFETI-DOc and AFETI-DOci methods show a rather larger number of iterations but show a comparable convergence rate against the AFETI-1 and FETI-DP methods. For the heterogeneous material (𝜒= 103and 106), the oscillatory convergence behavior of AFETI-1 and FETI-DP methods, shown in Fig. 7, is alleviated for the three-dimensional example. However, the solution accuracy of AFETI-1 is still inferior. FETIDP demonstrates better accuracy and convergence rate than AFETI-1, which may be due to the significantly increased number of coarse space DOFs and the continuity enforcement at the corner points. However, such a larger number of coarse DOFs results in increased computational cost, as depicted in Table 10. AFETI-1, FETI-DP, and the present AFETI-DOc and AFETI-DOci methods maintain the iteration count for the increased number of subdomains. The difference is that the initial residual of the present methods is significantly lower than that by AFETI-1 and FETI-DP methods for the case of the heterogeneity, and such initial accuracy Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 19 S.-H. Kang et al. Fig. 10. Plane stress example with heterogeneous channel: Global force residual history of AFETI-DOc with 6 deformational modes per subdomain (𝑁𝑠 𝛷) for the case of 𝜒= 106,𝐻∕ℎ= 8, and 𝑁𝑠= 82. Fig. 11. Three-dimensional solid with element-by-element heterogeneity: Cube for the case of 𝑁𝑠= 23. The gray- and black-colored elements are assigned with elastic modulus of 𝜒−0.5𝐸and 𝜒0.5𝐸, respectively. The corner point for FETI-DP is indicated by a sphere. Table 9 Three-dimensional solid with element-by-element heterogeneity: Performance comparison among the partitioning methods. N/A indicates that the method does not reach the convergence criterion. 𝜒𝐻 ℎ𝑁𝑠𝑁𝑟𝑁𝑐Iteration count for 𝜖𝑔 𝑘<10−6 (condition number) AFETI-1 FETI-DP AFETI-DO AFETI-DOc AFETI-DOci 100 443288 1,188 14 (1.1e1) 14 (5.2e0) 53 (2.4e2) 16 (1.8e1) 17 (2.0e1) 832688 12,348 22 (1.2e1) 16 (5.2e0) 114 (1.0e3) 27 (2.0e1) 28 (2.2e1) 843288 1,188 20 (1.8e1) 18 (7.5e0) 83 (5.7e2) 31 (4.6e1) 30 (4.8e1) 832688 12,348 29 (1.9e1) 20 (7.6e0) 174 (2.5e3) 40 (4.9e1) 39 (5.1e1) 1034832688 12,348 28 (1.5e1) 21 (1.4e1) 186 (1.6e3) 33 (2.3e1) 35 (2.5e1) 8832688 12,348 33 (2.1e1) 25 (1.3e1) 283 (4.6e3) 47 (5.1e1) 49 (5.3e1) 1064832688 12,348 N/A 22 (1.4e1) 187 (1.6e3) 33 (2.3e1) 35 (2.5e1) 8832688 12,348 N/A 26 (1.4e1) 284 (4.6e3) 47 (5.2e1) 49 (5.3e1) contributes to the computational efficiency. Furthermore, as shown in Table 10, the ratio of the fully parallelizable portion of the present methods ranges from 94% to 97%, which is higher than that of the FETI-DP (65%). This implies the advantages of the present methods for parallel architecture. In addition, all of the present methods exhibit superior solution accuracy to those of AFETI-1 and FETI-DP methods. The AFETI-DO method, which lacks the coarse space but is most efficient per iteration, ensures the solution accuracy regardless of the complexity of the situation to be solved. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 20 S.-H. Kang et al. Fig. 12. Three-dimensional solid with element-by-element heterogeneity: Convergence history of global force residual for the case of 𝐻∕ℎ= 4 and 𝑁𝑠= 43. Table 10 Three-dimensional solid with element-by-element heterogeneity: Computational time for the case of 𝜒= 103, 𝐻∕ℎ= 8, and 𝑁𝑠= 83. AFETI-1 FETI-DP AFETI-DO AFETI-DOc AFETI-DOci Iteration count 33 25 283 47 49 Computational time 55.8 s 115.9 s 70.6 s 61.2 s 52.9 s (fully parallelizable portion) 54.5 s 75.7 s 66.5 s 58.1 s 51.1 s 6.2. Fourth-order elasticity Fig. 13 illustrates the last example comprised of the discrete Kirchhoff triangular (DKT) plate element where 1 translational and 2 rotational DOFs are assigned per a node [42]. Such mixed DOFs of the fourth-order elasticity induce structural heterogeneity even for the homogeneous material. A plate of unit length is cantilevered at 𝑋= 0. In addition, the uniformly distributed 𝑍-directional Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 21 S.-H. Kang et al. Fig. 13. Fourth-order elastostatic example: Square plate for the case of 𝑁𝑠= 32. The gray- and black-colored elements are assigned with elastic modulus of 𝜒−0.5𝐸and 𝜒0.5𝐸, respectively. The corner point for FETI-DP is indicated by a circle. Table 11 Fourth-order elastostatic example: Performance comparison among the partitioning methods. N/A indicates that the method does not reach the convergence criterion. 𝜒𝐻 ℎ𝑁𝑠𝑁𝑟𝑁𝑐Iteration count for 𝜖𝑔 𝑘<10−6 (condition number) AFETI-1 FETI-DP AFETI-DO AFETI-DOc AFETI-DOci 100 84236 27 68 (3.3e2) 31 (1.7e1) 727 (2.6e6) 115 (1.2e3) 107 (1.2e3) 82168 147 N/A 39 (1.7e1) 3261 (5.0e7) 203 (1.4e3) 190 (1.4e3) 16 4236 27 124 (1.4e3) 39 (2.3e1) 1897 (2.1e7) 332 (1.0e4) 280 (9.9e3) 82168 147 N/A 50 (2.3e1) 9022 (4.2e8) 612 (1.4e4) 514 (1.1e4) 103882168 147 139 (3.3e2) 109 (1.5e2) 4711 (1.0e8) 248 (1.7e3) 228 (1.7e3) 16 82168 147 N/A 188 (7.6e2) 12,760 (8.6e8) 1173 (1.5e6) 629 (1.4e4) 106882168 147 144 (3.3e2) 145 (4.0e2) 4741 (1.0e8) 249 (1.7e3) 228 (1.7e3) 16 82168 147 N/A 269 (1.9e3) 12,826 (8.7e8) 790 (7.1e4) 627 (1.4e4) Table 12 Fourth-order elastostatic example: Computational time for the case of 𝜒= 106,𝐻∕ℎ= 16, and 𝑁𝑠= 82. FETI-DP AFETI-DO AFETI-DOc AFETI-DOci Iteration count 269 12,826 790 627 Computational time 3.47 s 63.5 s 7.72 s 5.25 s (fully parallelizable portion) 3.34 s 56.3 s 6.28 s 4.66 s force is imposed at the edge of 𝑋=𝑙. The thickness and Poisson’s ratio are 𝑡= 2 mm and 𝜈= 0.3. The element-by-element material heterogeneity is considered, and the elastic modulus is 𝐸= 100 MPa for the homogeneous case (𝜒= 1). For the FETI-DP method, Definition D2 in Ref. [19] is selected as follows: cross points are defined to be located at the points belonging to three or more subdomains. As shown in Table 11, the present AFETI-DOc and AFETI-DOci methods are not numerically scalable for the fourth-order elastostatic example. This induces the relatively degraded computation efficiency of the present method compared to the FETI-DP method, as shown in Table 12. Therefore, the future study will include improving the present methods for the fourth-order elasticity. It is found that the convergence threshold of AFETI-1 and FETI-DP methods is less sensitive to material heterogeneity, which is different from the previous two-dimensional second-order elastostatic example. The effect of heterogeneous DOFs is more dominant than that of the material heterogeneity in the current example. Nevertheless, as depicted in Fig. 14, an initial residual jump of AFETI-1 and FETI-DP methods is still significant for the heterogeneous material, while the present methods show a significantly smaller initial residual. In addition, despite the larger number of iterations, the present methods eventually converge to the accurate solution, while the AFETI-1 method exhibits a limited convergence threshold. Considering that the present and AFETI-1 methods originate from the same variational framework utilizing the localized Lagrange multiplier [30,31], such accuracy difference will highlight the benefits of the DO (displacement-only) approach. 7. Discussions Three iterative algorithms have been proposed for the solution of static structural problems employing the DO partitioned formulation [29], labeled as AFETI-DO, AFETI-DOc, and AFETI-DOci methods. These methods are primal alternatives to the dual Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 22 S.-H. Kang et al. Fig. 14. Fourth-order elastostatic example: Convergence history of global force residual for the case of 𝐻∕ℎ= 16 and 𝑁𝑠= 42. AFETI (AFETI-1) method based on the localized Lagrange multipliers, and employ the partitioned displacement (𝐝) as a fundamental unknown without introducing any interfacial variables. There is no direct assembly of interfacial constraints that are enforced by projector preconditioning or deflation with an additional stiffness regularization. While offering simpler implementation, the proposed methods yield significantly improved solution accuracy, especially for heterogeneous (stiff) problems. We first examined the AFETI-DO method as the basic simple PCG algorithm, parallel implementation of which only requires partitioned stiffness, its block-diagonal preconditioner, and negligible computations involving the coupling projector  𝑑for each processor. Several specific preconditioners are tested, including lumped inverse, pseudo-inverse, and regularized pseudo-inverse. Although its iteration count may increase with the number of subdomains, AFETI-DO yields converged solutions for the examples, while both AFETI-1 and FETI-DP methods fail to converge in some heterogeneous examples. Next, the AFETI-DOc method is shown as a two-level primal method to possess the accelerated convergence speed by introducing a coarse space based on the d’Alembert–Lagrange principal equation. Despite being somewhat more complicated than the AFETIDO method, it preserves the convergence rate for an increasing number of subdomains and allows the enrichment with numerically computed eigenvectors, while maintaining the accuracy. For the homogeneous material, a comparable convergence rate to the traditional iterative methods (there are many) is observed, except for the fourth-order elasticity example. For the highly heterogeneous material, the AFETI-DOc method exhibits superior accuracy compared to some of the existing methods. As an alternative approach to the AFETI-DOc method, the AFETI-DOci method is a primal method with additional projector preconditioner based on the partial decomposition of displacement into deformational and rigid-body motions, introducing the relatively simpler coarse Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 23 S.-H. Kang et al. space projector. For the numerically scalable examples, AFETI-DOci method shows almost equivalent convergence performance to that of the AFETI-DOc method. However, the present methods exhibit limitation for the irregular coefficient discontinuity and the fourth-order elasticity problems. The remedy for such weakness is currently being investigated. We have evaluated the performance of the AFETI-DO (and its extensions, viz., AFETI-DOc and AFETI-DOci methods) method with a focus on their accurate convergence properties as we believe efficiency matters only if an algorithm yields accurate solutions. Based on the performance of AFETI-DO (and its initial variations) for the limited yet representative heterogeneous problems, we submit that AFETI-DO is far simpler to implement than nearly all existing domain decomposition methods, and yields accurate converged solutions for all the test problems, while its iterations numbers may increase for some problems. In no case, the AFETI-DO algorithm failed to converge during iterations for all of the test problems. The accuracy advantage, implementation simplicity, and computational efficiency of the proposed AFETI-DO methods, we believe, originate from the properties of the unassembled block diagonal stiffness matrix (e.g., each of the partitioned stiffness matrices can be regularized individually without affecting the adjacent partitions). This and additional properties of the DO partitioned equations will be further investigated as applied not only to static problems but other problems such as flexibility-based solution algorithms [33], model-order reduction [43], and localized damage detection techniques [44]. Encouraged by the numerical evidence of its iteration robustness, several efficiency improvements of the AFETI-DO algorithm are presently being evaluated and we hope to report our findings in the near future. CRediT authorship contribution statement Seung-Hoon Kang: Writing – original draft, Software, Methodology, Formal analysis, Conceptualization. K.C. Park: Writing – review & editing, Methodology, Conceptualization. José A. González: Writing – review & editing, Validation. SangJoon Shin: Writing – review & editing, Supervision. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability No data was used for the research described in the article. Acknowledgments This work was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF), funded by the Ministry of Science, ICT and Future Planning, Republic of Korea, South Korea (2021R1A2C1007352). This work was also supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT), South Korea (No. 2021R1A5A1031868). Appendix A. Characteristics of the projector ⊤ 𝒅 This section will explain the range space of the projector ⊤ 𝑑and the enforcement of interfacial compatibility condition by it. A.1. Range space of the projector ⊤ 𝑑 After multiplying ⊤ 𝑑by 𝐈−𝐁𝐁⊤, only the column vectors corresponding to the internal DOFs will be preserved, and the rest will become zero. Such multiplied matrix is denoted as ⊤ 𝑑,1, and will be derived as follows: ⊤ 𝑑,1=⊤ 𝑑[𝐈−𝐁𝐁⊤]=𝐈−𝐁𝐁⊤.(A.1) In a similar way, the matrix ⊤ 𝑑,2, in which the column vectors of interfacial DOFs will remain, will be derived as in Eq. (A.2). ⊤ 𝑑,2=⊤ 𝑑𝐁𝐁⊤=𝐁𝐁⊤−𝐌−1𝐁𝐌𝑏⊤ 𝑓𝐁⊤(A.2) Based on these two matrices, the range space of ⊤ 𝑑will be as below: (⊤ 𝑑)=(⊤ 𝑑,1)∪(⊤ 𝑑,2).(A.3) The matrix 𝑑,1is the diagonal matrix in which the entries corresponding to internal and interfacial DOFs are 1 and 0. Therefore, its range space and rank will be as follows: (⊤ 𝑑,1)=(𝐁⊤),rank (⊤ 𝑑,1)= rank ((𝐁⊤))=𝑁𝑖.(A.4) Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 24 S.-H. Kang et al. Herein, 𝑁𝑖is the number of the internal DOFs in the partitioned system. The matrix ⊤ 𝑑,2will be further divided into ⊤ 𝑑,12 and ⊤ 𝑑,22, as below: (⊤ 𝑑,2)=(⊤ 𝑑,12)∪(⊤ 𝑑,22),(A.5) ⊤ 𝑑,12 =[𝐈−𝐁𝐁⊤]⊤ 𝑑,2,(A.6) ⊤ 𝑑,22 =𝐁𝐁⊤⊤ 𝑑,2=𝐁[𝐈−𝑓]𝐁⊤=𝐁𝐋𝑓𝐌−1 𝑓𝐋⊤ 𝑓𝐌𝑏𝐁⊤.(A.7) The range space of ⊤ 𝑑,12 will be restricted as follows: (⊤ 𝑑,12)⊆(𝐁⊤).(A.8) From the interface–frame connectivity and the full-ranked property of 𝐌𝑓and 𝐌𝑏, the range space and rank of ⊤ 𝑑,22 will be as below: (⊤ 𝑑,22)=(𝐁𝐋𝑓),rank (⊤ 𝑑,22)= rank (𝐁𝐋𝑓)=𝑁𝑓.(A.9) Herein, 𝑁𝑓is the number of frame DOFs in the partitioned system. Considering Eqs. (A.3)–(A.5),(A.8), and (A.9), the range space and rank of ⊤ 𝑑will finally become as follows: (⊤ 𝑑)=(𝐁⊤)∪(𝐁𝐋𝑓)=(𝐋𝑔),rank (⊤ 𝑑)=𝑁𝑖+𝑁𝑓=𝑁𝑔.(A.10) Herein, 𝑁𝑔is the number of DOFs of the assembled system. Because Eq. (A.10) is satisfied for ⊤ 𝑑with an arbitrary SPD matrix 𝐌, the range space of  ⊤ 𝑑in Section 3.2 will also become ( ⊤ 𝑑)=(𝐋𝑔). A.2. Interface compatibility condition by the projector ⊤ 𝑑 The interfacial compatibility equation will be recalled from Eq. (7), as follows: 𝐁⊤𝐝−𝐋𝑓𝐝𝑓=𝟎.(A.11) By multiplying the idempotent projector ⊤ 𝑓in Eq. (9), the frame component 𝐋𝑓𝐝𝑓will be filtered out, and another form of displacement-only interfacial relationship will be obtained as in Eq. (A.12). ⊤ 𝑓𝐁⊤𝐝=𝟎(A.12) Eq. (A.12) demonstrates that the partitioned interfacial displacement 𝐁⊤𝐝will satisfy the interfacial compatibility condition if it is included in 𝐁⊤𝐝∈(⊤ 𝑓). Based on such observation, the projector ⊤ 𝑑in Eq. (9) will be divided as below: ⊤ 𝑑=𝐈−𝐌−1𝐁𝐌𝑏⊤ 𝑓𝐁⊤ =𝐈−𝐌−1𝐁𝐌𝑏𝐁⊤ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟ ⊤ 𝑑,𝐼 +𝐌−1𝐁𝐌𝑏[𝐈−⊤ 𝑓]𝐁⊤ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ⊤ 𝑑,𝐵 .(A.13) Then, the interfacial components will be extracted by multiplying 𝐁⊤, as follows: 𝐁⊤⊤ 𝑑,𝐼 =𝟎,(A.14) 𝐁⊤⊤ 𝑑,𝐵 =[𝐈−⊤ 𝑓]𝐁⊤.(A.15) Therefore, multiplying ⊤ 𝑑will enforce the arbitrary vector to satisfy the interfacial compatibility condition, which will be also described as below: ⊤ 𝑓𝐁⊤⊤ 𝑑=𝑓𝐁⊤⊤ 𝑑,𝐼 +𝑓𝐁⊤⊤ 𝑑,𝐵 =𝟎+⊤ 𝑓[𝐈−⊤ 𝑓]𝐁⊤=𝟎.(A.16) As shown, the projector ⊤ 𝑑,𝐼 eliminates the interfacial components. Meanwhile, the projector ⊤ 𝑑,𝐵 complements the correct interfacial components by filtering out the components in (𝑓)using the projector 𝐈−⊤ 𝑓. Appendix B. AFETI-DOc and AFETI-DOci algorithms This section describes the algorithmic implementation of AFETI-DOc and AFETI-DOci methods contained in Algorithms 2and 3, respectively. Computer Methods in Applied Mechanics and Engineering 431 (2024) 117273 25 S.-H. Kang et al. Algorithm 2 AFETI-DOc algorithm. 1: 𝐝0= 𝐑𝑑[ 𝐑⊤ 𝑑𝐊 𝐑𝑑]−1  𝐑⊤ 𝑑𝐟,𝐟=𝐋+ 𝑔𝐟𝑔 2: 𝐫0= 𝑑{𝐟−𝐊𝐝0} 3: 𝐳0=𝐀+𝐫0 4: 𝐬0= ⊤ 𝑑 ⊤ 𝑐𝐳0 5: 𝐩0=𝐬0 6: for 𝑘= 0,⋯, 𝑚 do 7: 𝐪𝑘=𝐊𝐩𝑘⊳Fully parallelizable stage 8: 𝛼𝑘=𝐫⊤ 𝑘𝐬𝑘 𝐩⊤ 𝑘𝐪𝑘 9: 𝐝𝑘+1 =𝐝𝑘+𝛼𝑘𝐩𝑘 10: 𝐫𝑘+1 =𝐫𝑘−𝛼𝑘 𝑑𝐪𝑘⊳Inter-processor stage 11: if the solution converges then 12: exit the loop 13: end if 14: 𝐳𝑘+1 =𝐀+𝐫𝑘+1 ⊳Fully parallelizable stage 15: 𝐬𝑘+1 = ⊤ 𝑑 ⊤ 𝑐𝐳𝑘+1 ⊳Inter-processor stage with the coarse space projection 16: 𝛽𝑘=𝐫⊤ 𝑘+1𝐬𝑘+1 𝐫⊤ 𝑘𝐬𝑘 17: 𝐩𝑘+1 =𝐬𝑘+1 +𝛽𝑘𝐩𝑘 18: end for Algorithm 3 AFETI-DOci algorithm. 1: 𝐝0∈(𝐋𝑔),𝐟=𝐋+ 𝑔𝐟𝑔 2:  𝐫0= 𝑐𝑖  𝑑{𝐟−𝐊𝐝0} 3: 𝐳0=𝐀+ 𝐫0 4: 𝐬0= ⊤ 𝑑 ⊤ 𝑐𝑖𝐳0 5: 𝐩0=𝐬0 6: for 𝑘= 0,⋯, 𝑚 do 7: 𝐪𝑘=𝐊𝐩𝑘⊳Fully parallelizable stage 8: 𝛼𝑘= 𝐫⊤ 𝑘𝐬𝑘 𝐩⊤ 𝑘𝐪𝑘 9: 𝐝𝑘+1 =𝐝𝑘+𝛼𝑘𝐩𝑘 10:  𝐫𝑘+1 = 𝐫𝑘−𝛼𝑘 𝑐𝑖  𝑑𝐪𝑘⊳Inter-processor stage with the coarse space projection 11: if the solution converges then 12: exit the loop 13: end if 14: 𝐳𝑘+1 =𝐀+ 𝐫𝑘+1 ⊳Fully parallelizable stage 15: 𝐬𝑘+1 = ⊤ 𝑑 ⊤ 𝑐𝑖𝐳𝑘+1 ⊳Inter-processor stage with the coarse space projection 16: 𝛽𝑘= 𝐫⊤ 𝑘+1𝐬𝑘+1  𝐫⊤ 𝑘𝐬𝑘 17: 𝐩𝑘+1 =𝐬𝑘+1 +𝛽𝑘𝐩𝑘 18: end for References [1] J.W. Demmel, M.T. Heath, H.A. van der Vorst, Parallel numerical linear algebra, Acta Numer. 2 (1993) 111–197, http://dx.doi.org/10.1017/ S096249290000235X. [2] P.R. Amestoy, I.S. Duff, J.Y. L’Excellent, Multifrontal parallel distributed symmetric and unsymmetric solvers, Comput. Methods Appl. Mech. Engrg. 184 (2–4) (2000) 501–520, http://dx.doi.org/10.1016/S0045-7825(99)00242-X. [3] P. D’Ambra, F. Durastante, S. Filippone, AMG preconditioners for linear solvers towards extreme scale, SIAM J. Sci. Comput. 43 (5) (2021) S679–S703, http://dx.doi.org/10.1137/20M134914X. [4] S. Nikolopoulos, I. Kalogeris, G. Stavroulakis, V. Papadopoulos, AI-enhanced iterative solvers for accelerating the solution of large-scale parametrized systems, Internat. J. Numer. Methods Engrg. 125 (2) (2024) e7372, http://dx.doi.org/10.1002/nme.7372. [5] A. Toselli, O. Widlund, Domain Decomposition Methods - Algorithms and Theory, Springer-Verlag, Berlin, 2005. [6] C. Farhat, F.X. Roux, A method of finite element tearing and interconnecting and its parallel solution algorithm, Internat. J. Numer. Methods Engrg. 32 (6) (1991) 1205–1227, http://dx.doi.org/10.1002/nme.1620320604. [7] C. Farhat, J. Mandel, F.X. Roux, Optimal convergence properties of the FETI domain decomposition method, Comput. Methods Appl. Mech. Engrg. 115 (3–4) (1994) 365–385, http://dx.doi.org/10.1016/0045-7825(94)90068-X. [8] J. Mandel, Balancing domain decomposition, Commun. Numer. Methods Eng. 9 (3) (1993) 233–241, http://dx.doi.org/10.1002/cnm.1640090307.