Antonio Huerta
email: antonio.huerta@upc.edu

Enrique Nadal

| Francisco Chinesta

Proper Generalized Decomposition solutions within a Domain Decomposition strategy †

Keywords: Reduced Order Models, Proper Generalized Decomposition, Domain Decomposition, Parameterized solutions

Domain Decomposition strategies and the Proper Generalized Decomposition are efficiently combined to obtain a fast evaluation of the solution approximation in parameterized elliptic problems with complex geometries. The classical difficulties associated to the combination of layered domains with arbitrarily oriented midsurfaces, which may require in-plane-out-of-plane techniques, are now dismissed. More generally, solutions on large domains can now be confronted within a Domain Decomposition approach. This is done with a reduced cost in the offline phase. Because, the Proper Generalized Decomposition gives an explicit description of the solution in each subdomain in terms of the solution at the interface. Thus, the evaluation of the approximation in each subdomain is a simple function evaluation given the interface values (and the other problem parameters). The interface solution can be characterized by any a priori user-defined approximation. Here, for illustration purposes, hierarchical polynomials are used. The repetitiveness of the subdomains is exploited to reduce drastically the offline computational effort. The online phase requires to solve a nonlinear problem to determine all the interface solutions. But this problem only has degrees of freedom on the interfaces and the Jacobian matrix is explicitly determined. Obviously, other parameters characterizing the solution (material constants, external loads, geometry) can also be incorporated in the explicit description of the solution.

INTRODUCTION

Proper Generalized Decomposition (PGD) [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids[END_REF][START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF][START_REF] Chinesta | PGD for solving multidimensional and parametric models[END_REF][START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF] has proven its advantages and applicability in many parameterized problems. PGD reduces the computational complexity induced by a large number of dimensions (the sum of number of spatial dimensions plus the number of parameters) to the iterative resolution of low dimensional problems, usually 1-D or 2-D. Moreover, the PGD approach provides an explicit expression of the approximated solution. Thus, the online phase, that is the evaluation of the approximation given the parameters, is very fast and does not require any extra interpolation or another solve (typical, for instance, in Proper Orthogonal Decomposition strategies).

However, PGD has never been combined with Domain Decomposition strategies. The methodology proposed here makes it possible. The resolution is divided into two phases. The first one builds (offline) a PGD model (i.e. an explicit solution) of the Dirichlet-to-Neumann operator for each subdomain as a function of given parameters (materials, geometry...) and of the displacement/temperature along the interface. The imposition of displacement/temperature along the interface is done with Nitsche's method in order to grant the interface values its own discretization. Then, the second phase, the online one, imposes the weak balance of fluxes at the interface. It results in a small nonlinear system to be solved at the interface, which advantageously uses the explicit solutions in each subdomain to obtain the Jacobian matrix for the Newton-Raphson method. Since the PGD generalized approximation in each subdomain is explicit there are no subdomain solves.

This approach, which can be valuable in many sizable problems, is crucial when confronted to large problems with complex geometries. Splitting the original problem into different subdomains can be beneficial in many occasions. An obvious case occurs when the problem domain can be built as a non-overlapping union of a small set of subdomains geometrically similar, or when material parameters only affect one subdomain or are different from one subdomain to another, or with complex non-cartesian geometries. Recall that PGD separation in space works best in cartesian domains (generalized cubes) whereas its extensions to non-cartesian domains, although possible, see [START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF], may suffer as the complexity of the geometry increases.

Implementing a domain decomposition strategy with PGD is nonetheless mandatory when the domain is composed of subdomains where the PGD separation strategy is different from subdomain to subdomain. This is typical in problems obtained by the union of subdomains having non collinear local coordinates. For instance problems combine plates and shells with arbitrarily oriented mid-surfaces or mixed with solids areas.

In recent papers, PGD has shown its capabilities to solve shell (plate) structures with a 3D formulation at a 2D computational cost and without typical simplifications of shell (plate) elements, see [START_REF] Ghnatios | 3D Modeling of squeeze flows occurring in composite laminates[END_REF][START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF][START_REF] Bognet | Separated representations of 3D elastic solutions in shell geometries[END_REF]. They use an in-plane-out-of-plane strategy for the spatial PGD separation. This has even been extended to flow problems [START_REF] Ibáñez | Simulating squeeze flows in multiaxial laminates: towards fully 3D mixed formulations[END_REF]. Such advantages can only be exploited if all shells are aligned. Here, the formulation proposed overcomes these limitations.

Domain Decomposition (DD) techniques [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF] allow to overcome the difficulties of standard FE in order to: i) decrease the size of the problem, ii) combine different physics or iii) use different approximation methods in each subdomain. The subdomains generated with the DD technique are coupled via artificial boundaries. The problem coupling all subdomains can be solved with a wide range of alternatives [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF]. In this work, only non-overlapping methods are considered. Because, even the simplest implementation of overlapping DD with an in-plane-out-of-plane separated representation presents too many geometrical restrictions. In particular, the requirement of separability on the overlapping domains imposes major restrictions. A first attempt to combine PGD and DD strategies can be found in [START_REF] Nazeer | Arlequin based PGD domain decomposition[END_REF] where the Arlequin coupling strategy was considered. This reference shows the difficulties associated to the treatment of Lagrange multipliers, the geometrical transformation associated with general quadrilateral domains and the solution procedure for computing monolithically the in-plane and then the out-of-plane problems. Moreover, only homogeneous materials were considered to avoid the difficulties of arbitrarily oriented laminates and overlapping DD. For all these reasons, non-overlapping DD techniques are preferred in the presence of space separated representations for describing complex laminate structures.

In this context, it is possible to distinguish three different approaches: i) Dirichlet-Neumann [START_REF] Krause | A Dirichlet-Neumann type algorithm for contact problems with friction[END_REF], ii) Neumann-Neumann and iii) Dirichlet-Dirichlet methods. A popular variant of the latter is the Finite Element Tearing and Interconnecting (FETI) method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] or its extensions. The methodology proposed here is, however, closer to the concept introduced in [START_REF] Brezzi | A three-field domain decomposition method[END_REF] of interfaces with a discretizations independent of the subdomain and using Lagrange multipliers for coupling the subdomains. Here Nitsche's method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF][START_REF] Freud | On weakly imposed boundary conditions for second order problems[END_REF][START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF][START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] is used instead.

More recently, the idea of mixing DD techniques with Reduced Order Models (ROM) has show increasing interest. The aim being problems in which the domain can be built as a non-overlapping union of a small set subdomains geometrically similar [START_REF] Iapichino | A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks[END_REF][START_REF] Eftang | Port reduction in parametrized component static condensation: approximation and a posteriori error estimation[END_REF][START_REF] Eftang | A port-reduced static condensation reduced basis element method for large component-synthesized structures: approximation and A Posteriori error estimation[END_REF][START_REF] Vallaghé | The static condensation reduced basis element method for a mixed-mean conjugate heat exchanger model[END_REF][START_REF] Martini | Reduced basis approximation and a-posteriori error estimation for the coupled Stokes-Darcy system[END_REF][START_REF] Smetana | A new certification framework for the port reduced static condensation reduced basis element method[END_REF][START_REF] Smetana | Optimal local approximation spaces for component-based static condensation procedures[END_REF][START_REF] Iapichino | Reduced basis method and domain decomposition for elliptic problems in networks and complex parametrized geometries[END_REF]. In this sense, [START_REF] Iapichino | A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks[END_REF] proposes the Reduced Basis (RB) hybrid method to solve Stokes flow in a domain composed by different non-overlapping subdomains. These subdomains are grouped into geometrically similar subdomains, that is, subdomains with the same number of sides in 2D. For each set of similar subdomains a parametric RB technique is used. The parameters are geometrical ones. Therefore, by only generating a reduced basis for a small set of subdomain types, the method is able to generate a global reduced basis for the whole problem domain. This approach requires at the online phase to solve a global problem that includes the reduced basis (standard in this technique) plus a coarse finite element mesh (to capture the normal fluxes on the interfaces) and Lagrange multipliers (to impose the continuity of the primal function on the interface). This is further improved in [START_REF] Iapichino | Reduced basis method and domain decomposition for elliptic problems in networks and complex parametrized geometries[END_REF] by introducing (apart from the geometry) a set of parameters that characterize the profile of the function on the internal interfaces of the computational domain. This implies to solve, at the online phase, a system with has as many unknowns as parameters characterizing the solution on the interfaces plus all the reduced basis for each subdomain of the DD. In [START_REF] Iapichino | Reduced basis method and domain decomposition for elliptic problems in networks and complex parametrized geometries[END_REF] the proposed approach is further compared with other reduced basis strategies [START_REF] Maday | A reduced-basis element method[END_REF][START_REF] Iapichino | A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks[END_REF][START_REF] Huynh | A static condensation reduced basis element method: approximation and a posteriori error estimation[END_REF][START_REF] Martini | Reduced basis approximation and a-posteriori error estimation for the coupled Stokes-Darcy system[END_REF].

In this paper, in order to combine domains with different tensor product geometries, a DD approach within a PGD framework is proposed. Thus obtaining generalized solutions for each subdomain and reducing drastically the cost of the online phase. This approach is inspired by the original Schur complement method, which is considered the simplest and first version of nonoverlapping DD methods. The important feature is that it evaluates a parametric solution in each subdomain. The rationale proposed here requires first to obtain a parameterized PGD approximation for each set of similar subdomains. This offline phase is similar to [START_REF] Iapichino | A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks[END_REF][START_REF] Iapichino | Reduced basis method and domain decomposition for elliptic problems in networks and complex parametrized geometries[END_REF] but computes a generalized solution instead of a reduced basis. The parameters characterize the Dirichlet boundary conditions along the boundaries of these subdomains that are interfaces of the DD. These parameters describe the approximations used along the interface. Here for simplicity hierarchical polynomials are considered along the interfaces, without loss of generality and for illustration purposes.

Then, the online phase requires to solve a global problem, which in this case is a nonlinear system of equations, but only with unkowns along the interfaces. Moreover, the Newton-Raphson method can be readily applied because the Jacobian is easily computed. In fact, this global problem is established to guarantee equilibrium over the interfaces (transmission conditions for the fluxes). Note that the proposed approach allows for non-conforming discretizations between subdomains. More precisely, in the next sections it will be clear that this approach imposes weakly both the continuity of the solution and the equilibrium of the normal fluxes along the interfaces. Except for one numerical example, here, for simplicity, the parameters are restricted to the boundary conditions along the interface. But, in practice, this PGD approach can be extended easily to a more general framework with material or geometrical parameters [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF][START_REF] Ammar | Parametric solutions involving geometry: a step towards efficient shape optimization[END_REF][START_REF] Zlotnik | Proper generalized decomposition of a geometrically parametrized heat problem with geophysical applications[END_REF].

Note also that another advantage of the present approach is that the coupling problem is independent of the particular numerical method used in each subdomain. For instance, the proposed methodology allows to solve structures that are composed of degenerated subdomains, in which the space separation (in-plane-out-of-plane) is a must, and others in which that space separation is not required or even when in some subdomains the solution is governed by a nonlinear equation.

PROBLEM STATEMENT

The heat equation is used here for presentation purposes, because the same strategy can be readily extended to a more complex mechanical behavior. Moreover, also for presentation purposes, the geometry is simplified from the previously discussed shell reinforced structures to an open bounded domain in sd dimensions, that is Ω ⊂ ℝ sd , with boundary Ω split in the Dirichlet, Γ , and Neumann, Γ , parts. The strong form of the steady heat equation is well known:

⎧ ⎪ ⎨ ⎪ ⎩ -⋅ ( ) = in Ω, ⋅ ( ) = on Γ ⊂ Ω, = on Γ ∶= Ω ⧵ Γ , (1)
where is the unknown temperature, is the thermal conductivity matrix, the source/sink body term, the imposed heat flux on the Neumann boundary and the imposed temperature on the Dirichlet boundary. Note that Robin boundary conditions, which are typical in thermal problems, can also be considered. They are not included here to simplify the presentation.

Furthermore, Ω is partitioned in dd disjoint subdomains Ω with boundaries Ω , such that

Ω = dd ⋃ =1 Ω , Ω ∩ Ω = ∅ for ≠ ,
and the union of all the DD interfaces Γ , for = 1, … , fc , between the subdomains is denoted as The discontinuous setting induces a new problem equivalent to [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids[END_REF]. That is, the exact solution of (1) also satisfies the following problems,

fc ⋃ =1 Γ = Γ = dd ⋃ =1 Ω ⧵ Ω.
-⋅ ( ) = in Ω , ⋅ ( ) = on Ω ∩ Γ , = on Ω ∩ Γ , ⎫ ⎪ ⎬ ⎪ ⎭ for = 1, … , dd (2a)
= on Γ, (2b)
⋅ ( ) = 0 on Γ. (2c)
Note that the jump ⋅ operator has been introduced, it follows the definition in [START_REF] Montlaur | Discontinuous Galerkin methods for the Stokes equations using divergence-free approximations[END_REF]. That is, along each portion of the interface Γ it sums the values from the left and right of Γ (say, Ω and Ω ), namely

⊚ = ⊚ + ⊚ ,
and always must contain the normal vector .

THE DOMAIN DECOMPOSITION STRATEGY

The rationale proposed here is based on defining local problems for each subdomain and a global problem to determine the interface parameters characterizing each local problem. Note that the domain decomposition strategy presented here is not an algebraic domain decomposition strategy, but a physical-based strategy adapted for mechanical problems. The local problems will be solved offline using the PGD approach, as shown in the next section, and the global problem, which is nonlinear, is solved online. The unknowns of the global problem are new variables characterizing the primal unknown along the interfaces. This approach is described as a hybrid method [START_REF] Ciarlet | The finite element method for elliptic problems Classics in Applied Mathematics[END_REF] because a new variable is introduced on the interface, which coincides with the trace of the unknown. Figure 1 introduces schematically the notation used in what follows.

The local problem

By introducing a new variable ̂ along the interface Γ problem (2) becomes

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ -⋅ ( ) = in Ω , ⋅ ( ) = on Ω ∩ Γ , = on Ω ∩ Γ , = ̂ on Ω ∩ Γ, (3)
for = 1, … , dd , which accounts for (2a) and (2b), leaving the flux continuity imposed by (2c) for the global equation. Note that the set of problems defined in (3) are similar to those used by the first phase of a Neumann-Neumann based DD method [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF].

Solving problem [START_REF] Chinesta | PGD for solving multidimensional and parametric models[END_REF] in each subdomain determines the solution as a function of the unknown ̂ on the interface Γ. The weak problem associated to (3) is defined such that the Dirichlet boundary conditions are imposed weakly following a Nitsche rationale [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF][START_REF] Freud | On weakly imposed boundary conditions for second order problems[END_REF][START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF][START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF]. Thus, it consists in finding ∈  1 (Ω ), for all ∈  1 (Ω ) such that

∫ Ω ⋅ ( ) Ω -∫ Ω ⧵Γ ⋅ ( ̂ ) Γ = ∫ Ω Ω + ∫ Ω ∩Γ Γ,
where the traces of the fluxes, denoted by ̂ , are imposed as

⋅ ( ̂ ) = ⋅ ( ) -( -) on Ω ∩ Γ , ⋅ ( ) -( -̂ ) on Ω ∩ Γ,
and is a stabilization parameter, whose selection has an important effect on the stability and accuracy of the results and is proportional to 1∕ℎ where ℎ is the characteristic element size. More precisely, in the examples shown below, is taken as

= 2 max ∕ min ℎ .
That is, the maximum value of the diffusion coefficient at subdomain divided by the minimum characteristic size of the element at subdomain. Finally, the symmetrized local weak problem becomes: find ∈  1 (Ω ) such that

( , ) = ( ) ∀ ∈  1 (Ω ) (4a)
with

( , ) = ∫ Ω ⋅ ( ) Ω -∫ Ω ⧵Γ ⋅ ( ) Γ -∫ Ω ⧵Γ ⋅ ( ) Γ + ∫ Ω ⧵Γ Γ, (4b)
and

( ) = ∫ Ω Ω + ∫ Ω ∩Γ Γ -∫ Ω ∩Γ ⋅ ( ) Γ -∫ Ω ∩Γ ̂ ⋅ ( ) Γ + ∫ Ω ∩Γ Γ + ∫ Ω ∩Γ ̂ Γ. (4c)
It is important to remark that a similar strategy has been implemented imposing all Dirichlet boundary conditions (i.e. on Γ and Γ) as essential ones, see APPENDIX A:. However, imposing on Γ and ̂ along Γ as essential boundary conditions, in spite of being feasible, has two drawbacks. On one hand, as usual on PGD approaches, essential boundary conditions require and ad hoc initial mode. And, on the other hand, weak imposition of these conditions allow to use non conforming meshes on each side of the interface.

The interface problem

The discrete counterpart of local problem described in the previous section for each subdomain allows to evaluate approximations of as functions of ̂ . This unknown ̂ is determined with a interface problem or global problem, which is constructed imposing continuity of fluxes along the DD interface Γ, see (2c). These are the transmission conditions for the fluxes. The corresponding global weak problem can be stated as:

find ̂ ∈  2 (Γ) for all ∈  2 (Γ) such that dd ∑ =1 ∫ Ω ⧵ Ω ⋅ ( ) - Γ + ∫ Ω ⧵ Ω ̂ Γ = 0. (5)
In fact, the function ̂ ∈  2 (Γ) is defined by multiple subfunctions, with each subfunction applying to an interface of the domain decomposition, namely ̂ ( ) ∶= ̂ ( ) for ∈ Γ ,

and ̂ ∈  2 (Γ ) for = 1, … , fc .

THE PGD DOMAIN DECOMPOSITION STRATEGY

The previous section allows to split the original problem into many local problems, which determine the primal solution as a function of the interface ̂ , and one global problem to solve precisely for ̂ . As usual, the unknow function ̂ is characterized by a set of coefficients (û 1 , … , û ̂ ) = ̂ . Consequently, each local problem is now a parameterized problem and the parameters are precisely the coefficients in ̂ . The PGD strategy proposed here, corresponding to the classical PGD formulation, solves all the parameterized local problems in the offline phase in order to obtain the generalized solution in terms of ̂ , i.e. the Dirichlet boundary conditios as parameters. In case of similar subdomains, the number of offline PGD problems is always much smaller than the number of subdomains, see the examples in Section 5. Subsections 4.1 to 4.3 describe how to define and obtain generalized solutions for each subdomain during the offline phase of the PGD strategy. Then Subsection 4.4 describes the global problem that implies solving a nonlinear system of equations to determine ̂ . Note that the Jacobian matrix needed to implement a Newton-Raphson scheme can be easily computed. This global solve is the only computational burden of the online phase because once ̂ is determined the primal variable is determined on-the-fly given the explicit expressions for each subdomain obtained by the PGD.

The parameterized subdomain problem

To proceed with this strategy a parametric characterization of ̂ along Γ is required, for instance

̂ ( ) = ̂ ∑ =1 ̂ ( ) û = ̂ ( ) ̂ , for ∈ Γ, (6)
where ̂ ∈ ℝ ̂ are the coefficients and functions ̂ ( ) are the basis used to represent the solution along the interfaces. As noted earlier, this approximation is constructed from the subfunctions on each interface Γ , for = 1, … , fc . In fact, for each interface

̂ ( ) = ̂ ∑ =1 ̂ ( ) û for ∈ Γ , (7)
where ̂ are hierarchic shape functions derived from the set of integrated Legendre polynomials as introduced in [START_REF] Szabó | Finite element analysis: formulation, verification, and validation[END_REF]. To simplify the presentation the approximation degree in every interface is uniform ̂ = ̂ . Consequently, ̂ = ̂ fc but, as commented in Remark 1, this is not necessary. The parametrization of ̂ defines a finite dimensional space  ⊂  2 (Γ), with dim  = ̂ . The global problem ( 5) is solved precisely in . But this space is completely independent of the trial spaces  ⊂  1 (Ω ) used in each local problem (each subdomain). That is, ̂ ( ) is not related to the shape functions of the mesh in Ω for = 1, … , dd . And consequently, nonconforming meshes between subdomains are naturally taken into consideration.

The strategy proposed here, which follows precisely the computational vademecum rationale proposed in [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF], is to solve each local problem (4) using parameters ̂ , which characterize ̂ , as new coordinates. Thus, the objective is to compute ( , ̂ ) with ( , ̂ ) ∈ Ω ×  ⊂ ℝ sd × ℝ ̂ . Where  ⊂ ℝ ̂ characterizes the range for parameters ̂ and can be defined as the Cartesian combination of the range for each parameter, namely,  ∶=  1 ×  2 × ⋯ ×  ̂ with û ∈  for = 1, … , ̂ . This range is user-defined depending on the expected values, a priori physical knowledge of the problem (for instance, avoiding non-physical temperatures), and the imposed boundary conditions (and source terms).

The solution ( , ̂ ) formally lies in a tensor product space, namely,

∈  ⊗  2 ( 1 ) ⊗  2 ( 2 ) ⊗ ⋯ ⊗  2 ( ̂ ).
A standard weighted residuals approach, with integrals in Ω ×  and the usual integration by parts only in Ω produces a weak form in this multi-dimensional setup.

For simplicity, in order to avoid cumbersome notation, the subdomain indices will not be explicitly indicated unless it is crucial to understand the development. The weak problem for each subdomain becomes: find

∈  ⊗  2 ( 1 ) ⊗  2 ( 2 ) ⊗ ⋯ ⊗  2 ( ̂ ) such that ( , ) = ( ), ∀ ∈  ⊗  2 ( 1 ) ⊗  2 ( 2 ) ⊗ ⋯ ⊗  2 ( ̂ ), (8a)
with the following definitions of the bilinear and linear forms

( , ) ∶= ∫  1 ∫  2 ⋯ ∫  ̂ ( , ) û ̂ ⋯ û2 û1 and ( ) ∶= ∫  1 ∫  2 ⋯ ∫  ̂ ( ) û ̂ ⋯ û2 û1 , ( 8b
)
that require the previously defined spatial forms (4).

It is important to recall that the resolution of (8a) in a PGD context requires an affine parameter dependence of the different forms. This as usual in reduced order methods and it is very well discussed in [START_REF] Patera | Reduced Basis Approximation and A-Posteriori Error Estimation for Parametrized Partial Differential Equations[END_REF][START_REF] Rozza | Fundamentals of reduced basis method for problems governed by parametrized PDEs and applications[END_REF]. In fact, finding the affine parameter dependence of (8b) is discussed in subsequent sections.

The PGD approximation in each subdomain

As usual in a PGD strategy, a separated representation PGD is imposed to approximate the solution of (8) in each subdomain, namely

( , ̂ ) ≈ PGD ( , ̂ ) = ∑ =1 1 ( 1 ) 2 ( 2 ) ̂ ∏ =1 ̂ û (û ) = -1 PGD ( , ̂ ) + 1 ( 1 ) 2 ( 2 ) ̂ ∏ =1 ̂ û (û ) (9)
The function PGD stands for the PGD approximation with terms (or modes) and is defined as a sum of separable terms. Each term (mode) is the product of functions depending only on one of the arguments. This separation as a sum of products of unknown functions is standard in PGD, see for instance [START_REF] Chinesta | PGD for solving multidimensional and parametric models[END_REF][START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF]. Obviously, this implies that the unknown interface parameters are arguments of unknown functions ̂ û . Note that in some of the PGD implementations the separated modal functions are normalized and therefore a scalar coefficient affects each mode and characterizes its amplitude.

Moreover, it is important to note that expression (9) already considers a spatial separated representation of the solution. This space separation can precisely coincide with an in-plane ( 1 ) out-of-plane ( 2 ) strategy, which is useful when the geometry of the corresponding subdomain is a degenerated geometry such as a laminate [START_REF] Bognet | Separated representations of 3D elastic solutions in shell geometries[END_REF]. Of course, all developments that follow in this section can also be applied when no space separation is required.

As usual, the separated approximation proposed by ( 9) is determined using a greedy algorithm. This assumes that -1 PGD is known and the last mode,

1 ( 1 ) 2 ( 2 ) ∏ ̂ =1 ̂ û (û )
, is computed replacing ( 9) into [START_REF] Bognet | Separated representations of 3D elastic solutions in shell geometries[END_REF]. Then, the fixed-point solver for this nonlinear problem is applied, details can be found elsewhere see for instance [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF][START_REF] Zlotnik | Proper generalized decomposition of a geometrically parametrized heat problem with geophysical applications[END_REF].

Affine parameter dependence

As indicated previously, the greedy strategy implies replacing ( 9) into (8) to determine the last mode. This entails solving a nonlinear problem to compute

1 ∈  1 , 2 ∈  2 , ̂ û1 ∈  2 ( 1 ), … , ̂ û ̂ ∈  2 ( ̂ ) for all 1 ∈  1 , 2 ∈  2 , ̂ û1 ∈  2 ( 1 ), … , ̂ û ̂ ∈  2 ( ̂ ) such that 1 2 ̂ ∏ =1 ̂ û , = ( ) -( -1 PGD , ) with = 1 2 ̂ ∏ =1 ̂ û + 1 2 ̂ ∏ =1 ̂ û + ̂ ∑ =1 1 2 ̂ û ̂ ∏ =1 ≠ ̂ û .
It is important to realize that the definition of these forms, see (8b), and the separated structure of the approximation allow to clearly see the affine parameter dependence, namely

1 2 ̂ ∏ =1 ̂ û , = ∫ 1 ̂ û1 ∫ 2 ̂ û2 ⋯ ∫ ̂ ̂ û ̂ 1 2 , û1 ⋯ û ̂ , ( ) = ∫ 1 ∫ 2 ⋯ ∫ ̂ ̃ ( ) + ̂ ∑ =1 û ̂ ( ) û1 ⋯ û ̂ , ̃ ( ) = ∫ Ω Ω + ∫ Ω ∩Γ Γ -∫ Ω ∩Γ ⋅ ( ) Γ + ∫ Ω ∩Γ Γ, ̂ ( ) = ∫ Ω ∩Γ ̂ Γ -∫ Ω ∩Γ ̂ ⋅ ( ) Γ.
The details on the greedy approach and the alternate direction method to solve this nonlinear problem can be found elsewhere, see for instance [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF][START_REF] Chinesta | PGD for solving multidimensional and parametric models[END_REF][START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF][START_REF] Zlotnik | Proper generalized decomposition of a geometrically parametrized heat problem with geophysical applications[END_REF].

Implementation aspects of the global problem

The generalized solution of every local problem provides an explicit expression of the solution in each subdomain for any value of the primal variable on the interfaces. Recall that only a small set of subdomains are actually solved, only those that are not similar and present variations in topology, boundary conditions on non-interface borders, different in-plane-out-of-plane strategies, etc.

A trivial Schur complement strategy can now be implemented by means of imposing the flux transmission conditions on the interface Γ. First, the global problem defined in (5) must be solved. As noted in Section 4.1, the parametrization of ̂ defines a finite dimensional space  ⊂  2 (Γ), with dim  = ̂ . Since the global problem ( 5) is solved in , the test functions are in  = span{ ̂ 1 , … , ̂ ̂ }. Consequently, the nonlinear system of ̂ equations and ̂ ∈ ℝ ̂ unknowns is obtained after replacing the approximation PGD ( , ̂ ) ≈ ( , ̂ ) in (5), namely

(̂ ) = dd ∑ =1 ∫ Ω ⧵ Ω ̂ ( ) ⋅ PGD ( , ̂ ) - PGD ( , ̂ ) Γ + ∫ Ω ⧵ Ω ̂ ( ) ̂ ( ) ̂ Γ = . ( 10
)
Note that the nonlinearity comes from the fact that the expression of PGD ( , ̂ ), see [START_REF] Ibáñez | Simulating squeeze flows in multiaxial laminates: towards fully 3D mixed formulations[END_REF], involves sums of products of functions,

̂ û (û ), of the unknowns ̂ = [û 1 , … , û ̂ ] .
Second, it is important to note that once the local problems are solved, the following expression can be readily computed

PGD ( , ̂ ) û = ∑ =1 1 ( 1 ) 2 ( 2 )
̂ û (û ) û ̂ ∏ =1 ≠ ̂ û (û ) .
And, consequently, the chain rule allows to explicitly compute the Jacobian matrix needed to implement a Newton-Raphson scheme for solving [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF]. Finally, it is important to observe that in all the linear examples discussed in the following sections the number of iterations never exceeded three. The quadratic convergence of the Newton-Raphson method allows to obtain the desired accuracy in this reduced number of iterations.

Remark 1. Hierarchical shape functions. Here, a polynomial basis is used to describe ̂ along each interface, see [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF]. Obviously other options are possible, for instance, piecewise polynomials or any other user-defined linear approximation along each interface. To show the applicability of the approach, without loss of generality, polynomials are used along each interface moreover to simplify the presentation the approximation degree in every interface is uniform. In fact, the user can decide, using or not a priori information on the problem, to use continuous functions with discontinuous derivatives on heterogeneous interfaces or to adapt, on each interface in Γ, the degree of the continuous polynomial depending on the problem at hand. See [START_REF] Smetana | Optimal local approximation spaces for component-based static condensation procedures[END_REF] for more details.

NUMERICAL EXAMPLES

Academic example with two subdomains

This is a simple example where PGD is not needed and, consequently, applying the PGD with domain decomposition is unnecessary. It is used to check convergence to the FEM solution as the approximation on the interface is enriched. The material is considered isotropic = , where is the identity matrix.The diffusion coefficient in the domain is = 1. The computational domain Ω is partitioned in Ω = [ -2, -1] × [-1, 1] for = 1, 2. Dirichlet and Neumann boundary conditions are imposed according to the schematic presentation shown in Figure 2 . The prescribed source term in (2) is defined as

( 1 , 2 ) = 50 exp -10( 1 -0.5) 2 -10 2 2 for ( 1 , 2 ) ∈ [-1, 1] × [-1, 1].
A uniform FEM mesh of 60 × 30 or 120 × 60 bilinear elements is used in complete domain Ω as a reference. For comparison, the PGD separation in and uses linear elements of the same characteristic size in each subdomain. Therefore, in each subdomain a mesh of 30 × 30 or 60 × 60 elements is used. Along the interface Γ, as noted earlier, a normalized hierarchical polynomial basis is used to describe the unknown ̂ , see [START_REF] Ghnatios | 3D Modeling of squeeze flows occurring in composite laminates[END_REF]. The parameters û ( = 1, … , ̂ ) are allowed to range in  ∶= [-10, 10]. Each domain  is discretized in 100 uniform elements. Figure 3 shows the results obtained with the proposed method with two different approximations along the interface: constant, i.e. ̂ = 1, (left) and a 6th order hierarchical polynomial approximation, i.e. ̂ = 7, (right). Note that, due to the symmetry of the solution, only even-order polynomial are used (the odd ones behave exactly as the corresponding lower even approximation). The spatial mesh of 120 × 60 bilinear elements is used in the whole domain for the FEM reference solution. It corresponds to discretizations of 60 linear elements in the PGD separation. Figure 4 shows error maps between the Domain Decomposition PGD and the FEM approximation for the two interface parameterizations (constant and polynomial of degree 6). The absolute error is normalized by the maximum value of the solution in the domain (approximately 4.39). It is important to note that the FEM approximation on the 120×60 bilinear mesh is used to evaluate the error and it corresponds to an equivalent spatial discretization. Thus, this error map includes both the standard PGD truncation error, see [START_REF] Zlotnik | Effect of the separated approximation of input data in the accuracy of the resulting PGD solution[END_REF], and the error induced by the approximation imposed along the interface. The second one is due to approximation of ̂ along the interface and must be controlled. In fact, the influence of this error can be clearly seen comparing the two images in Figure 4 . Note that a crude constant approximation for ̂ induces relative errors below 13.1%, whereas a 6th order polynomial approximation reduces the maximum relative error to 1.86%.

The computational cost of the offline and online phases and the number of modes required are shown in Table 1 . The offline phase implies solving two PGD problems since the source term is different in each subdomain. Since the offline computations are independent for each subdomain, the offline computational cost corresponds to the maximum cost for each subdomain.

The discretization is 60×60 elements per subdomain. Obviously, increasing the polynomial degree along the interface implies more modes of the PGD solution and more offline computational cost. The online cost is associated to the generation of the Jacobian matrix needed to implement a Newton-Raphson scheme for solving [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF] and the number of iterations. Obviously, this problem, which does not need a PGD separation in the subdomains and the decomposition, can be solved with FEM on a 120×60 mesh very fast, namely 0.077 . However, even in this case, the online computational cost is smaller than the FEM computation. Next examples will show situations where the online solve is orders of magnitude smaller than the FEM resolution.

In fact, the influence of the polynomial order is further analyzed in Figure 5 where the  2 relative error along the interface is depicted as the polynomial order in the approximation of ̂ is increased, at left. At right, the  2 relative error of the flux is also evaluated, considering the flux outwards subdomain Ω 1 . The reference solution here is the one obtained with an overkilled FEM calculation on a mesh of 600 × 300 bilinear elements over the whole domain Ω. The dotted and solid lines correspond to DD PGD computations with spatial meshes of characteristic size ℎ = 1∕30 and ℎ = 1∕60 respectively. Additionally the influence of the truncation of the greedy algorithm used for the local (PGD) problem is analysed through the squared and diamond markers. The truncation tolerance is set to 10 -4 and 10 -5 , and is evaluated as the 2-norm of the last mode with respect to the first one. In the rest of the problems this tolerance is set to 10 -5 .

There are three sources of error included in the curves of Figure 5 . The first two are standard in PGD approximations and correspond to the FE error and the PGD truncation, see [START_REF] Zlotnik | Effect of the separated approximation of input data in the accuracy of the resulting PGD solution[END_REF]. The standard FE error is (ℎ 2 ), namely (1∕30 2 ) ≈ 10 -3 and (1∕60 2 ) ≈ 3 10 -4 respectively. Note that it includes the weak imposition of the Dirichlet boundary condition on the interface using Nitsche's approach. The third source of error is associated to the approximation of ̂ along the interface. The temperature on the interface instead of being piecewise linear (with the corresponding characteristic size of ℎ = 1∕30 or ℎ = 1∕60) is approximated by a polynomial of degree ̂ -1. The error of approximating a piecewise linear function by a polynomial of increasing order is bounded below by (ℎ 2 ). Thus, the saturation observed in Figure 5 , which is one order of magnitude larger than ℎ 2 , is due to the PGD truncation and can be reduced by increasing the number of PGD modes in [START_REF] Ibáñez | Simulating squeeze flows in multiaxial laminates: towards fully 3D mixed formulations[END_REF] and the proper approximation of input data, see [START_REF] Zlotnik | Effect of the separated approximation of input data in the accuracy of the resulting PGD solution[END_REF]. As you can appreciate, the increase in the accuracy when a smaller truncation error is considered, i.e. increasing the number of PGD modes, has only a marginal effect. A considerable higher number of modes should be added to be able appreciate the influence. However, this is not required since the accuracy levels that the proposed methodology shows are enough for engineering applications. On the other hand, the right graph in Figure 5 shows that the accuracy in the fluxes keeps increasing when the degree of approximation increases along the interface. As a difference from the solution error, the error in the flux does not seem to depend on the mesh size until a certain level of approximation, i.e. degree 6, is reached. Finally, it is important to recall that as noted in Section 4.4 a nonlinear global problem is solved to determine the parameters characterizing the approximation of ̂ along the interface. Quadratic convergence of the residual of the Newton-Raphson solver is always obtained since the Jacobian can be explicitly evaluated. Moreover, only three iterations are enough to obtain convergence (residual smaller than 10 -8 ). In fact, the next examples show that the convergence of the Newton-Raphson scheme is independent of the number of parameters and subdomains or interfaces.

Layered academic example with two subdomains

The previous example (same source term and boundary conditions) is repeated here but with layered subdomains. Two cases are considered depending if layers are collinear or orthogonal in both subdomains, see Figure 6 . Each subdomain is composed by laminates: 5 fiber layers (black ones) and 6 matrix layers (gray ones). The thermal conductivity of the fiber layers ( = 10), is one order of magnitude larger than the one for the matrix layers ( = 1). The layered structure suggests for each subdomain an in-plane-out-of-plane formulation with 110 linear elements in the direction across the laminates and only 22 linear element in the layer direction. Note that the nonconformity along the interface for the example with orthogonal layers is not a problem in the proposed formulation. Because, as noted in section 2, along the interface temperature is imposed in weak form. Thus, the approach proposed here allows for non-conforming meshes at the interfaces. However, when the standard FEM is used in the whole domain for the orthogonal case, in order to have conformity along the interface Γ the mesh on the left subdomain must have also 110 linear element in the layer direction.

The influence of the laminates can be appreciated in Figure 7 . This figure shows the solutions obtained with a polynomial of degree 6 along the interface Γ. The in-plane-out-of-plane strategy combined with the domain decomposition produces reasonable results when compared with FE computations on an overkilled mesh, always conforming, of 2200 × 1100 bilinear elements. Note however from Figure 8 that as the degree on the interface increases the error along the interface decreases in the collinear case but remains constant in the orthogonal one. In this figure, the error is measured with respect to the previously indicated overkilled mesh. As expected, the hierarchical polynomials employed along the interface are not able to reproduce along the interface Γ the slope discontinuities of the temperature field across the laminates. Because as shown in Figure 9 double discontinuities (in both Cartesian directions) for the orthogonal layered case are difficult to capture by continuous polynomials. Consequently, if the interface area is of interest other approximations along the interface, such as piecewise polynomials, could be used instead of hierarchical shape functions.

In order to further analyze the quality of the proposed approach, error maps with respect to equivalent FE computations are plotted in Figure 10 . The equivalent FE computations require conforming meshes, thus a mesh of 110 × 220 bilinear elements is used as reference in Figure 10 . As expected from Figure 8 , the maximum error at the interface is not influenced as the polynomial degree is increased. However, as the order of the polynomial approximation on the interface increases, the error is clearly more localized and decreases away from the interface Γ.

Obviously, this issue is not present in the case of collinear laminates. Figure 11 shows the error with respect to an equivalent FE mesh of 110×44 bilinear elements. It can be observed that the error decreases drastically everywhere as the polynomial order increases on the interface. Finally, for the collinear configuration, Table 2 shows computational times for the offline, which implies solving two PGD problems because as noted in the previous example the source term is different in each subdomain (recall that only the maximum computational cost for each subdomain is indicated). Obviously, increasing the polynomial degree along the interface implies more modes of the PGD solution and more offline computational cost. As noted earlier, the online cost is associated to the generation of the Jacobian matrix needed to implement a Newton-Raphson scheme and the number of iterations. Evidently, this academic problem can also be solved with FEM on a 110 × 44 mesh very fast, namely 0.06 . As in the previous example, the online computational cost of the proposed technique is always smaller than the FEM computation. In summary, the domain decomposition PGD strategy allows combining subdomains with different in-plane-out-of-plane orientations. The lack of mesh conformity is not a problem for this formulation and the approximation along the interface must be oriented to the desired quantities of interest.

Composite thermal fin problem

This problem, inspired by [START_REF] Maday | The Reduced Basis Element Method: Application to a Thermal Fin Problem[END_REF][START_REF] Iapichino | Reduced basis method and domain decomposition for elliptic problems in networks and complex parametrized geometries[END_REF], models the heat diffusion in the fin depicted in Figure 12 with Dirichlet imposed temperatures at Γ in ( = 1) and Γ out ( = 0), and homogeneous Neumann elsewhere. In contrast with the cited previous references, here the material is composed by laminates with different orientations in each part of the geometry, see Figure 12 . As in the previous example, these laminates are composed of 5 fiber layers (black ones) and 6 matrix layers (gray ones). The thermal conductivity of the fiber layers ( = 10), is one order of magnitude larger than the one for the matrix layers ( = 1). To capture the features of the solution in the direction perpendicular to the laminate, an in-plane-out-of-plane strategy is used. Ten elements are considered for each layer, i.e. a total of 110 elements are used in the out-of-plane direction. Whereas, only 22 elements are employed along the in-plane direction.

The interface Γ, represented by shadowed lines, is now the union of all 15 sides joining these subdomains. In each one of these sides ̂ is approximated by a polynomial of degree one or three. Given the temperatures imposed at Γ in and Γ out , the parameters characterizing temperature at the interfaces are taken in [-2, 2]. All these 1D domains for each parameter are discretized with a uniform mesh of 100 linear elements.

As observed in Figure 12 , all the fins are geometrically similar. In fact, all of then can be transformed into a unit square with a simple scaling of each direction, see Figure 13 . The physical coordinates ̃ are scaled by the physical dimensions of the subdomain ( 1 , 2 ) = ( ̃ 1 ∕ 1 , ̃ 2 ∕ 2 ). Consequently, the thermal problem that originally was defined in the physical For this problem the online computational cost is orders of magnitude smaller. Note that for this problem the variation of the dimensions of any subdomain will only require to perform another online evaluation and not the full process, considerably reducing now the computational cost with respect to FEM since the offline computation is performed only once. The extension for any other parameters is straightforward. degree offline ( ) online ( ) modes Sub: 1/2 3 Computational cost on a laptop and number of modes for the composite thermal fin example with different degrees of approximation along the interface. A mesh of 22 × 110 elements per subdomain Ω is used.

FIGURE 16

Laminate skin and stringers. Model of the problem and boundary conditions. The interface between subdomains is represented by a shadowed line.

Laminate skin and stringers

This problem is devoted to show the capabilities of the proposed approach to a more complex and realistic geometry, see Figure 16 . The in-plane-out-of-plane PGD approach is need for the skin but adding the stringers precludes its use under classical implementations.

As shown in Figure 16 the domain is decomposed in 24 subdomains and its corresponding 24 interfaces (represented by shadowed lines). As in the previous example, the exterior annulus is formed by 11 layers (5 fibers depicted in black and 6 matrix in gray). The thermal conductivity in the fiber layers ( = 10) is one order of magnitude larger than the thermal ones ( = 1).

Given the geometry, the boundary conditions and the number of interfaces for each subdomain, five different typologies are identified. This imply that the PGD will only be applied to five different problems. These typologies are indicated in Figure 16 . Typologies 1 (skin connected to stringer), 2 (skin subjected to the external heat source on Γ source ) and 3 (skin with both external boundaries corresponding to Γ warm ) correspond to the skin, whereas 4 (stringer with one boundary Γ hot ) and 5 (stringer with one boundary Γ cold ) are associated to the stringers. The prescribed temperature at each boundary is

= ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 0 on Γ cold , 1 on Γ hot , 1∕2 on Γ warm , (1 -2 ) exp(-200 2 ∕9) + 1∕2 on Γ source ,
where ∈ [-1, 1] is the curvilinear coordinate along Γ source (shown as a dotted line in Figure 16 ). For subdomains corresponding to the skin (typologies 1, 2 and 3) and in order to capture the richness of the solution in the radial direction 110 1D linear elements are employed (10 elements per layer). Whereas only 20 1D linear elements are used along the circumferential direction for each subdomain. Stringers are discretized with a standard FE linear triangular elements for a total of 1105. The parametric space for each one of the parameters is Ω = [-2, 2] and it is discretized with 100 linear elements.

In this problem the temperature field along the interfaces is approximated linearly and quadratically. Figure 17 shows the PGD solution considering a quadratic approximation along the interface between all subdomains. The temperature field for the linear case is not shown because it is difficult to see any differences. But the improvement of the quadratic approximation is clearly seen if the error map is presented, see Figure 18 . These maps depict the error (normalized by the temperature imposed at Γ hot , i.e. = 1) between the PGD approximation and a FE computation on an equivalent mesh. For the linear approximation, the error is concentrated close to the right stringer between elements of typology 1 and 2 (8.6%). But just one extra degree of freedom on the interfaces reduces considerably the error (below 3.5%), which is now localized at the interfaces between stringers and annulus. As introduced in the previous examples, the convergence of the nonlinear solver used for the global problem along the interfaces does not depend on the number of interfaces neither the number of parameters used to define the solution along the interfaces since it converges in only three iterations for a tolerance lower than 10 -8 .

Finally, table 4 shows the computational time of the offline and online phases and the number of modes considered for the different subdomain types, considering a mesh of 20 × 110 elements per subdomain, with linear and quadratic approximation along the interfaces. As in the previous case, the equivalent FEM problems is also solved using an equivalent mesh with a computational cost of 43 . As in the previous problems, the online computational cost is smaller than the FEM computational cost.

Nonlinear problem

In this section we extend the proposed methodology to material nonlinear models. To illustrate its capabilities, the first academic example is revisited, see Section 5.1 and Figure 2 . The thermal source term is ignored but now the thermal conductivity is dependent on the temperature, in particular, it is assumed that

( ) = 10 + 7 ,
where is the temperature field. This nonlinearity does not affect to the general methodology. The only difference is that the local parametric solution is obtained by solving the corresponding nonlinear parametric problem with the PGD. Among the different choices for addressing the nonlinearity, the simplest option is considered here: a fixed point strategy. In this case, the affine parameter dependence of the bilinear form is ensured because the conductivity can always be expressed in separated form given the polynomial dependence of and the separated expression for the temperature ≈ -1 PGD . Thus, the nonlinear solution does not implies a major difficulty.

Along the interface a linear approximation of the temperature field is considered. From the a-priori knowledge provided by the linear model solution discussed in Section 5.1, the parametric space for each one of the parameters was chosen in  ∶= [-0.5, 2.5]. A uniform mesh of 100 linear elements is used in each of theirs segments. The space coordinates are discretized using 30 linear elements in each direction.

Figure 19 shows the solution obtained with the proposed method. The influence of the material nonlinearity induces a slight qualitative difference with respect to the linear case where the temperature field is expected to increase linearly along the 1 -coordinate. Figure 19 also shows the error with respect to the FE solution using an equivalent mesh and the same nonlinearity. The error on the interface is almost negligible (because the approximation can exactly represent the exact solution on the interface), and the error in each subdomain is mainly due to the truncation of the finite sum decomposition.

Finally, figure 20 shows the convergence of the Newton-Raphson strategy used for solving the global problem at the interfaces. In this case, and despite the fact of having a final solution even more simple that the solution related to the linear problem (having neglected here the localized thermal source term) 5 iterations are needed to reach the same precision that the one considered in the linear case. This fact proves that now, the precomputed parametric solution have a more complex dependence with respect to the parameters related to the boundary conditions.

CONCLUSIONS

A new approach has been proposed to obtain an explicit description of the solution in parameterized elliptic problems with complex geometries. In this context, complex geometries include multi-domains decomposed in arbitrarily oriented repetitive topologies. The proposed approach allows to combine different formulations, among them and as shown in the examples, standard finite elements in non-conforming subdomains and/or in-plane-out-of-plane techniques for layered domains with arbitrarily oriented mid-surfaces. For this purpose, domain decomposition strategies and the Proper Generalized Decomposition are combined in an efficient manner. The repetitiveness of the subdomains is exploited to reduce drastically the offline computational effort. Moreover, an efficient global domain decomposition solve is also included in this offline phase to drastically reduce the cost of the online phase, which is now reduced to the explicit evaluation of the solution with no additional problem solving.

The numerical examples show the potential of the proposed methodology for practical engineering applications.

FIGURE 1

 1 FIGURE 1 Splitting of the whole domain Ω into two domains Ω 1 and Ω 2 .

FIGURE 2

 2 FIGURE 2 Academic example model problem and boundary conditions with an interface between two subdomains.

4 FIGURE 3

 43 FIGURE 3 Academic example solution with a constant (left) and order 6th (rigth) approximation of ̂ along the interface. 60 × 60 elements per subdomain are used.

3 FIGURE 4 TABLE 1

 341 FIGURE 4 Academic example relative error distribution (with respect to the 120 × 60 bilinear FEM) for a constant (left) and order 6th (right) approximation of ̂ along the interface. Note that the scales are different. 60 × 60 elements per subdomain are used.

30 FIGURE 5

 305 FIGURE 5Academic example  2 relative error, with respect to a reference FEM solution obtained with 600 × 300 elements in Ω, along the interface versus polynomial order in the approximation of ̂ .

FIGURE 6

 6 FIGURE[START_REF] Ghnatios | 3D Modeling of squeeze flows occurring in composite laminates[END_REF] Academic example model with two layered subdomains. The interface between subdomains is in the middle and is represented by a grey vertical line.

2 FIGURE 7 FIGURE 8

 278 FIGURE 7 Solution for the collinear (left) and orthogonal (right) laminates.

1 FIGURE 9

 19 FIGURE 9 Layered subdomains with collinear (left) and orthogonal (orthogonal) layers FE reference solution.

2 FIGURE 10

 210 FIGURE 10Orthogonal layered subdomains error distribution for a constant (left) and order 6th (right) approximation of ̂ along the interface.

25 FIGURE 11 TABLE 2

 25112 FIGURE[START_REF] Nazeer | Arlequin based PGD domain decomposition[END_REF] Collinear layered subdomains error distribution for a constant (left) and order 6th (right) approximation of ̂ along the interface.

FIGURE 15

 15 FIGURE[START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] Composite thermal fin problem. Error map with respect to a FE mesh for a linear (top) and cubic (bottom) approximation of the solution along the interfaces.

FIGURE 17

 17 FIGURE 17 Laminate skin and stringers. Solution of the problem for a quadratic approximation of the solution along the interfaces.

FIGURE 18

 18 FIGURE[START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] Laminate skin and stringers. Error map with respect to a FE mesh for a linear (top) and quadratic (bottom) approximation of the solution along the interfaces.

3 FIGURE 19

 319 FIGURE[START_REF] Iapichino | A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks[END_REF] Nonlinear problem. DD-PGD solution (left) and relative error with respect to the FEM solution for a mesh of 30 × 30 bilinear elements (right).

FIGURE 20

 20 FIGURE[START_REF] Eftang | Port reduction in parametrized component static condensation: approximation and a posteriori error estimation[END_REF] Convergence of the Newton-Raphson strategy used for solving the global problem in the nonlinear thermal heat equation.

TABLE 4

 4 Computational cost on a laptop and number of modes for the laminate skin and stringers example with different degrees of approximation along the interface. A mesh of 20 × 110 elements per subdomain Ω is used.

	3/4/5

FIGURE 12

Composite thermal fin problem statement.

FIGURE [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] Scaling the physical domain into the unit square. rectangle with layers in the ̃ 1 direction an a thermal conductivity = ( ̃ 2 ) is now solved in the unit square (reference domain ]0, 1[×]0, 1[) with a thermal conductivity

Obviously, for each subdomain the scaling may be changed. In summary, the distinct geometrical domains are solved in the unit square as an orthotropic material, and, consequently, the geometry is simply parametrized as material constants. Thus, the PGD approximation in each subdomain includes now as parameter the dimensions 1 and 2 , namely

The dimensions of the subdomains are allowed to range two orders of magnitude, that is ∈  ∶= [0.1, 10] for = 1, 2. Each of these parametric domains  are discretized with a uniform mesh of 100 linear elements. Details on the PGD resolution of parameterized problems can be found, for instance, in [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF][START_REF] Chinesta | PGD for solving multidimensional and parametric models[END_REF][START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF]. Two cases are studied depending on the approximation along the interfaces: linear ( ̂ = 2) or cubic ( ̂ = 4). As observed earlier, here for simplicity the same approximation is used along all the interfaces but nothing precludes using different approximations along each interface depending on the richness of the temperature field. In fact, the hierarchical polynomials allow for an easy implementation of this approach.

Only the temperature field for the cubic interface approximation is depicted in Figure 14 , because the differences with the linear case are not distinguishable. Note that, in this problem the PGD approximation includes as extra parameters the dimensions of the subdomain (to account for the geometry variations) apart from the temperature on the interfaces.

The error map for the linear and cubic interface approximations are shown in Figure 15 . These errors are normalized by the maximum value of the solution (precisely = 1 at Γ in ) and computed with respect to a FE solution on a mesh with the same resolution in the out-of-plane direction (orthogonal to the layers) and an overkill discretization in the in-plane-direction. This is done to avoid, in the FE computation, the use of nonconforming meshes along the interfaces where fibers orthogonal on each side. Note that the DD PGD proposed can handle these nonconforming interfaces naturally. With a linear approximation of the temperature on the interfaces, the maximum local error is near 11.4%. As expected, if approximation degree on the interfaces is increased, the error decreases and its maximum drops below 5.6%. Note also that, as commented previously, only three Newton-Raphson iterations suffice to converge with a tolerance of 10 -8 and this is independent on the number of interfaces.

Finally, table 3 shows the computational time of the offline and online phases and the number of modes considered for the different subdomain types, with a mesh of 22 × 110 elements per subdomain. As in the previous examples, when the degree at the interfaces increases, the online computational cost slightly increases. The FEM solution for this problem requires a mesh of 110 × 110 elements for each of the 16 domains, with a total amount of 193600 elements, with a computational cost of 132 .

APPENDIX A: DIRICHLET CONDITIONS AS ESSENTIAL BOUNDARY CONDITIONS

The local problem in each subdomain type can also be solved imposing Dirichlet boundary conditions as essential ones. In this case, the weak problem associated to [START_REF] Chinesta | PGD for solving multidimensional and parametric models[END_REF]

with

where

and  1 is the standard Sobolev space of square integrable functions and first derivatives.

This approach, however, is a little more cumbersome for the PGD technique since it requires the evaluation of the first mode in a more specific manner. It is standard in PGD [START_REF] Chinesta | PGD for solving multidimensional and parametric models[END_REF][START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations[END_REF] to impose the essential Dirichlet boundary conditions in the first mode and then solve the other modes with homogeneous Dirichlet boundary conditions.

Here, however, the Dirichlet boundary conditions on the interfaces are characterized as a set of parameters that PGD sees as new coordinates. This implies also a specific strategy for the first mode whereas all the others are, as usual, solved with homogeneous conditions.

In fact the strategy only affects the first spatial modes. In fact, the same separated representation proposed in ( 9) is used. Two extra terms are added, one to take care of the imposed Dirichlet conditions in Γ and another to account for the parameters û. The former is standard and requires solving the separated problem with homogeneous conditions everywhere except for Ω ∩ Γ to compute 1 ( 1 ) and 2 ( 2 ). The second term is the one associated to each parameter û for = 1, … , ̂ , see [START_REF] Ghnatios | 3D Modeling of squeeze flows occurring in composite laminates[END_REF]. It consists in solving for the spatial modes 1 ( 1 ) and 2 ( 2 ) as before but with homogeneous conditions everywhere except for Ω ∩ Γ ∩ supp( ̂ ) where the condition = ̂ is imposed.

Recalling the expression (9) the PGD approximation for each subdomain introduces these new terms and becomes

This strategy, requires to interpolate the shape functions ̂ for the interface variables ̂ on the finite element mesh used in the subdomain and special care must be taken of the corners of the subdomain in order to avoid duplicity in the imposed boundary condition.