The domain interface method in non-conforming domain decomposition multifield problems

The Domain Interface Method (DIM) is extended in this contribution for the case of mixed fields as encountered in multiphysics problems. The essence of the non-conforming domain decomposition technique consists in a discretization of a fictitious zero-thickness interface as in the original methodology and continuity of the solution fields across the domains is satisfied by incorporating the corresponding Lagrange Multipliers. The multifield DIM inherits the advantages of its irreducible version in the sense that the connections between non-matching meshes, with possible geometrically non-conforming interfaces, is accounted by the automatic Delaunay interface discretization without considering master and slave surfaces or intermediate surface projections as done in many established techniques, e.g. mortar methods. The multifield enhancement identifies the Lagrange multiplier field and incorporates its contribution in the weak variational form accounting for the corresponding consistent stabilization term based on a Nitsche method. This type of constraint enforcement circumvents the appearance of instabilities when the Ladyzhenskaya–Babuška–Brezzi (LBB) condition is not fulfilled by the chosen discretization. The domain decomposition framework is assessed in a large deformation setting for mixed displacement/pressure formulations and coupled thermomechanical problems. The continuity of the mixed field is studied in well selected benchmark problems for both mixed formulations and the objectivity of the response is compared to reference monolithic solutions. Results suggest that the presented strategy shows sufficient potential to be a valuable tool in situations where the evolving physics at particular domains require the use of different spatial discretizations or field interpolations.


Introduction
The growing demand of industrial complex simulations, together with the continuous improvement of computer technology, have stimulated the popularity of algorithms capable of processing large systems of equations. This is the case of algorithms based on domain decomposition techniques which basically stem from the need of decomposing a large discretization into a number of smaller discretizations, e.g. using finite elements (FEs), or from the assembly of independently generated meshes which form part of a whole system. The latter scenario can be found for instance during simulations of an aircraft where both wings and fuselage discretizations are designed independently. In these situations, the type of algorithms capable of gluing non-matching meshes are crucial and represent the main focus in the present contribution.
Examples of non-conforming interfaces arising from the assembly of non-matching grids can be found in [24,25] where particular structural components are reused in evolv-  (1) Ω (2) Γ I Ω (1) Ω (2) Γ I Non-matching limits Ω (1) Γ I Ω (1) Ω (2) Γ I refinement h/p-General non-conforming Ω (2) Geometrically compatible Geometrically incompatible Piece-wise linear curves ing designs such as wings or rotor blades among diverse aircraft fuselages. The most complex and time-evolving situations for non-conforming interfaces can be definitely found in the field of contact mechanics [16,28,31,41] which has boosted the appearance of a considerable number of strategies to handle non-matching grids. It is also necessary to resort to these type of algorithms when the nature of the physics involved in each of the studied domains is remarkably different, e.g. fluid structure interaction [10] or because the spatial resolution for the adjacent discretizations is significantly different, e.g. multiscale and multiresolution techniques [21,22]. The most general scenario at a non-conforming interface I arises when the geometrical compatibility between the non-matching meshes is not satisfied, i.e. when the boundaries of the domains at the common discretized interface are not identical in the undeformed configuration (cf. Fig. 1).
Many Domain Decomposition techniques assume that the interface is geometrically compatible [33] but more general situations may arise from the assembly of independent discretizations and curved interfaces with piece-wise linear distinct discretizations. In such scenarios, powerful techniques are required such as mortar techniques [5,[30][31][32]40] or the recently introduced DIM [8].
In mortar methods, the coupling constraints are enforced in a weak sense, i.e. in an integral or average sense along portions of the interface, and are often referred to as segmentto-segment techniques. In general, a displacement field at the interface, u I , is taken as reference in order to generate the constraints. A frequently adopted choice is to consider the reference field based on one of the domain discretizations at the common interface.
One significant drawback of such classical mortar techniques arises when the coarsest discretization is selected for the reference solution field. In such scenarios, the performance of mortar methods might be less optimal as argued in [8,18]. More advanced mortar techniques [30][31][32]40] consider a more convenient Lagrange multiplier space or even a third auxiliary surface with an optimal node distribution in order to minimize the gap between domains (cf. Fig. 2). This may circumvent these flaws although the associated computational cost due to the extra unknowns may increase as well as the complexity of the resulting algorithm.
In contrast with established mortar methods, the DIM [8] is based on a continuous Delaunay discretization of the interface obtained after a fictitious contraction of the adjacent domains (cf. Figs. 2,3). This technique has proven to achieve equivalent optimal performances without the drawback of selecting a projection surface since this is automatically accounted for by the Delaunay interface discretization. Moreover, the adjacent interface patch discretizations has been utilized to remove the singularities that arise when the decomposed domain is floating, i.e. exhibits rigid body modes (RBMs). This is performed in a non-intrusive fashion compared to other established dual formulations.
Lagrange multipliers [1] are typically employed to enforce the necessary constraints at general interfaces, e.g. possibly non-conforming. To this end, the term concerning the mechanical work at the non-conforming interface is added to the variational statement. Special care needs to be taken when selecting the shape functions associated to the mesh discretizations for both adjacent domains and the distribution functions of the Lagrange multipliers. In other words, the selection of such functions can not be arbitrary since they need to fulfill the Ladyzhenskaya-Babuška-Brezzi (LBB) condition (also known as the inf-sup condition) [2] in order to guarantee that both discretizations converge to the right solution upon mesh refinement. In the current approach a Lagrange multiplier method is adopted to enforce the interface constraints but a Nitsche method [27] is considered to stabilize the system. Such method essentially modifies the weak form by adding a term with a positive parameter that depends on the particular mesh size but does not lead to an ill-conditioned system even for small values of the parameter. A Nitsche method to account for the interface constraints derived from domain decomposition methods has been introduced in [4] and the extra term added at the variational principle basically couples the multipliers with the stress fields at the interface [17]. A main advantage is that the discretization of the hybrid solution field contain- ing the Lagrange multipliers can be performed without any constraints since the stability is accounted for by the extra Nitsche-type term. In the present contribution, this feature is specially valuable since the multifield nature of the problem renders a particularly hybrid solution field. Some generalities of the DIM are revisited in Sect. 2 and the multifield extensions to the method in order to tackle mixed incompressible u/ p and coupled thermomechanical u/ p/θ problems are detailed in Sects. 3 and 4. The reader that is not familiar with the constitutive model details for the mixed formulations is referred to Appendices 1 and 1 and the references therein for a complete overview. The domain decomposition framework is validated considering benchmark examples for both mixed type formulations in Sect. 6.

General description
The geometrical aspects of the DIM introduced in [8] are summarized in this section for the case of an irreducible formulation. The concept of geometrical gap and the Lagrange multipliers identification are revisited and the basic notation is introduced. Finite strain theory is considered in our developments, although an infinitesimal deformation approach can be recovered by considering small displacements (i.e. compared to the domain dimensions) and infinitesimal displacement gradients. Compact notation is utilized for tensor quantities in this contribution.

Geometrical description of the DIM
Consider the union of two non-overlapping domains (1) and (2) as depicted in Fig. 3  An independent FE discretization of the two bodies, leads to a number of N (s) λ vertices at the domain boundary ∂ (s) around (s) I . These vertices are involved in the interface discretization that characterizes the DIM (cf. step (i) in Fig. 3). The interface generation is summarized in the chart sequence in Fig. 3 and basically starts with a fictitious contraction of the vertices V i in the direction −ν ν ν (s) by a magnitude k as indicated in step (ii). As discussed in [8], the magnitude of k does not interfere with the results although k ≈ h e is commonly adopted in our analyses where h e denotes an average equivalent FE size.
A Delaunay triangulation defining the interface domain is obtained based on the resulting fictitious coordinates x i of the new vertices. The interface patches are denoted by D (q) and the interface surface Note that for D the case of a geometrically compatible interface, (s) I = (s) D since the trace of the interface coincides with the interface discretizations of the adjacent domains. Triangular linear elements D (q) are considered in our analyses for the discretization of the interface. In this situation, all integrals over a geometrically incompatible interface converge to a bounded value when the distance h → 0 Note that the error due to the piece-wise linear approximation of the interface is O(h). For the case of geometrically compatible interfaces O(h) → 0 and The main strength of the method is that both adjacent domain interfaces (s) D are automatically taken into account in the construction of the interface links since the integration is performed over D (q) when h → 0. This represents an obvious advantage from the user point of view considering that no decision has to be made on which interface discretization is selected for the projection of the interface quantities, e.g. in the mortar method [5].
The corresponding deformation map for each domain (s) is denoted by φ (s) where material points X (s) ∈ t . Such current configuration x (s) can be also obtained in terms of the total displacement field U (s) t (X (s) ) as A pseudo-time domain t ∈ [0, T ] is considered in our computations with subdivisions in discrete intervals [t n , t n+1 ] of incremental time length t = t n+1 −t n . Configurations at the previous t n and current t n+1 times are denoted as As shown in [8], an expression of the incremental motion φ (s) can be found by substitution of the current and previous configurations leading to Considering the incremental motion in (5), the incremental field The incremental motion of the interface domain can be expressed as where D (q) n denotes the interface domain D at time t n . The current and previous domain interfaces are related by D n+1 = φ D (D n ) and γ D = φ D ( D ) sets the link between current and previous interface surfaces, respectively (cf. Fig. 4).
A linear interpolation of the displacement incrementsû D i corresponding to the interface element vertices is employed to compute the incremental displacement field at the interface domain u D as n , (8) where N u,D i constitute the linear shape functions for threenode triangular finite elements employed to interpolate the interface displacements. The incremental gradient deformation tensor f D can be expressed as 1 being the second order unity tensor and∇ ∇ ∇ denoting the material gradient with respect to the reference previous configuration n. Due to the linear interpolation of the incremental displacements in (8), the incremental gradient deformation tensor f D is constant at every patch D (q) . As introduced in [8] and illustrated in Fig. 4, the normal to the base-line of D (q) in the sense of the normal to the adjacent domain ν ν ν (s) is denoted by N (q) and T (q) =ê×N (q) , whereê denotes the out-of-plane unit vector. With these definitions in hand and considering the incremental motion φ D , the current tangential and normal unit vectors read It should be stressed that these vectors are constant within every patch (cf. (10)) but discontinuous across the interface patches since their definition depends only on the local baseline for every interface patch. The initial normal gap g 0 N depicted in Fig. 4 is obtained at the previous configuration n for a given point x n and its normal projection to the base-linex n as Consequently, the final gap vector where x n+1 andx n+1 stand for the convected points x n and x n , respectively. Expressing g(x n ) as a sum of its normal and tangential projections onto the current base-side, equation (12) can be rewritten as The normal gap g N (x n ) is interpreted as a projection in the direction of n (q) and, consequently, denotes penetration when g N (x n ) < 0. In the same spirit, the tangential gap g T (x n ) stands for the slid distance in the sense of t (q) . The gap is conveniently expressed in terms of the displacement field as detailed in [8]. It essentially considers a Taylor series expansion of φ D (x n ) aroundx n up to second order terms in order to obtain the dependency on the displacements u (q) . A dimensionless measure of the gap, i.e. the gap intensityḡ(x n ) is defined as Finally, one obtains the expression of the gap ready to be utilized in the variational statements as The interface traction vector t I defined at the surface (s) D is expressed in terms of the first Piola-Kirchoff stress tensor P (s) at time n + 1 with respect to the configuration at time n and the normal vector The normal and tangential components of t I are obtained w.r.t. the current normal and tangential vectors as and are depicted in Fig. 5 inside the surface γ (s) D . The normal and tangential components of the interface traction vector implicitly define the Lagrange multipliers that enforce displacement compatibility between domains at each interface patch D (q) as In this scenario, the Lagrange multipliers λ N and λ T are identified as the normal and tangential stresses that connect the adjacent domains (cf. Fig. 5).

Remark 2.1
The identification of the Lagrange multipliers λ N and λ T stem from the variational statement corresponding to the weak form of the governing equations where δu andV V V u denote the displacement variations and their corresponding function space (cf. Appendices 1 and 1 for a more detailed formulation). The energy functionals δ u , δ u int,ext and δ u I denote the total mechanical work, the mechanical work due to external and internal forces and the mechanical work performed at the interface I , respectively, and can be written as: wheret denote the imposed tractions. The term concerning the interface work δ u I can be further elaborated as 1 : where integration by parts is employed considering the expression of the gap as a function of the displacements in (15) and the constant character of the Lagrange multipliers λ λ λ u . In the limit where the h (q) tends to zero, the expression in (22) can be written as 1 In the sequel, all summations of integrals over the interface patches n , i.e. with q from 1 to the total number of interface patches N q , are denoted as integrals over D n skipping the summation over q in order to simplify the notation.
The total work in (19) can be rewritten considering that ∇ ∇ ∇ · P (s) · δu (s) = ∇ ∇ ∇ · P (s) · δu (s) + P (s) :∇ ∇ ∇ δu (s) and employing the divergence theorem taking into account all boundary contributions, i.e. = σ ∪ u ∪ I , as Note that the boundary integral over u vanishes according to the definition of the displacement variations δu (cf. (168) in Appendix 1). From the variational principle in (24) one can identify the corresponding Euler-Lagrange equations and natural boundary conditions as which concludes the proof that identifies the Lagrange multipliers λ λ λ u with the tractions t I that connect the adjacent domains.
The strong and weak forms of the irreducible problem employing the DIM are introduced in detail in [8] and skipped in the present approach in order to focus only on the additional formulation content. In this spirit, the following sections are devoted to present a complete variational formulation of the DIM for the mixed u/ p and thermo-mechanical mixed u/ p/θ formulations.

DIM for mixed u/ p formulations
The DIM method introduced in the previous sections, and described in detail in [8], is modified in this Section in order to tackle incompressible problems using mixed u/ p formulations (cf. Appendix 1). The kernel of the extension resides in the identification of the Lagrange multiplier λ p for the correct pressure transference between domains which forces the pressure gap g p to nullify. Figure 6 illustrates the Lagrange multiplier for both displacement λ λ λ u (u) and pressure λ p ( p) fields inside a triangular interface patch D (q) with nodal quantities u 1 to u 3 and p 1 to p 3 . The pressure p (q) at the interface patch D (q) is calculated by linearly interpolating the pressure unknowns p D i as n . (28) The geometrical pressure gap g p is defined for each discretized interface patch D (q) as and is understood as the pressure gradient projected onto the normal to the base side of the triangular patch N (q) . The effective pressure gapḡ p is found normalizing the pressure gap g p as

Strong and weak forms of the mixed u/ p DIM formulation
The strong form of the u/ p DIM problem is obtained considering the equations shown in Appendix 1 at each domain (s) and accounting for the interface connections between domains as: FIND: FULFILLING: Constitutive pressure equation: Dirichlet's boundary conditions: Neumann's boundary conditions: Lagrange multiplier identification: Compatibility constraints: ∇ ∇ ∇ being the material Nabla operator containing derivatives with respect to the previous configuration, P the first Piola-Kirchhoff stress tensor andt andû the prescribed tractions and displacements at the corresponding boundaries σ and u , respectively. The bulk modulus is denoted by κ and J = ||f|| is the determinant of the deformation gradient tensor. A proper justification of the pressure Lagrange multiplier identification in Eq. (36) is provided in Remark 3.2 after the introduction of the weak form of the problem.
As mentioned in Appendix 1, Eq. (32) is solved considering the splitting of the Cauchy stress tensor into the spheric and deviatoric counterparts.
The strong form of the problem stated in (32) to (37) is expressed through the equivalent variational principles. To this end, the solution V V V and weightingV V V spaces for the displacements and pressure fields together with the Lagrange Fig. 6 Lagrange multiplier identification at the discretized interface for the mixed u/ p DIM formulation multiplier space L L L for the solution and weighting functions read: H 1 ( ) and L 2 (D) being the Sovolev space of functions with square integrable derivatives and the Lebesgue space of square integrable functions, respectively. The variational statement for the mixed u/ p formulation reads: FIND: FULFILLING: AND The part of the mechanical work corresponding to internal and external forces reads and The work performed at the interface by the virtual gap reads where δḡ denotes the gap intensity variations. The interface normal and tangential work contributions can also be written in terms of the displacement variations following the expressions of the gap intensity variations developed in [16,28]. The work performed by the pressure Lagrange multipliers at the interface can be expressed as: The variational statements for the displacement constraint equations are introduced in detail in [8] and read: and force the interface work to nullify in an average sense along the domain interface D n . Accordingly, the new variational statement for the interface pressure constraint is expressed as: Remark 3.1 In order to stabilize the pressure field for the case of compressible and nearly compressible materials the projection term in (183) is added to the weak form of the constitutive equation concerning the incompressibility condition in (46). Such projection term arises from the Polynomial Pressure Projection (PPP) method briefly introduced in Appendix 1.

Remark 3.2
The identification of the pressure Lagrange multipliers λ p stem from the variational statement corresponding to the weak form of the governing equations in (46). The expression in (54) concerning the interface work δ p I developed by the pressure Lagrange multipliers can be further elaborated as: where integration by parts is employed considering the expression of the pressure gap in (29) and the constant character of the pressure Lagrange multipliers λ p . In the limit where h (q) tends to zero, the expression in (58) can be written as The total virtual work performed by the pressure variations reads and, therefore, the associated Euler-Lagrange equations with the corresponding natural boundary conditions can be identified as: which concludes the proof that identifies the pressure Lagrange multipliers λ p at the interface.

Discretization using FE and lambda-solvability of the mixed u/ p system
A Galerkin-based spatial discretization is considered in which the solution field components and their variations at each domain (s) are interpolated as indicated in (174) to (177) (cf. Appendix 1). The Lagrange multipliers and their variations connecting displacement and pressure fields are interpolated as where the subscript b stands for the the discrete nodes corresponding to the Lagrange multipliers interpolation using the domain-wise constant shape functions u and p which read: The corresponding FE approximations of the mixed variational statements in (45) to (49) are denoted as follows: Considering that (68) to (72) hold for any virtual quantity δû, δp, δ N , δ T and δ p , the resulting discrete problem consists in finding the nodal quantitiesû,p, u and p such that the following residual is nullified: The residuals concerning the mechanical R u int,ext and incompressible part R p are detailed in (180) and (181) (cf. Appendix 1). The discrete contributions of the normal and tangential Lagrange multipliers R u I as well as the contri-butionsR λ u has been already introduced in [8] for the irreducible formulation. The contribution of the pressure Lagrange multipliers at the interface and the corresponding residual of the pressure Lagrange multiplier restriction The discrete equations (73) to (76) are linearized and solved incrementally via a Newton-Raphson procedure leading to the following system: As argued in [8], the component ∂R λ ∂ u differs from zero due to the introduction of a stabilization term added to the constraint variational statements (70) and (71). Following an analogous procedure, a new stabilization term needs to be added to the variational statement related to the pressure constraint δ p λ p in order to provide a dependency of the residual component R λ p with respect to the Lagrange multipliers p . The modified residual component and essentially consists in the addition of a weak format of the pressure Lagrange multiplier identification in (37). It should be noted that, since the penalized term is part of the Euler-Lagrange equations of the variational principle (37), it will tend to zero upon mesh refinement. For this reason the stabilization procedure described in (80) is understood as a consistent penalty method in which, unlike other nonconsistent penalty methods, the parameter τ p > 0 can be made significantly small without affecting the quality of the obtained results. Obviously, arbitrarily large values of τ p may lead to an ill-conditioning of the resulting system. However, as noted in the example reported in Sect. 6.1, results are remarkably insensitive to the value of τ p since a large difference in such values barely affects the resulting pressure field.
The new modified residual and the linearized system in (79) can be rewritten as Remark 3.3 Note that the stabilization parameter τ p has no dimensions associated with it as opposed to the stabilization parameter τ utilized for the component R λ and conveniently justified in [8].
The dual assembly in (82) is rewritten considering the discretized quantities for each domain (s) as (1) . . . where and (s) and ∀ = (s) .
(89) Fig. 7 Lagrange multiplier identification at the discretized interface for the coupled thermomechanical u/ p/θ DIM formulation Different parallel solution strategies for solving such a system are briefly discussed in Sect. 5 together with an iterative scheme for a general non-linear problem.

DIM for coupled thermomechanical u/ p/θ formulations
The coupled thermomechanical formulation utilized in this manuscript is explained in detail in Appendix 1 for the case of a monolithic approach. In this manner, the present section is exclusively focused on the necessary extensions to the DIM in order to tackle such mixed problems. The kernel of the coupled thermomechanical DIM consists in the identification of the appropriate Lagrange multipliers concerning the connection of the displacement, pressure and temperature fields across the interface (cf. Fig. 7). The main addition with respect to the mixed formulation detailed in Sect. 3 is the introduction of the temperature unknowns θ which are transferred through the corresponding Lagrange multipliers λ θ located at each interface patch D (q) . The temperature θ (q) at the interface patch D (q) is calculated by linearly interpolating the temperature unknowns θ D i as The heat flux defined at the interface patch D (q) is denoted as Q (q) and Q (q) e refers to the heat flux at the FE adjacent to the base-side of the interface patch. In this manner, the Lagrange multiplier λ θ is defined as the the projection of the adjacent flux Q (q) e to the normal N (q) of the base-line as: The thermal gap g θ is defined for each discretized interface patch D (q) as and is understood as the temperature gradient projected onto the normal to the base side of the triangular patch N (q) . In order to provide a continuous heat flux across the interfaces the thermal gap is forced to be null. The effective temperature gapḡ θ is found normalizing the temperature gap g θ as

Strong and weak forms of the coupled thermomechanical u/ p/θ DIM formulation
The strong form of the mixed u/ p/θ formulation introduced in detail in Appendix 1 is considered at each domain (s) accounting for the interface restrictions. In this view, the strong form of the DIM thermomechanical formulation can be stated as: FIND: FULFILLING: Constitutive pressure equation: Energy balance: ρ 0U (s) +∇ ∇ ∇ · Q (s) = D (s) Dirichlet's boundary conditions: Neumann's boundary conditions: Lagrange multiplier identification: Compatibility constraints: where the domain energy balance equation, the identification of the temperature Lagrange multiplier and the temperature gap restriction are the new ingredients compared to the strong form presented in Sect. 3. The strong form in (95) to (101) is expressed through the equivalent variational principles in which the weak form presented in Sect. 3 is extended accounting for equations concerning the temperature field.
To this end, the solution V V V and weightingV V V spaces for the displacements, pressure and temperature fields together with the Lagrange multiplier space L L L for the solution and weighting functions read: H 1 ( ) and L 2 (D) being the Sovolev space of functions with square integrable derivatives and the Lebesgue space of square integrable functions, respectively. The variational statement for the mixed u/ p/θ formulation reads: FIND: FULFILLING: AND The part of the thermal work corresponding to internal and external contributions and whereR (s) and D (s) stand for the heat sources and the outward normal heat flux term applied at the interface n The work performed at the interface by the Lagrange multipliers λ θ reads where δg θ denotes the thermal gap variations.

Remark 4.1
The identification of the temperature Lagrange multipliers λ θ stem from the variational statement corresponding to the weak form of the governing equations in (114). The expression in (122) concerning the interface work δ θ I developed by the temperature Lagrange multipliers can be further elaborated as: where integration by parts is employed considering the expression of the temperature gap in (92) and the constant character of the temperature Lagrange multipliers λ θ . In the limit where h (q) tends to zero, the expression in (123) can be written as The total virtual work performed by the temperature variations reads and, therefore, the associated Euler-Lagrange equations with the corresponding natural boundary conditions can be identified as: which concludes the proof that identifies the pressure Lagrange multipliers λ θ at the interface.
The variational statement for the thermal constraint equations reads: and basically enforces in a weak sense the restriction to the thermal gap g θ = 0 (cf. equation (118)). In other words, the variational statement in (128) forces overall continuity of the heat flux across the interface. As argued in [8] and briefly commented in Sect. 3.2, the independence of the restriction statement with respect to the corresponding Lagrange multiplier, e.g. independence of (128) with respect to λ θ , leads to a system prone to exhibit instabilities if the adopted solution field discretization does not satisfy the LBB condition [2]. To this end, an extra term is added to the thermal restriction variational statement in order to provide the dependence with λ θ as The extra term is based in the expression in (100) and enforces the temperature Lagrange multiplier λ θ to be equal to the heat flux vector Q of the adjacent element projected on to the normal N (q) of the base-side of the corresponding patch D (q) (cf. Fig. 7).
The stabilization parameter τ θ > 0 is defined as α > 0 being an non-dimensional parameter, L stands for the length of the base-side of D (q) and k denotes the thermal conductivity coefficient.

Discretization using FE and lambda-solvability of the thermal system
As seen in the monolithic formulation (cf. Appendix 1), the coupled thermomechanical problem is resolved in a staggered way and, for this reason, only the discretized temperature equations are outlined in this section. The final discretized system needs to be solved together with the mechanical system in (83) outlined in Sect. 3. A full iterative procedure for the solution of the coupled thermomechanical domain decomposition problem is outlined in Sect. 5. A Galerkin-based spatial discretization is considered in which the solution field components and their variations at each domain (s) are interpolated as indicated in (232) to (237) (cf. Appendix 1). The temperature Lagrange multipliers are interpolated as where the subscript b stands for the the discrete nodes corresponding to the Lagrange multipliers interpolation using the shape functions θ = p defined at each interface patch D (q) as described in (67). The corresponding FE approximations of the thermal variational statements in (119), (122) and (128) are denoted as follows: Considering that (133) and (134) hold for any virtual quantity δθ and δ θ , the resulting discrete problem consists in finding the nodal quantitiesθ and θ such that the following residual is nullified: where The discrete equations (135) and (136) are linearized and solved incrementally via a Newton-Raphson procedure leading to the following system: Assembling the contributions for each domain (s) the global system of equations reads: where and Different parallel solution strategies for solving such a system are briefly discussed in Sect. 5 together with an iterative scheme for the general non-linear thermomechanical problem.

Parallel system solution strategies and non-linear iterative scheme
Both mechanical (83) and thermal systems (142) are expected to be large and suitable to be solved using parallel solvers. As argued in [8], both sparse systems can be resolved in a robust manner using direct parallel solvers such as multi-frontal or block-LU methods [34]. Considering that multifield problems may lead to ill conditioned systems, i.e. due to the contrast between stiffness coefficients, this choice might be a suitable one since they are based on independent simultaneous factorizations of the domain matrices and specific, i.e. problem dependent, preconditioners are not required. However, a considerable amount of memory is required and scalability is only expected in a moderate number of processors. As pointed out in [12,23], block-LU solvers account for automatic load-balancing and multi-threading which is optimal when dealing with highly heterogeneous domains in terms of the number of DOF.
In order to obtain scalabilities in a massive number of processors it is convenient to recast the global systems (83) or (142) in terms of the interface flexibility problem 2 considering that the matrices K (s) dd can be inverted, i.e. they do not exhibit singularities arising from the existence of rigid body modes (RBMs). The flexibility of the interface F I can be understood as the condensation of the domain stiffness matrices at the interface whereas g I is, in turn, the condensed residual force vector at the interface I .
After solving the flexibility problem in (147) and obtaining the Lagrange multiplier field , the domain displacement increments d (s) can be independently calculated for each domain as A blend of direct solvers are employed to independently compute the factorizations of the domain stiffness and an iterative solver is utilized for the solution of the interface problem in (147) which is never assembled in practice and, thus, requiring a considerably low memory profile [13,14]. Since the resolution of the interface problem and the computation of the domain solution fields are inherently parallel tasks, the methodology scales well in massively parallel computers.
If one or more subdomains (s) exhibit rigid body modes (RBMs), the corresponding matrices K The new expressions basically extend the unknown field d accounting for the rigid body mode intensities α α α as explained in [8,20,22]. This methodology is general for the parallel processing of any dual system (83) or (142) but its implementation in commercial FE packages is regarded highly intrusive.
In this view, exploiting the novel features of the DIM method, i.e. a continuous fictitious discretization of the interface, a new methodology to handle rigid body modes has been proposed in [8] and adopted in this contribution. The method is described in detail in [8] and essentially adds a new term to the variational statements in (45) to (46) or (112) to (114) which penalizes a function of the type 1 2 ||ḡ|| 2 . It is important to mention that the new term is nullified in the solution since it is related with the Euler-Lagrange equations of the constraint variational principle, this giving rise to the term "consistent penalty". The new variational statements in (45) to (46) and (112) to (114) read: where (r ) D corresponds to the interface segments such that r ) , N r representing the number of patches utilized to handle the RBMs (cf. Fig. 8) and c denotes the penalty coefficient utilized to enforce the new condition. A more detailed formulation considering the weak and discretized contributions of the stabilizing interface patches to the global system can be found in [8] for the irreducible formulation.
Remark 5.1 Note that the term cθ in (137) allows for the computation of K −1 θθ since it behaves similarly than a mass matrix in a dynamic setting. However, for the choice of null specific heat c = 0 the RBMs need to be eliminated with the addition of stabilizing interface contributions as indicated in (153) leading to a modified matrix of the type c θ > 0 corresponding to the penalization term.
The linearized set of equations in (83) and (142) obtained with a Newton-like scheme is solved iteratively for each load/time step t (cf. Newton-Krylov-Schur methods [9,15]). In this view, two types of iterations can be identified. A first type refer to the solution of the non-linear problem with successive linear approximations and second type arise from the solution of the flexibility problem in (147) where usually Conjugate Gradient or GMRES iterates are considered. The Schur complements are utilized for the local solutions at each domain (s) as indicated in (150).
Assuming a fix domain decomposition and a given Delaunay interface discretization, the iterative scheme for the coupled thermomechanical non-linear DIM framework is summarized in Box 1.

Framework validation through representative simulations
Two representative tests are selected to assess the performance of the multifield domain decomposition framework. Attention is focused on the continuity of the hybrid solution field at the interface, the convergence rate upon mesh refinement and the influence of the stabilization parameters as well as the activation of the Lagrange multipliers for the different solution fields at the interface. Finite strain theory and plane strain conditions are considered in all two dimensional examples.

Mixed u/ p formulation for incompressible problems
The mixed u/ p formulation is tested in the context of the DIM by solving the Cook's membrane problem. This test constitutes a well known benchmark problem for assessing the performance of the mixed u/ p element in both compressible and incompressible conditions (cf. [26,37]). The setup is basically sketched in Fig. 9 where a vertical force F is applied at one end of the membrane while the opposite end is fixed in both directions. The membrane is considered elastic with Young's modulus E = 70 N/mm 2 and the Poisson's ratio ν is set to 1/3 and nearly 0.5 for the compressible and incompressible case, respectively.
In order to assess the convergence of the proposed methodology, five different discretization scenarios are considered. These concern a monolithic approach and the DIM on the body partitioned in two domains, i.e. (1) and (2) depicted in Fig. 9, with both conforming and non-conforming interfaces (cf. Table 1; Fig. 10).
Parallel solver: • Factorize local stiffness and compute the flexibility components F I,i and g I,i • Solve for i+1 = F I,i −1 g I,i using an iterative solver. The vertical displacement of point A (cf. Fig. 9) is monitored for both monolithic and domain decomposition analyses employing the five different FE discretizations. Figure 11 shows the convergence behaviour in terms of the vertical displacement for the compressible and incompressible cases upon different mesh refinements.  Table 1 Number of finite elements for the five discretization scenarios considered at the Cook's membrane test. Each discretization is particularized for the monolithic approach on the whole domain and for the domain interface method (DIM) with conforming and non-conforming interfaces at partitions (1) and (2)  It is observed that all analyses provide an analogous convergent behaviour and the differences between the monolithic approach and the domain decomposition strategy are hardly visible beyond the third discretization, i.e. around 200-250 FE. In both compressible and incompressible analyses the DIM with non-conforming interfaces shows a slightly slower convergence with respect to the domain decomposition with conforming interfaces.
The influence of the activation of the pressure Lagrange multiplier λ p is studied for the incompressible case employing both conforming and non-conforming interfaces (cf. Fig. 12).
It is observed that the activation of λ p hardly influences the displacement solution field and the monitored vertical displacement is almost identical in all cases. Consequently, the Lagrange multiplier λ p does not seem to play a relevant role in the correct transference of the displacements across the interfaces.
The influence of the activation of λ p on the correct transference of the pressure variable across the interface is studied by monitoring the norm of the pressure jump p at the interface 1,2 between domains (1) and (2) as where P 2 ( p 1 ) denotes the projection of the coarse domain interface pressure p 1 onto the fine interface discretization 2 . It can be observed that the error p L 2 reduces with decreasing element size h (cf. Fig. 13).
Moreover, the convergence rate is increased and the magnitude of the error is diminished when the pressure Lagrange multiplier λ p is activated. The benefits of the activation of λ p are illustrated in the pressure distribution plots across the interface 1,2 for different mesh discretizations (cf. Fig. 14).
The deactivation of λ p leads to a pressure jump across the interface 1,2 . This effect is obviously less visible for the fine discretizations but the benefits of activating the pressure Lagrange multiplier are clear among all tested FE meshes.
The influence of the stabilization parameter τ p is studied for the incompressible Cook membrane with a nonconforming interface. The pressure distribution along the segment AA (cf. Fig. 15) is plotted for τ p = 10 +4 and τ p = 10 −4 employing the FE discretization (3) (cf. Table1; Fig. 10). Results depicted in Fig. 15 (left) clearly show that large values of τ p lead to a pressure jump at the interface, located at x = 25 mm, whilst low values of the stabilization parameter lead to a continuous pressure distribution along the intersection between AA and 1,2 .
As shown in Fig. 15 (right), the error p L 2 diminishes upon mesh refinement for both values of the stabilization parameter. Moreover, for large values of τ p the convergence rate decreases and the overall magnitude of the error increases. It is noted that a difference in eight orders of magnitude for τ p leads to a pressure jump of around 0.02 Mpa.

Coupled thermomechanical formulation applied to "bulk-forming" problems
In this example, a non-isothermal elastoplastic test is considered, in which a cylindrical metallic specimen is compressed until it experiences a 60% reduction of its initial height as depicted in Fig. 16a. The test is known as 'upsetting of a billet' in the bulk metal forming field and is considered appropriate for the assessment of the correct transference of both temperature and pressure fields in the present framework. Adiabatic conditions are considered, i.e., no heat exchange is allowed between the body and its surrounding environment, and, therefore, the temperature increments are solely related to self heating induced by dissipation during large (3) DIM conforming DIM non-conforming FE discretization Monolithic approach plastic deformation. Perfect stick is considered between the rigid tool and the specimen leading to 'barreling' of the compressed billet.
As indicated in Fig. 16b, only a quarter of the specimen is considered due to the symmetry of the geometry and boundary conditions. In order to assess the correct transference of the pressure and temperature fields across the domain bound-aries, the solution obtained with a monolithic approach, i.e. considering only one domain, is compared to the one obtained with four and twelve partitions with non-conforming interfaces as depicted in Fig. 17. A FE discretization of around 700 three-node axisymmetric linear elements is considered in all analyses. The thermo-plastic model detailed in Appendix 1 is utilized with the material parameters summarized in Table 2.
The yield stress evolution is given by with The relations (156) to (158) describe a linear thermo-plastic softening widely used in metal plasticity for most steels in a temperature range from 300 to 400 K [36]. The temperature distribution for the three discretizations depicted in Fig. 17 are plotted in Fig. 18 at the final time step (t = 1 s). Note that the temperature fields of the domain decomposition analysis are remarkably comparable and equivalent to the one of the reference monolithic analysis. The close-ups at the non-conforming interface regions show Ω (1) Discretization (5) Discretization (5) Without activating λ p Activating λ p Discretization (1) Discretization (1) Discretization (3) Discretization (3) a continuous temperature field which evidences its correct transference throughout the domains. Figure 19 shows the pressure distribution for the three discretizations at the final time step (t = 1 s). As it was observed for the temperature field, the distribution of the pressure is continuous in all domain decomposition analyses across the non-conforming interfaces. It is observed that the pressure field at the top right corner of the specimen could be improved by using a finer discretization in both monolithic and domain decomposition approaches. However, the aim of the analyses is to proof that the results obtained with the domain decomposition framework are practically equivalent between them, and that no differences are observed when compared to the reference monolithic solution which is actually reflected in Fig. 15 Sensitivity of the stabilization parameter τ p . Left Pressure distribution along the AA segment across the non-conforming interface 1,2 considering the FE discretization (3). Right convergence of the pressure transference across the interface upon different FE sizes h at the coarsest domain (1) . Results are plotted for the incompressible case employing non-conforming interfaces  The mechanical response of the billet, considering the applied force and displacement at the top of the specimen, is shown in Fig. 20. The reference responses reported by Ponthot [38], Taylor and Becker [29] are considered together with the monolithic and domain decomposition tests with four and twelve decompositions, respectively. The mechanical response of the domain decomposition method does not depend of the number of decompositions and conformity of the interfaces and is in agreement with the monolithic response which, in turn, shows that the selected constitutive model matches the results provided in literature, particularly the ones by Ponthot [29].

Conclusions
The Domain Interface Method (DIM) presented in [8] is extended in this contribution for the case of mixed formulations encountered in multiphysics problems. Particularly, the domain decomposition framework is assessed in a large deformation setting for the particular case of mixed displacement/pressure formulations in incompressible materials and coupled thermomechanical problems with applications to bulk forming processes.
The extended DIM method for multiphysics problems essentially relies in an explicit discretization of the interface by means of a zero-thickness Delaunay triangulation as in the original methodology. This strategy allows setting a full tied contact in the most demanding scenarios, e.g. between geometrically incompatible interfaces. Continuity of the mixed solution field is satisfied across domains through the incorporation of Lagrange Multipliers. A Nitsche method is adopted in which a stabilization term is added at the constraint equations of the mixed field, i.e. one stabilization term per field. The incorporation of the stabilization term avoids the appearance of zero diagonal terms at the global system and instabilities are avoided if the LBB condition is not fulfilled by the chosen discretization. It is important to remark that, since the penalized term is part of the Euler-Lagrange equations of the variational principle, it will vanish upon mesh refinement (consistent penalty method).
The formulation presented in this manuscript is based on the identification of all Lagrange multipliers to connect the mixed field and their consideration in the weak variational form together with the corresponding consistent stabilization term. This technology is first assessed for a mixed displacement/pressure formulation to tackle incompressible materials in which a Polynomial Pressure Projection method is employed to account for the stability of the resulting systems within each domain. It is demonstrated that the pressure field is transferred with continuity across the non-conforming interfaces only if the Lagrange multiplier concerning the pressure field is taken into account, i.e. the continuity of the displacement field at the interface, enforced by the displacement Lagrange multipliers, does not guarantee the continuity of the pressure field at the interface. In fact, this effect vanishes upon mesh refinement since the stabilization term is no longer needed and vanishes too. It should be highlighted that the sensitivity of the results to the stabilization parameter is, however, very small. For the mixed displacement/pressure formulation, a variation of eight orders of magnitude in the pressure stabilization parameter would cause an insignificant jump in the pressure field.
The assessment of the DIM formulation for multiphysics problems is completed with a coupled thermomechanical problem with applications to bulk metal forming processes. In this case, a mixed displacement/pressure/temperature field is considered together with the corresponding Lagrange multipliers at the interface. Continuity of the three fields is obtained across all interfaces and the distribution of the mixed field is in very good agreement with the reference monolithic solutions. The mechanical response resulting from the upsetting of a billet studied in this example is remarkably close to those reported in literature which validates both the adopted constitutive model and domain decomposition technique. An alternative non-intrusive methodology to handle RBMs was introduced in the original contribution presenting the DIM method. This methodology essentially adds an extra stabilization term to the energy functional with contributions of all adjacent interface patches avoiding, in this manner, the calculation of a pseudo-inverse at floating domains. As a result, the band structure of the global system is preserved and, for this reason, possible parallel solution strategies employed in other dual domain decomposition methods can still be employed. In the analyses considering the mixed displacement/pressure problem no singular modes may arise from the pressure field. However, in the thermomechanical problem, singular temperature modes with constant temperature distribution may arise at floating domains. An extra stabilization term analogous to the term restricting the RBMs should be considered but in the present formulation the subdomain matrices with diagonal temperature coefficients can still be inverted due to the introduction of a transitory term which acts in a similar way than the mass matrix in a dynamics problem.
The DIM method for mixed fields inherits all advantages outlined for the irreducible formulation, i.e. the interface connections are set up as a result of the connectivity in the automatic Delaunay triangulation and do not rely in user criteria for the master and slave DOFs or intermediate interface projections as it is done in established techniques such as the mortar approach. For this reason, the methodology shows potential to be used in cases where completely different discretizations are needed among domains since different The weak form of the problem stated in (160) to (163) is expressed through the corresponding variational principle. The solution V V V and weightingV V V spaces for the displacements and pressure fields read: H 1 ( ) and L 2 (D) being the Sovolev space of functions with square integrable derivatives and the Lebesgue space of square integrable functions, respectively. The variational statement for the mixed u/ p formulation reads: FULFILLING: AND

Discretization using FE
A Galerkin-based spatial discretization is considered in which the solution field components and their variations are interpolated as: where the subscript a denotes the discrete nodes corresponding to the displacement interpolation and N mec refer to the shape functions utilized for the displacement and pressure fields.
The FE approximation of the variational statements in (172) and (173) considering the additional stabilization term in (183) can be written using the above interpolations (174) to (177) as follows: Considering that (178) and (179) hold for any virtual displacement δu and pressure δp, the residual vector of the variational principle accounting for the stabilization term is composed by The discrete equations (180) and (181) are linearized and solved incrementally via a Newton-Raphson procedure leading to the following system:

Stabilization of the pressure field
The mixed displacement/pressure (u/ p) formulation presents oscillations for incompressible and nearly incompressible materials. Such oscillations arise when the same interpolation order utilized in both the displacements and pressure fields does not satisfy the LBB condition [2]. In the present study the Polynomial pressure projection (PPP) [6] is employed since, in contrast with other methods, it does not require a stabilization parameter depending on the mesh size or the calculation of higher-order derivatives.
The following projection term is added to the weak form of the constitutive equation involving the incompressibility condition (161) β and μ being the stabilization parameter and shear modulus, respectively, and p 1 to p 3 denote the corresponding nodal pressure values of the mixed element. Application of the projection operator to the weak form of the constitutive equation involving the pressure term removes the instability arising when pressures and displacements are interpolated employing shape functions of the same order [6].

Appendix 2: Coupled thermomechanical formulation (mixed u/ p/θ )
In the following equations the strong and weak format of the coupled termomechanical problem for a single domain is outlined with the corresponding adopted constitutive relations and the staggered solution algorithm of the discretized FE quantities.

Strong form of the coupled thermomechanical problem
The strong form of the coupled thermomechanical problem neglecting the body forces can be expressed as FIND: FULFILLING: Equilibrium equation:∇ ∇ ∇ · P = 0, in n (187) Constitutive pressure equation: Energy balance: ρ 0U +∇ ∇ ∇ · Q = D int + ρ 0 R Dirichlet's boundary conditions: u =û, in u Neumann's boundary conditions: t =t, in σ , with ρ 0 being the initial density,∇ ∇ ∇ being the material Nabla operator containing derivatives with respect to the original configuration, P being the first Piola-Kirchhoff stress tensor, U is the internal energy, Q denotes the nominal heat flux vector and R stands for the internal generated heat.

Remark 7.1
Since incompressibility is assumed in all constitutive relations, the splitting of the Cauchy stress tensor in the current configuration is considered together with the constitutive pressure equation (188) in the remaining of the formulation.

Dissipation inequalities
The momentum and energy balance equations appearing in the strong form (187) to (191) are restricted by the the second law of thermodynamics from which the following dissipation inequalities can be written [35]: where D int and D con denote the internal dissipation and dissipation due to conduction, respectively and F stands for the deformation gradient tensor. The entropy is denoted by η and U and Q denote the internal energy and nominal heat flux, respectively. The expression in (192), also referred to as the Clausius Plank form of the second law of thermodynamics, indicates that the internal energy can never increase at the expense of decreasing stress power. The restriction of the dissipation by conduction in (193) indicates that the heat flux adopts the opposite sense of the temperature gradient which can be trivially verified by considering the expression of the heat flux according to Fourier's law: with k representing the positive thermal conductivity coefficient.

Thermomechanical model for J 2 flow theory
The following form of the free energy , i.e. the part of the internal energy U amenable to be transformed into work at constant temperature and defined as = U − θη, is considered for the adopted thermomechanical model at finite strains [35,36] be obtained by manipulating the constitutive entropy equation [11] and reads cθ = −∇ ∇ ∇Q + R + D mech − H, c being the specific heat at constant strain and constant internal variables per unit of reference volume and H stands for the structural elastoplastic heat [36].