# Arlequin based PGD domain decomposition

S. Mohamed Nazeer, Adrien Leygue, Francisco Chinesta, Felipe Bordeu

► To cite this version:

S. Mohamed Nazeer, Adrien Leygue, Francisco Chinesta, Felipe Bordeu. Arlequin based PGD domain decomposition. Computational Mechanics, 2014, 54 (5), pp.1175-1190. 10.1007/s00466-014-1048-7 . hal-04325797

HAL Id: hal-04325797

https://hal.science/hal-04325797v1

Submitted on 11 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Arlequin based PGD domain decomposition

S. Mohamed Nazeer, Felipe Bordeu, Adrien Leygue, Francisco Chinesta

**Abstract** Problems defined in fully or partially separable domains can be solved by considering a space separated representation of the unknown fields. Thus three-dimensional problems can be solved from the solution of some one-dimensional problems in the case of fully separated representations involving the three space coordinates or as a sequence of 2D and 1D problems in the case of partially separated representations (plates, shells or extruded geometries). When the domains become more complex, sometimes they can be simplified by using appropriate mappings. When it is not possible or such a transformation becomes too complex, the use of domain decomposition could facilitate the use of separated representations. However, domain coupling in the context of space separated representations have never been analyzed. In this paper we propose a domain decomposition strategy based on the use of space separated representations and the Arlequin coupling strategy. First we consider separated representations of the physical space that will be then extended to address parametric solutions.

**Keywords** Separated representations · Domain decomposition · Proper generalized decomposition · Arlequin · Parametric analysis

## 1 Introduction

Problems defined in spaces of moderate dimension are encountered in many branches of engineering sciences (e.g. computational solid mechanics, computational fluid mechanics, thermal engineering, computational electromagnetism, etc.). Different discretization techniques have been proposed, analyzed and validated for addressing the solution of such models, even when they involves strong non linearities and multi-physics couplings (e.g. finite elements, finite volumes, finite differences, meshless techniques, etc.).

When the complexity of such models increases because of the geometrical complexity of the domains in which they are defined or because of their multi-scale character, efficient solutions require the use of advanced techniques and powerful computers. Multi-scale approaches, domain decomposition, techniques for accelerating the convergence of iterative solvers, multigrid algorithms, parallel computing, among many others, make possible enlarging more and more the frontiers of numerical simulation allowing the use of high fidelity models in most branches of engineering science.

However, despite the impressive progresses attained during the last decades the efficient solution of many problems remains a challenge. One of these challenging scenarios concerns the solution of models defined in high-dimensional spaces as the ones encountered in the fine description of mechanics of materials (quantum chemistry) or the kinetic theory description of nanostructured fluids where models are defined in spaces involving *D* dimensions with *D* ≫ 1. In this case the use of mesh-based discretizations fail because the complexity of an hypothetical mesh scales exponentially with the space dimension. In one dimension functions can be approximated from their values in a discrete number of points (nodes) *M*, from which they can be interpolated everywhere in the domain. In two dimensions when using *M* nodes in

<footnote>S. M. Nazeer · F. Bordeu · A. Leygue · F. Chinesta (✉)
GeM, UMR CNRS - Ecole Centrale de Nantes, 1 Rue de la Noe,
BP 92101, 44321 Nantes cedex 3, France
e-mail: francisco.chinesta@ec-nantes.fr</footnote>

<footnote>S. M. Nazeer
e-mail: mohamed.shahul-hameed@ec-nantes.fr</footnote>

<footnote>F. Bordeu
e-mail: felipe.bordeu@ec-nantes.fr</footnote>

<footnote>A. Leygue
e-mail: adrien.leygue@ec-nantes.fr</footnote>
each direction, it results $M^2$ nodes, and in $D$ dimensions the number of points in which the solution must be evaluated results $M^D$. When $D \gg 3$ mesh-based solutions are unaffordable despite the nowadays computational availabilities.

In those scenarios the use of separated representations allows sometimes circumventing the curse of dimensionality. Hartree-Fock and post-Hartree-Fock approaches considered such representations [12]. Later in [1] and [2] separated representations were considered for describing molecular systems involving numerous conformational dimensions in the framework of the so-called Proper Generalized Decomposition – PGD –. If we consider a field defined in a space of dimension $D$ the simplest separated approximation of it reads

$$u(x_1, \dots, x_D) \approx \sum_{i=1}^{N} \prod_{j=1}^{D} X_i^j(x_j) \quad (1)$$

Thus, the solution procedure consists of calculating the one-dimensional functions $X_i^j(x_j)$.

Even if the use of this kind of representations seems fully justified in the case of high dimensional models ($D > 3$) because in that case there is not valuable alternatives, separated representations were also considered for solving standard 3D models by separating the space coordinates. When the problem is defined in a fully separable domain $\Omega$, i.e. $\Omega = \Omega_x \times \Omega_y \times \Omega_z$ the most appealing decomposition consists of a fully space separated representation

$$u(x, y, z) = \sum_{i=1}^{N} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \quad (2)$$

Thus the 3D solution is performed from a sequence of one-dimensional solutions for calculating functions $X_i(x)$, $Y_i(y)$ and $Z_i(z)$, $\forall i = 1, \dots, N$.

One can also proceed by separating the coordinates into several groups. For example, consider a three-dimensional extruded domain $\Omega = \Omega_{xy} \times \Omega_z$. Here, the domain $\Omega$ is the extrusion of the non-separable two-dimensional domain $\Omega_{xy}$ along the $z$ axis. In this case, the most natural decomposition reads

$$u(x, y, z) = \sum_{i=1}^{N} X_i(x, y) \cdot Z_i(z). \quad (3)$$

This approach is particularly well suited to problems involving plate and shell geometries. The calculations thus involve a series of decoupled two-dimensional problems in $\Omega_{xy}$ to compute functions $X_i(x, y)$, and one-dimensional problems in $\Omega_z$ to compute functions $Z_i(z)$. Thus the fully three-dimensional simulation has a numerical complexity typical of two-dimensional analyses.

This approach was considered in our former works in the analysis of extruded profiles, plates, shells and composites laminates [9, 10, 13, 26]. The computing time savings can be enormous when using such a separated representation. If we consider the 3D solution in a plate domain approximated from $M_x \times M_y \times M_z$ nodes, $M_z$ being the nodes distributed in the domain thickness, its calculation requires the solution of a linear system of size $M_x \times M_y \times M_z$, however when considering the separated representation (3) we must solve of the order of $N$ 2D problems of size $M_x \times M_y$ for determining functions $X_i(x, y)$ and $N$ one-dimensional problems of size $M_z$ for calculating functions $Z_i(z)$. As proved in [9] when $M_x \times M_y$ or $M_z$ increase the use of separated representations allows computing time savings of many orders of magnitude (see Fig. 4 in reference [9]).

As illustrated in [9–11] separated representations can be considered in more complex geometries as soon as the material parameters can be expressed in a separated form. For example when considering an elastic problem in the geometry depicted in Fig. 1, the solution (displacement field) can be expressed in a separated form if we consider the following separated form of the Young’s modulus

$$E(x, y, z) \approx E \cdot \sum_{i=1}^{2} e_i(x, y) \cdot \epsilon_i(z) \quad (4)$$

where the characteristic functions $e_i(x, y)$, $\epsilon_i(z)$ ($i = 1, 2$) are depicted in Fig. 2 [10].

Despite the interest of performing such a separated representation for alleviating the computational cost of 3D

<figure><img src="image_2.png" /><figcaption>Fig. 1 Complex domain</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 Material domain characteristic functions defining the Young’s modulus</figcaption></figure>
calculation, such decomposition fails when the geometry does not fit the coordinate axes, as is the case of the geometry depicted in Fig. 3. Thus we can conclude that the main difficulty of separated representations involving the space coordinates concerns the treatment of substructures assemblages, even if for each one of them a fully or partial separated representation works (e.g. plates or shells assembling). In this case domain decomposition strategies seem a valuable alternative for addressing complex structures in which separated representations are considered at the level of each substructure.

Domain decomposition techniques have been widely employed in a diversity of applications. Some techniques are based on the subdomains overlapping, other consider non-overlapping coupling strategies. Among the overlapping-based strategies we can mention the classical Schwartz [28] or the more recent and versatile Arlequin technique and its variants [5–7] (the interested reader can refer to the numerous works of Ben Dhia, the initiator of this powerful numerical technique). Among the non-overlapping ones we can mention the mortar-based techniques [8,27], the FETI strategy and its variants [16,17] and the LATIN-multiscale and LATIN-multidomain approaches widely developed by Ladeveze and coauthors [20–24]. These works considered multi domains coupled with space-time separated representations, however in our knowledge, never separated representations involving either space coordinates or parametric extra-coordinates were considered within the framework of domain decomposition.

In this work we are analyzing the use of domain decomposition combined with separated representations of the problem solution. These separated representations will involve both the space coordinates and also model parameters that within the PGD framework are considered as extra-coordinates. The Arlequin technique will be considered for performing the subdomains coupling. For this purpose, we summarize in Sect. 2 the space separated representation. In Sect. 3 we revisits the Arlequin-based coupling strategy on a simple model defined in a geometry representative of an structural assembling. This strategy will be extended to separated representations in Sect. 4. The numerical results described in Sect. 5 allow to conclude about the feasibility and potentiality of such a coupling methodology.

This work aims only proving the feasibility of considering space separated representations in complex structures when considering the decomposition at the level of each substructure and then coupling all them by using an appropriate domain decomposition technique. The numerical examples here addressed are enough for proving its feasibility, but they are not complex enough for quantifying the CPU time savings, that are expected to be significant in 3D as just mentioned. This analysis constitutes a work in progress.

## 2 Space separated representation and proper generalized decomposition – PGD –

Consider the solution of the Poisson equation

$$ \Delta u(x, y) = f(x, y), \qquad (5) $$

in a two-dimensional rectangular domain $\Omega = \Omega_x \times \Omega_y = (0, L) \times (0, H)$.

We specify homogeneous Dirichlet boundary conditions for the unknown field $u(x, y)$, i.e. $u(x, y)$ vanishes at the domain boundary $\Gamma = \partial \Omega$. Furthermore, we assume, without loss of generality, that the source term $f$ is constant over the domain $\Omega$.

For all suitable test functions $u^*$, the weighted residual form of (5) reads

$$ \int_{\Omega_x \times \Omega_y} u^* (\Delta u - f) \, dx \, dy = 0, \qquad (6) $$

or more explicitly

$$ \int_{\Omega_x \times \Omega_y} u^* \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} - f \right) \, dx \, dy = 0. \qquad (7) $$

Our goal is to obtain an approximate solution to (5) in the separated form

$$ u(x, y) = \sum_{i=1}^{N} X_i(x) \cdot Y_i(y). \qquad (8) $$

We shall do so by computing each term of the expansion one at a time, thus enriching the approximation until a suitable convergence criterion is satisfied.

At each enrichment step $n$ ($n \ge 1$), we have already computed the $n-1$ first terms of the approximation (8):

$$ u^{n-1}(x, y) = \sum_{i=1}^{n-1} X_i(x) \cdot Y_i(y). \qquad (9) $$

We now wish to compute the next term $X_n(x) \cdot Y_n(y)$ to obtain the enriched PGD solution

$$ \begin{aligned} u^n(x, y) &= u^{n-1}(x, y) + X_n(x) \cdot Y_n(y) \\ &= \sum_{i=1}^{n-1} X_i(x) \cdot Y_i(y) + X_n(x) \cdot Y_n(y) \end{aligned} \qquad (10) $$

Both functions $X_n(x)$ and $Y_n(y)$ are unknown at the current enrichment step $n$, and they appear in the form of a
product. The resulting problem is thus non-linear and a suitable iterative scheme is required. We shall use the index p to denote a particular iteration.

At enrichment step n, the PGD approximation $u^{n,p}$ obtained at iteration p thus reads

$$u^{n,p}(x, y) = u^{n-1}(x, y) + X_n^p(x) \cdot Y_n^p(y). \quad (11)$$

The simplest iterative scheme is an alternating direction strategy that computes $X_n^p(x)$ from $Y_n^{p-1}(y)$, and then $Y_n^p(x)$ from $X_n^p(x)$. An arbitrary initial guess $Y_n^0(y)$ is specified to start the iterative process. The non-linear iterations proceed until reaching a fixed point within a user-specified tolerance $\epsilon$, i.e.

$$\frac{\|X_n^p(x) \cdot Y_n^p(y) - X_n^{p-1}(x) \cdot Y_n^{p-1}(y)\|}{\|X_n^{p-1}(x) \cdot Y_n^{p-1}(y)\|} < \epsilon, \quad (12)$$

where $\| \cdot \|$ is a suitable norm.

The enrichment step n thus ends with the assignments $X_n^p(x) \to X_n(x)$ and $Y_n^p(y) \to Y_n(y)$.

The enrichment process itself stops when an appropriate measure of error $\mathcal{E}(n)$ becomes small enough, i.e $\mathcal{E}(n) < \tilde{\epsilon}$. Several stopping criteria are available [3, 15, 25]. Appendix 7 describes in more detail one particular alternating direction iteration at a given enrichment step.

### 3 Arlequin-based coupling strategy

In this section instead of considering the simple domain just addressed that can be fully decomposed as a cartesian product of one-dimensional domains, i.e. $\Omega = \Omega_x \times \Omega_y$ we are considering a more complex one that does not accept a simple direct decomposition, as the one depicted in Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Complex domain</figcaption></figure>

In [19] a technique for addressing complex geometries was proposed. It proceeds by embedding the real domain $\Omega$ into a parallelepiped one $\Theta$ fully separable, $\overline{\Omega} \subset \Theta$, and then applying an adequate penalty formulation in the complementary domain $\Theta \setminus \Omega$. There are different alternatives for enforcing boundary conditions as discussed in [19] or [18].

In this paper we are exploring more natural alternatives, as the one related to domain decomposition, because it could allow treating easily more complex geometries with the possibility of proceeding by using parallel computing platforms.

Obviously one could expect that non-overlapping domain decomposition techniques could be simpler than implement for addressing the geometry depicted in Fig. 3, however it is not the case. The simplest non-overlapping decomposition consists of two rectangles each one mapping one of the main quadrilaterals in domain $\Omega$. But after this simple geometrical transformation one must enforce the continuity conditions between the solutions living at each subdomain. Here we find two main difficulties: (i) the meshes are not necessarily compatible along the two domains boundaries along which the two domains are glued; and (ii) even in the case of using compatible meshes the separated representation implies a difficulty to express the solution continuity. There is not a direct way of enforcing the solution continuity (only the final solution must be continuous but not necessarily each mode involved in the separated representations). Thus, the use of an overlapping strategy seems as simple as the use of a non-overlapping one and moreover the first can be extended easily for addressing multi-models. For this reason we considered the Arlequin strategy, described in what follows, and even if the numerical examples here solved could seem quite simple, it is not the case when using separated representations. As soon as we define and validate the formulation on the overlapping patch, complex structures can be assembled as a sort of “lego”.

Within the Arlequin framework, the one considered in this work, the domain $\Omega$ depicted in Fig. 3 is decomposed in two domains $\Omega_1$ and $\Omega_2$ that overlap in the region $\Omega_o$ as depicted in Fig. 4.

The model equations must be enforced in both subdomains $\Omega_1$ and $\Omega_2$ with an appropriate coupling in the overlapping domain $\Omega_o$ that consists of enforcing the model equations while ensuring the equality of the solutions related to each subdomain $u_1(X, Y)$ and $u_2(X, Y)$. This equality will be enforced by introducing a Lagrange multiplier and considering a L2 or a H1 norm. In what follows we consider the L2 norm. At the same time the model enforcement in the overlapping domain $\Omega_o$ requires an adequate partition between the equations enforced in each domain that overlaps in $\Omega_o$ by using an appropriate weight function that can be, when applied on the weak form, continuous or discontinuous.
<figure><img src="image_5.png" /><figcaption>Fig. 4 Domain decomposition</figcaption></figure>

To analyze the coupling procedure we consider the problem

$$ \Delta u(X, Y) + 1 = 0 \qquad (13) $$

defined in the domain $\Omega$ depicted in Fig. 3. For the sake of simplicity and without loss of generality we consider homogenous Dirichlet boundary conditions on the whole boundary of $\Omega$, $\partial\Omega$. The unknown field in each subdomain will be noted by $u_1(\mathbf{X} \in \Omega_1)$ and $u_2(\mathbf{X} \in \Omega_2)$.

The weak form of Eq. (13) reads: Find $u(\mathbf{X}) \in \mathcal{H}_0^1(\Omega)$ such that

$$ \int_{\Omega} \nabla u^* \cdot \nabla u \, d\mathbf{X} = \int_{\Omega} u^* \, d\mathbf{X}, \quad \forall u^* \in \mathcal{H}_0^1(\Omega) \qquad (14) $$

where $\mathcal{H}_0^1(\Omega) = \{u(\mathbf{X}) \in \mathcal{H}^1(\Omega), u(\mathbf{X} \in \partial\Omega) = 0\}$.

When decomposing the domain into $\Omega_1$ and $\Omega_2$ the Arlequin method considers the following weak form based

on the use of a Lagrange multiplayer $\lambda(\mathbf{X} \in \Omega_o)$ for enforcing the equality of $u_1(\mathbf{X})$ and $u_2(\mathbf{X})$ in $\Omega_o$:

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega_1} \xi_1 (\nabla u_1^* \cdot \nabla u_1) \, d\mathbf{X} + \int_{\Omega_o} u_1^* \lambda \, d\mathbf{X} = \int_{\Omega_1} \xi_1 u_1^* \, d\mathbf{X} \\[1.5em] \displaystyle \int_{\Omega_2} \xi_2 (\nabla u_2^* \cdot \nabla u_2) \, d\mathbf{X} - \int_{\Omega_o} u_2^* \lambda \, d\mathbf{X} = \int_{\Omega_2} \xi_2 u_2^* \, d\mathbf{X} \\[1.5em] \displaystyle \int_{\Omega_o} \lambda^*(u_1 - u_2) \, d\mathbf{X} = 0 \end{array} \right. \qquad (15) $$

where $\lambda \in \mathcal{L}^2(\Omega_o)$, $u_1 \in \mathcal{H}^1(\Omega_1)$ with $u_1(\mathbf{X} \in \partial\Omega_1 \cap \partial\Omega) = 0$ and $u_2 \in \mathcal{H}^1(\Omega_2)$ with $u_2(\mathbf{X} \in \partial\Omega_2 \cap \partial\Omega) = 0$.

The weight functions $\xi_1(\mathbf{X} \in \Omega_1)$ and $\xi_2(\mathbf{X} \in \Omega_2)$ that can be continuous or not, must define a partition of unity on $\Omega_o$ where both energy contributions are weighted:

$$ \left\{ \begin{array}{ll} \xi_1(\mathbf{X}) = 1 & \text{if } \mathbf{X} \in \Omega_1 \setminus \Omega_o \\ \xi_2(\mathbf{X}) = 1 & \text{if } \mathbf{X} \in \Omega_2 \setminus \Omega_o \\ \xi_1(\mathbf{X}) + \xi_2(\mathbf{X}) = 1 & \text{if } \mathbf{X} \in \Omega_o \end{array} \right. \qquad (16) $$

## 4 Separated representation

The decomposition illustrated in Fig. 4 does not allow a direct separated representation involving the space coordinates $(X, Y)$ that would require regular rectangular domains with their faces parallel to the coordinates axes.

### 4.1 Domain mapping

We consider domains $\Omega_1$ and $\Omega_2$ decomposed respectively in the 6 and 2 quadrilaterals depicted in Fig. 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5 Domain mapping</figcaption></figure>
We denote these quadrilaterals by $Q_1^j, j = 1, \dots, 6$ for the ones covering $\Omega_1$ and by $Q_2^k, k = 1, 2$ for the ones covering $\Omega_2$. Two of them, one related to each subdomain, correspond to the overlapping domain $\Omega_o$.

The integral form on $\Omega_1$ in Eq. (15) can be rewritten as:

$$ \sum_{j=1}^{6} \int_{Q_1^j} \xi_1 (\nabla u_1^* \cdot \nabla u_1) \, d\mathbf{X} + \int_{\Omega_o} u_1^* \lambda \, d\mathbf{X} = \sum_{j=1}^{6} \int_{Q_1^j} \xi_1 u_1^* \, d\mathbf{X} \quad (17) $$

and the one on $\Omega_2$ as

$$ \sum_{k=1}^{k=2} \int_{Q_2^k} \xi_2 (\nabla u_2^* \cdot \nabla u_2) \, d\mathbf{X} - \int_{\Omega_o} u_2^* \lambda \, d\mathbf{X} = \sum_{k=1}^{k=2} \int_{Q_2^k} \xi_2 u_2^* \, d\mathbf{X} \quad (18) $$

Now, we can define the mapping between the real domain and the one depicted in Fig. 5 (right) in which the transformed domains of $\Omega_1, \Omega_2$ and $\Omega_o$, denoted by $\mathcal{T}_1, \mathcal{T}_2$ and $\mathcal{O}_o$ respectively, allow a fully space separation. Quadrilaterals $Q_1^j$ and $Q_2^k$ are mapped onto $\mathcal{V}_1^j$ and $\mathcal{V}_2^k$ respectively, all them rectangular and parallel to the coordinates axes $(x, y)$ (see Fig. 5 right).

Mappings $Q_1^j \to \mathcal{V}_1^j$ and $Q_2^k \to \mathcal{V}_2^k$, and their inverse transformations, can be easily obtained as well their associated Jacobian $\mathcal{J}_1^j$ and $\mathcal{J}_2^k$. We denote by $\mathcal{J}_o$ the Jacobian related to the overlapping domain mapping $\Omega_o \to \mathcal{O}_o$.

Using all this notation integrals forms (15) become:

$$ \sum_{j=1}^{6} \int_{\mathcal{V}_1^j} \xi_1 \nabla u_1^{*T} (\mathcal{J}_1^j)^{-T} (\mathcal{J}_1^j)^{-1} \nabla u_1 \det (\mathcal{J}_1^j) \, d\mathbf{x} + \int_{\mathcal{O}_o} u_1^* \lambda \det(\mathcal{J}_o) \, d\mathbf{x} = \sum_{j=1}^{6} \int_{\mathcal{V}_1^j} \xi_1 u_1^* \det (\mathcal{J}_1^j) \, d\mathbf{x} \quad (19) $$

$$ \sum_{k=1}^{2} \int_{\mathcal{V}_2^k} \xi_2 \nabla u_2^{*T} (\mathcal{J}_2^k)^{-T} (\mathcal{J}_2^k)^{-1} \nabla u_2 \det (\mathcal{J}_2^k) \, d\mathbf{x} - \int_{\mathcal{O}_o} u_2^* \lambda \det(\mathcal{J}_o) \, d\mathbf{x} = \sum_{k=1}^{2} \int_{\mathcal{V}_2^k} \xi_2 u_2^* \det (\mathcal{J}_2^k) \, d\mathbf{x} \quad (20) $$

and

$$ \int_{\mathcal{O}_o} \lambda^*(u_1 - u_2) \det(\mathcal{J}_o) \, d\mathbf{x} = 0 \quad (21) $$

where for the sake of notational simplicity we omitted the dependence of all the fields $\lambda, u_1$ and $u_2$ on the new coordinates $(x, y)$ as well as the definition of the gradient operator $\nabla$ in the new coordinates, that is $\nabla = (\partial/\partial x, \partial/\partial y)$.

4.2 Space separated representation

The separated representation of the different unknown fields is defined according to:

$$ \left\{ \begin{aligned} u_1(x, y) &\approx \sum_{i=1}^{N} X_1^i(x) \cdot Y_1^i(y) \\ u_2(x, y) &\approx \sum_{i=1}^{N} X_2^i(x) \cdot Y_2^i(y) \\ \lambda(x, y) &\approx \sum_{i=1}^{N} L_x^i(x) \cdot L_y^i(y) \end{aligned} \right. \quad (22) $$

This separated representation is constructed progressively, one term of the sum at each step. Thus if we assume the solution given at step $n-1$ by

$$ \left\{ \begin{aligned} u_1^{n-1}(x, y) &= \sum_{i=1}^{n-1} X_1^i(x) \cdot Y_1^i(y) \\ u_2^{n-1}(x, y) &= \sum_{i=1}^{n-1} X_2^i(x) \cdot Y_2^i(y) \\ \lambda^{n-1}(x, y) &= \sum_{i=1}^{n-1} L_x^i(x) \cdot L_y^i(y) \end{aligned} \right. \quad (23) $$

at step $n$ the next term of the sums will be calculated according to

$$ \left\{ \begin{aligned} u_1^n(x, y) &= \sum_{i=1}^{n-1} X_1^i(x) \cdot Y_1^i(y) + X_1^n(x) \cdot Y_1^n(y) \\ &= u_1^{n-1} + X_1^n(x) \cdot Y_1^n(y) \\ u_2^n(x, y) &= \sum_{i=1}^{n-1} X_2^i(x) \cdot Y_2^i(y) + X_2^n(x) \cdot Y_2^n(y) \\ &= u_2^{n-1} + X_2^n(x) \cdot Y_2^n(y) \\ \lambda^n(x, y) &= \sum_{i=1}^{n-1} L_x^i(x) \cdot L_y^i(y) + L_x^n(x) \cdot L_y^n(y) \\ &= \lambda^{n-1} + L_x^n(x) \cdot L_y^n(y) \end{aligned} \right. \quad (24) $$

Test functions within a Galerkin setting write:

$$ \left\{ \begin{aligned} u_1^*(x, y) &= X_1^*(x) \cdot Y_1^n(y) + X_1^n(x) \cdot Y_1^*(y) \\ u_2^*(x, y) &= X_2^*(x) \cdot Y_2^n(y) + X_2^n(x) \cdot Y_2^*(y) \\ \lambda^*(x, y) &= L_x^*(x) \cdot L_y^n(y) + L_x^n(x) \cdot L_y^*(y) \end{aligned} \right. \quad (25) $$

The transformed subdomains $\mathcal{T}_1, \mathcal{T}_2$ and $\mathcal{O}_o$ being rectangular and with their faces parallel to the coordinate axes
(x, y) they can be decomposed according to:

$$ \left\{ \begin{array}{l} \mathcal{T}_1 = \mathcal{T}_1^x \times \mathcal{T}_1^y \\ \mathcal{T}_2 = \mathcal{T}_2^x \times \mathcal{T}_2^y \\ \mathcal{O}_o = \mathcal{O}_o^x \times \mathcal{O}_o^y \end{array} \right. \qquad (26) $$

Functions $X_1^i(x)$ can be approximated in $\mathcal{T}_1^x$ by using any standard one dimensional approximation technique. The simplest one consists of using piecewise linear interpolants, like the ones used in finite elements. Thus if $M_1^x$ nodes located at positions $x_r$, $r = 1, \dots, M_1^x$, are considered on $\mathcal{T}_1^x$ the interpolation writes

$$ X_1^i(x) = \sum_{r=1}^{M_1^x} X_1^i(x_r) N_r^{x,1}(x) \qquad (27) $$

where $N_r^{x,1}(x)$ are the usual shape functions and where the superscript “1” refers to the considered domain.

If $M_1^y$ are the number of nodes located at positions $y_s \in \mathcal{T}_1^y$ used for defining the approximation of functions $Y_1^i(y)$ it results

$$ Y_1^i(y) = \sum_{s=1}^{M_1^y} Y_1^i(y_s) N_s^{y,1}(y) \qquad (28) $$

The approximation in the second transformed domain $\mathcal{T}_2$ defined from $M_2^x$ and $M_2^y$ nodes located on $\mathcal{T}_2^x$ and $\mathcal{T}_2^y$ respectively, reads:

$$ X_2^i(x) = \sum_{r=1}^{M_2^x} X_2^i(x_r) N_r^{x,2}(x) \qquad (29) $$

and

$$ Y_2^i(y) = \sum_{s=1}^{M_2^y} Y_2^i(y_s) N_s^{y,2}(y) \qquad (30) $$

Finally for approximating the Lagrange multiplier we consider $M_\lambda^x$ and $M_\lambda^y$ nodes on $\mathcal{O}_o^x$ and $\mathcal{O}_o^y$ respectively, allowing writing:

$$ L_x^i(x) = \sum_{r=1}^{M_\lambda^x} L_x^i(x_r) \bar{N}_r^{x,\lambda}(x) \qquad (31) $$

and

$$ L_y^i(y) = \sum_{s=1}^{M_\lambda^y} L_y^i(y_s) \bar{N}_s^{y,\lambda}(y) \qquad (32) $$

where the shape functions used for approximating the functions involved in the Lagrange multiplier approximation can be different to the ones used for approximating the functions involved in the expression of $u_1$ and $u_2$. In our numerical implementation we considered both, piecewise linear functions and piecewise constant functions. When two different mesh sizes are considered in $\mathcal{T}_1$ and $\mathcal{T}_2$, the one considered in the coupling domain $\mathcal{O}_o$ was the coarser one.

The efficient application of the separated representation constructor, the one allowing the calculation of the different unknown fields (here $u_1$, $u_2$ and $\lambda$) requires that all the functions involved in the weak form have a separated representation. In our case we must express functions $\xi_1(x, y)$, $\xi_2(x, y)$, $\mathcal{J}_1^j$, $\mathcal{J}_2^k$ as well as their inverse transformations $(\mathcal{J}_1^j)^{-1}$ and $(\mathcal{J}_2^k)^{-1}$ in a separated form.

There is a general tool for performing such separated representation. It relies in using a singular value decomposition –SVD– (in 2D) or its multidimensional counterpart, the so-called high order singular value decomposition – HOSVD – (when $D > 2$). However, in the case here addressed a separated representation can be directly obtained (see appendix 8)

Finally we must discuss the choice of function $\xi_1(\mathbf{x} \in \mathcal{O}_o)$ (function $\xi_2(\mathbf{x})$ being $1 - \xi_1(\mathbf{x})$ inside the overlapping domain $\mathcal{O}_o$).

A first possibility consists of defining a regular enough function ensuring the continuity with $\xi_1(\mathbf{x} \in \Omega_1 \setminus \Omega_o)$. For this purpose it suffices solving

$$ \Delta \xi(\mathbf{x}) = 0, \quad \mathbf{x} \in \mathcal{O}_o \qquad (33) $$

with

$$ \left\{ \begin{array}{l} \xi(\mathbf{x} \in \partial \mathcal{O}_o \cap \mathcal{T}_1) = 1 \\ \xi(\mathbf{x} \in \partial \mathcal{O}_o \cap \partial \mathcal{T}_1) = 0 \end{array} \right. \qquad (34) $$

The solution of this problem is searched in the separated form

$$ \xi(x, y) \approx \sum_{i=1}^{P} F_i(x) \cdot G_i(y) \qquad (35) $$

by using the standard PGD algorithm. As soon as this solution is available we define $\xi_1(\mathbf{x}) = \xi(\mathbf{x})$ and $\xi_2(\mathbf{x}) = 1 - \xi_1(\mathbf{x})$.

In this work we compared such continuous choice of function $\xi_1(\mathbf{x})$ and $\xi_2(\mathbf{x})$ with another less regular but computationally cheaper defined by

$$ \xi_1(x, y) = 1 - \frac{x - x_o^-}{x_o^+ - x_o^-}, \quad x \in [x_o^-, x_o^+] \qquad (36) $$

where $x_o^-$ and $x_o^+$ define the x-interval $\mathcal{O}_o^x$.
## 4.3 Separated representation constructor

With all the functions and fields involved in the weak form expressed in a separated form, we can proceed to apply the PGD technique in order to calculate the different functions involved in the unknown fields representation.

At iteration $n$, $u_1^{n-1}(\mathbf{x})$, $u_2^{n-1}(\mathbf{x})$ and $\lambda^{n-1}(\mathbf{x})$ are assumed known and we look for the new term according to Eq. (24). Because the problem remains nonlinear, we consider the simplest linearization whose robustness was proved in our former works. Thus assuming functions depending on the y-coordinate known we will compute those depending on the x-coordinate, and from the just computed functions we will update the ones depending on y. The iteration process continues until reaching the convergence (fixed point).

### 4.3.1 Computing functions depending on the x-coordinate

By introducing into the weak form (15) the separated representations of all the unknown fields and functions, as well as the expression of the test functions

$$ \left\{ \begin{array}{l} u_1^*(x, y) = X_1^*(x) \cdot Y_1^n(y) \\ u_2^*(x, y) = X_2^*(x) \cdot Y_2^n(y) \\ \lambda^*(x, y) = L_x^*(x) \cdot L_y^n(y) \end{array} \right. \qquad (37) $$

and performing all the integrals involving the y-coordinate (all the functions depending on the y-coordinate being known), it results

$$ \int_{\mathcal{T}_1^x} \left( \alpha_1^y(x) \frac{dX_1^*}{dx} \frac{dX_1^n}{dx} + \beta_1^y(x) X_1^* X_1^n \right) dx \\ + \int_{\mathcal{O}_o^x} \gamma_1^y(x) X_1^* L_x^n dx = \int_{\mathcal{T}_1^x} X_1^* f_1^y(x) dx \\ + \int_{\mathcal{T}_1^x} \frac{dX_1^*}{dx} g_1^y(x) dx \qquad (38) $$

where functions $\alpha_1^y(x)$, $\beta_1^y(x)$, $\gamma_1^y(x)$, $f_1^y(x)$ and $g_1^y(x)$ come from all the integrals in $\mathcal{T}_1^y$. Analogously,

$$ \int_{\mathcal{T}_2^x} \left( \alpha_2^y(x) \frac{dX_2^*}{dx} \frac{dX_2^n}{dx} + \beta_2^y(x) X_2^* X_2^n \right) dx \\ - \int_{\mathcal{O}_o^x} \gamma_2^y(x) X_2^* L_x^n dx = \int_{\mathcal{T}_2^x} X_2^* f_2^y(x) dx \\ + \int_{\mathcal{T}_2^x} \frac{dX_2^*}{dx} g_2^y(x) dx \qquad (39) $$

and

$$ \int_{\mathcal{O}_o^x} \delta^y(x) L_x^*(X_1^n - X_2^n) dx = \int_{\mathcal{O}_o^x} L_x^* h^y(x) dx \qquad (40) $$

By using a one-dimensional finite element approximation the problem above can be solved for obtaining the searched fields $X_1^n(x)$, $X_2^n(x)$ and $L_x^n(x)$.

### 4.3.2 Computing functions depending on the y-coordinate

By introducing now into the weak form (15) the separated representations of all the unknown fields and functions, as well as the expression of the test functions that now read

$$ \left\{ \begin{array}{l} u_1^*(x, y) = X_1^n(x) \cdot Y_1^*(y) \\ u_2^*(x, y) = X_2^n(x) \cdot Y_2^*(y) \\ \lambda^*(x, y) = L_x^n(x) \cdot L_y^*(y) \end{array} \right. \qquad (41) $$

and performing all the integrals involving the x-coordinate (all the functions depending on the x-coordinate being known), it results

$$ \int_{\mathcal{T}_1^y} \left( \alpha_1^x(y) \frac{dY_1^*}{dy} \frac{dY_1^n}{dy} + \beta_1^x(y) Y_1^* Y_1^n \right) dy \\ + \int_{\mathcal{O}_o^y} \gamma_1^x(y) Y_1^* L_y^n dy = \int_{\mathcal{T}_1^y} Y_1^* f_1^x(y) dy \\ + \int_{\mathcal{T}_1^y} \frac{dY_1^*}{dy} g_1^x(y) dy \qquad (42) $$

where functions $\alpha_1^x(y)$, $\beta_1^x(y)$, $\gamma_1^x(y)$, $f_1^x(y)$ and $g_1^x(y)$ come from all the integrals on $\mathcal{T}_1^x$. Analogously,

$$ \int_{\mathcal{T}_2^y} \left( \alpha_2^x(y) \frac{dY_2^*}{dy} \frac{dY_2^n}{dy} + \beta_2^x(y) Y_2^* Y_2^n \right) dy \\ - \int_{\mathcal{O}_o^y} \gamma_2^x(y) Y_2^* L_y^n dy = \int_{\mathcal{T}_2^y} Y_2^* f_2^x(y) dy \\ + \int_{\mathcal{T}_2^y} \frac{dY_2^*}{dy} g_2^x(y) dy \qquad (43) $$

and

$$ \int_{\mathcal{O}_o^y} \delta^x(y) L_y^*(Y_1^n - Y_2^n) dy = \int_{\mathcal{O}_o^y} L_y^* h^x(y) dy \qquad (44) $$
By using a one-dimensional finite element approximation the problem above can be solved for obtaining the searched fields $Y_1^n(y)$, $Y_2^n(y)$ and $L_y^n(y)$.

When considering the separated representation constructor just described, the usual way to enforce Dirichlet boundary conditions consists of choosing some couples of separated functions regular enough and that verifies the boundary conditions, and add those terms to the first terms of the separated representation of $u(x, y)$. This procedure works very well and it is extremely easy to implement when we consider a single domain as was the case in our former works. Here the situation is a bit different because we should enforce that the field $u(x, y)$ vanishes in $\partial\Omega$ and then only in the parts of $\partial\Omega_1$ and $\partial\Omega_2$ on $\partial\Omega$. For this purpose the simplest way to proceed consists of introducing Dirichlet boundary conditions from a penalty term in the weak form (15). Then all the procedure just described remains valid and it includes an extra-term in the $x$ and $y$ problems (38)–(39) and (42)–(43) respectively.

### 4.4 Extension to parametric solutions

In the previous sections we considered for the sake of simplicity domain decomposition of solutions involving a separated representation of the space coordinates.

In this section we generalize the procedure just described for considering model parameters as extra-coordinates and then applying the domain decomposition technique just described for computing a quite general parametric solution in complex domains composed of simpler substructures.

For this purpose we consider the simplest parametric heat equation

$$k \cdot \Delta u - f = 0 \qquad (45)$$

with homogeneous boundary conditions.

Here $(x, y, k) \in \Omega_x \times \Omega_y \times \Omega_k$, with $\Omega_k \subset \mathbb{R}$ the interval in which the thermal conductivity $k$ takes its values. For the sake of simplicity and without loss of generality the source term $f$ is assumed constant. In this section the conductivity $k$ is viewed as a new coordinate defined in the interval $\Omega_k$. Thus, instead of solving the thermal model for different values of the conductivity parameter, we wish to solve at once a more general problem, the price to pay being an increase of the problem dimensionality. However, as the complexity of the PGD scales only linearly (and not exponentially) with the space dimension, consideration of the conductivity as a new coordinate still allows one to efficiently obtain an accurate solution.

The weighted residual form related to Eq. (45) reads:

$$\int_{\Omega_x \times \Omega_y \times \Omega_k} u^* \cdot (k \cdot \Delta u - f) dx dy dk = 0 \qquad (46)$$

for all test functions $u^*$ selected in an appropriate functional space.

The PGD solution is sought in the form:

$$u(x, y, k) \approx \sum_{i=1}^{N} X_i(x) \cdot Y_i(y) \cdot K_i(k) \qquad (47)$$

The separated representation constructor was deeply described in our former works (see [14] and the references therein). When considering domain decomposition, and in particular the Arlequin strategy, it suffices considering the procedure previously described but now applied on the parametric solutions of $u_1, u_2$ and $\lambda$:

$$\left\{ \begin{array}{l} u_1(x, y, k) \approx \sum_{i=1}^{N} X_1^i(x) \cdot Y_1^i(y) \cdot K_1^i(k) \\ u_2(x, y, k) \approx \sum_{i=1}^{N} X_2^i(x) \cdot Y_2^i(y) \cdot K_2^i(k) \\ \lambda(x, y, k) \approx \sum_{i=1}^{N} L_x^i(x) \cdot L_y^i(y) \cdot L_k^i(k) \end{array} \right. \qquad (48)$$

All the procedure just described can be easily generalized. The calculation of space functions require the solution of problems defined in the $x$ and $y$ coordinates, (38)–(39) and (42)–(43) respectively, where the coefficients now also include the integration in $\Omega_k$.

The equations to be solved for calculating the functions depending on the conductivity $k$ at the enrichment iteration $n$, $K_1^n(k)$, $K_2^n(k)$ and $L_k^n(k)$, will be

$$\int_{\Omega_k} (\alpha_1^{xy} K_1^* k K_1^n) dk + \int_{\Omega_k} \gamma_1^{xy} K_1^* L_k^n dk \\ = \int_{\Omega_k} K_1^* f_1^k(k) dk \qquad (49)$$

where functions $\alpha_1^{xy}$, $\gamma_1^{xy}$ and $f_1^k(k)$ come from integrals in $x$ and $y$. Analogously,

$$\int_{\Omega_k} (\alpha_2^{xy} K_2^* k K_2^n) dk - \int_{\Omega_k} \gamma_2^{xy} K_1^* L_k^n dk \\ = \int_{\Omega_k} K_2^* f_2^k(k) dk \qquad (50)$$

and

$$\int_{\Omega_k} \delta^{xy} L_k^* (K_1^n - K_2^n) dk = \int_{\Omega} L_k^* h^k(k) dk \qquad (51)$$
<figure><img src="image_7.png" /><figcaption>Fig. 6 FEM reference solution</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 PGD-Arlequin solution</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8 FEM and PGD-Arlequin solutions difference</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 PGD solution in domain $\Omega_1$</figcaption></figure>

## 5 Numerical results

In this section we present the results obtained when proceeding as described in the previous section. In order to quantify the error and the convergence and being unknown the exact solution, we computed the solution of the same problem in the whole domain by using the finite element method on a fine enough mesh. The reference solution is given in Fig. 6. The solution computed by using the separated representation combined with the Arlequin coupling strategy is depicted in Fig. 7. The difference of both solutions is depicted in Fig. 8 where the highest errors as expected locate in the coupling region. Figures 9 and 10 show the solution in each subdomain. The convergence in both domains is given in Figs. 11 and 12 that depict the evolution of the equations residual with the number of modes considered in separated representation. We can notice than tens modes are roughly enough for approximation the solutions with a reasonable degree of accuracy. For calculating more accurate solution meshes should be refined as discussed in [15]. Finally Figs. 13 and 14 group the main functions $(X_1^i(x), Y_1^i(y))$ and $(X_2^i(x), Y_2^i(y))$ respectively.

We considered two more complex geometries consisting of some substructures. The solution of these structures cannot be carried out by applying a direct space separated representation, being the only possibility performing a separated representation at the level of each subdomain and then couple all them by using an appropriate domain decomposition technique, the Arlequin technique in the present work. Again ten modes were enough for computing an accurate solution. The solutions accuracy was checked by solving the problems using standard finite element discretizations and considering these solutions as the reference ones. Figure 15 depicts the computed solutions in these two more complex geometries.

Finally we are addressing the solution of a parametric model related to the heat equation

$$k \cdot \Delta u + 1 = 0 \qquad (52)$$

defined in the domain $\Omega$ depicted in Fig. 16.
<figure><img src="image_11.png" /><figcaption>**Fig. 10** PGD solution in domain $\Omega_2$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Convergence analysis in domain $\Omega_1$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12** Convergence analysis in domain $\Omega_2$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13** Main modes involved in the PGD solution in $\Omega_1$</figcaption></figure>

As soon as the domain is decomposed in three subdomains with two overlapping regions, the procedure previously described can be applied, in order to finally obtain the parametric solution of the problem, that is, the temperature fields for each value of the conductivity in the interval $\Omega_k$. As soon as this solution is available, we can particularize it for different values of the conductivity as illustrated in Fig. 17.

We compared these particularized solutions with the references solutions obtained by solving the thermal model by using the finite element method in the whole domain for the four values of the thermal conductivity. Reference solutions and the ones resulting from the space-parameter separated representation coupled with the Arlequin domain decomposition method agree in minute.

## 6 Conclusion

We proved in our former works that problems defined in plate or shell geometries and needing a fully 3D solution can be efficiently solved by considering an in-plane-out-of-plane separated representation of the unknown fields. The use of such solution approximation allowed the calculation of fully 3D solutions by preserving a computational complexity

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 14 Main modes involved in the PGD solution in $\Omega_2$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15 Solution in two more complex domains</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16 Domain related to the parametric solution</figcaption></figure>

characteristic of 2D models. Even if such representation was extended for addressing slightly more complex geometries by invoking a kind of domain penalization, it definitively failed for addressing general plate or shell assembling.

A natural route for circumventing the just referred difficulties consists of considering domain decomposition strategy and then applying a separated representation to each substructure. The technique here considered, analyzed and validated in the one making use of the Arlequin strategy that can be coupled easily with a separated representation at the level of each substructure.

The coupling of parametric substructures where the separated representation involves the space and other material (or

geometrical) parameters, as just proved, opens new routes in the efficient simulation of complex structures.

## 7 Separated representation constructor: alternating direction strategy

Each iteration of the alternating direction scheme consists in the following two steps:

- Calculating $X_n^p(x)$ from $Y_n^{p-1}(y)$
In this case, the approximation reads

$$u^{n,p}(x, y) = \sum_{i=1}^{n-1} X_i(x) \cdot Y_i(y) + X_n^p(x) \cdot Y_n^{p-1}(y), \quad (53)$$

where all functions are known except $X_n^p(x)$. The simplest choice for the weight function $u^*$ in the weighted residual formulation (7) is

$$u^*(x, y) = X_n^*(x) \cdot Y_n^{p-1}(y), \quad (54)$$

which amounts to select the Galerkin weighted residual form of the Poisson equation.
<figure><img src="image_19.png" /><figcaption>**Fig. 17** Parametric solution particularized at four values of the thermal conductivity: $k = 1$ (top-left); $k = 3$ (top-right); $k = 5$ (down-left) and $k = 6$ (down-right)</figcaption></figure>

Injecting (53) and (54) into (7), we obtain

$$ 
\begin{aligned}
& \int_{\Omega_x \times \Omega_y} X_n^* \cdot Y_n^{p-1} \cdot \left( \frac{d^2 X_n^p}{dx^2} \cdot Y_n^{p-1} + X_n^p \cdot \frac{d^2 Y_n^{p-1}}{dy^2} \right) dx \cdot dy = \\
& - \int_{\Omega_x \times \Omega_y} X_n^* \cdot Y_n^{p-1} \cdot \sum_{i=1}^{n-1} \left( \frac{d^2 X_i}{dx^2} \cdot Y_i + X_i \cdot \frac{d^2 Y_i}{dy^2} \right) dx \cdot dy \\
& + \int_{\Omega_x \times \Omega_y} X_n^* \cdot Y_n^{p-1} \cdot f \, dx \cdot dy.
\end{aligned} 
\quad (55) $$

Here comes a crucial point: since all functions of y are known in the above expression, we can compute the following one-dimensional integrals over $\Omega_y$:

$$ 
\left\{ 
\begin{aligned}
\alpha^x &= \int_{\Omega_y} (Y_n^{p-1}(y))^2 dy \\
\beta^x &= \int_{\Omega_y} Y_n^{p-1}(y) \cdot \frac{d^2 Y_n^{p-1}(y)}{dy^2} dy \\
\gamma_i^x &= \int_{\Omega_y} Y_n^{p-1}(y) \cdot Y_i(y) dy \\
\delta_i^x &= \int_{\Omega_y} Y_n^{p-1}(y) \cdot \frac{d^2 Y_i(y)}{dy^2} dy \\
\xi^x &= \int_{\Omega_y} Y_n^{p-1}(y) \cdot f dy
\end{aligned} 
\right. 
\quad (56) $$

Equation (55) becomes

$$ \int_{\Omega_x} X_n^* \cdot \left( \alpha^x \cdot \frac{d^2 X_n^p}{dx^2} + \beta^x \cdot X_n^p \right) dx = $$

$$ 
\begin{aligned}
& - \int_{\Omega_x} X_n^* \cdot \sum_{i=1}^{n-1} \left( \gamma_i^x \cdot \frac{d^2 X_i}{dx^2} + \delta_i^x \cdot X_i \right) dx \\
& + \int_{\Omega_x} X_n^* \cdot \xi^x dx.
\end{aligned} 
\quad (57) $$

We have thus obtained the weighted residual form of a one-dimensional problem defined over $\Omega_x$ that can be solved (e.g. by the finite element method) to obtain the function $X_n^p$ we are looking for. Alternatively, we can return to the corresponding strong formulation

$$ \alpha^x \cdot \frac{d^2 X_n^p}{dx^2} + \beta^x \cdot X_n^p = - \sum_{i=1}^{n-1} \left( \gamma_i^x \cdot \frac{d^2 X_i}{dx^2} + \delta_i^x \cdot X_i \right) + \xi^x, \quad (58) $$

and solve it numerically by means of any suitable numerical method (e.g. finite differences, pseudo-spectral techniques, ...). The strong form (57) is a second-order ordinary differential equation for $X_n^p$. This is due to the fact that the original Poisson equation involves a second-order x-derivative of the unknown field $u$.

With either the weighted residual or strong formulations, the homogeneous Dirichlet boundary conditions $X_n^p(x = 0) = X_n^p(x = L) = 0$ are readily specified.

Having thus computed $X_n^p(x)$, we are now ready to proceed with the second step of iteration $p$.

### Calculating $Y_n^p(y)$ from the just-computed $X_n^p(x)$

The procedure exactly mirrors what we have done above. Indeed, we simply exchange the roles played by all relevant functions of $x$ and $y$.

The current PGD approximation reads

$$ u^{n,p}(x, y) = \sum_{i=1}^{n-1} X_i(x) \cdot Y_i(y) + X_n^p(x) \cdot Y_n^p(y), \quad (59) $$
where all functions are known except $Y_n^p(y)$.

The Galerkin formulation of (7) is obtained with the particular choice

$$u^*(x, y) = X_n^p(x) \cdot Y_n^*(y). \quad (60)$$

Then, by introducing (59) and (60) into (7), we get

$$\begin{aligned} & \int_{\Omega_x \times \Omega_y} X_n^p \cdot Y_n^* \cdot \left( \frac{d^2 X_n^p}{dx^2} \cdot Y_n^p + X_n^p \cdot \frac{d^2 Y_n^p}{dy^2} \right) dx \cdot dy \\ &= - \int_{\Omega_x \times \Omega_y} X_n^p \cdot Y_n^* \cdot \sum_{i=1}^{n-1} \left( \frac{d^2 X_i}{dx^2} \cdot Y_i + X_i \cdot \frac{d^2 Y_i}{dy^2} \right) dx \cdot dy \\ &+ \int_{\Omega_x \times \Omega_y} X_n^p \cdot Y_n^* \cdot f \, dx \cdot dy. \end{aligned} \quad (61)$$

As all functions of $x$ are known, the integrals over $\Omega_x$ can be computed to obtain

$$\left\{ \begin{aligned} \alpha^y &= \int_{\Omega_x} (X_n^p(x))^2 dx \\ \beta^y &= \int_{\Omega_x} X_n^p(x) \cdot \frac{d^2 X_n^p(x)}{dx^2} dx \\ \gamma_i^y &= \int_{\Omega_x} X_n^p(x) \cdot X_i(x) dx \\ \delta_i^y &= \int_{\Omega_x} X_n^p(x) \cdot \frac{d^2 X_i(x)}{dx^2} dx \\ \xi^y &= \int_{\Omega_x} X_n^p(x) \cdot f \, dx \end{aligned} \right. . \quad (62)$$

Equation (61) becomes

$$\begin{aligned} & \int_{\Omega_y} Y_n^* \cdot \left( \alpha^y \cdot \frac{d^2 Y_n^p}{dy^2} + \beta^y \cdot Y_n^p \right) dy = \\ & - \int_{\Omega_y} Y_n^* \cdot \sum_{i=1}^{n-1} \left( \gamma_i^y \cdot \frac{d^2 Y_i}{dy^2} + \delta_i^y \cdot Y_i \right) dy \\ & + \int_{\Omega_y} Y_n^* \cdot \xi^y dy. \end{aligned} \quad (63)$$

As before, we have thus obtained the weighted residual form of an elliptic problem defined over $\Omega_y$ whose solution is the function $Y_n^p(y)$. Alternatively, the corresponding strong formulation of this one-dimensional problem reads

$$\begin{aligned} & \alpha^y \cdot \frac{d^2 Y_n^p}{dy^2} + \beta^y \cdot Y_n^p \\ &= - \sum_{i=1}^{n-1} \left( \gamma_i^y \cdot \frac{d^2 Y_i}{dy^2} + \delta_i^y \cdot Y_i \right) + \xi^y. \end{aligned} \quad (64)$$

This again is an ordinary differential equation of the second order, due to the fact that the original Poisson equation involves second-order derivatives of the unknown field with respect to $y$. With both the weighted residual and strong formulations, the homogeneous Dirichlet boundary conditions $Y_n^p(y=0) = Y_n^p(y=L) = 0$ are readily specified.

We have thus completed iteration $p$ at enrichment step $n$.

## 8 Separated representation of the transformation Jacobians

The different Jacobian $\mathcal{J}_1^j$ related to the mapping $\mathcal{D}_1^j \to \mathcal{V}_1^j$ can be in the present example easily obtained for $j = 1, \dots, 6$:

$$\left\{ \begin{aligned} \mathcal{J}_1^1 &= \begin{pmatrix} 1 & 0 \\ 0 & 0.5 \end{pmatrix}, & \mathbf{x} \in (0, 1.5) \times (0, 1) \\ \mathcal{J}_1^2 &= \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, & \mathbf{x} \in (0, 1.5) \times (1, 2) \\ \mathcal{J}_1^3 &= \begin{pmatrix} 1 & 0 \\ 0 & 1.5 \end{pmatrix}, & \mathbf{x} \in (0, 1.5) \times (2, 3) \\ \mathcal{J}_1^4 &= \begin{pmatrix} 1 & y \\ 0 & x-1 \end{pmatrix}, & \mathbf{x} \in (1.5, 2) \times (0, 1) \\ \mathcal{J}_1^5 &= \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}, & \mathbf{x} \in (1.5, 2) \times (1, 2) \\ \mathcal{J}_1^6 &= \begin{pmatrix} 1 & 3-y \\ 0 & 3-x \end{pmatrix}, & \mathbf{x} \in (1.5, 2) \times (2, 3) \end{aligned} \right. \quad (65)$$

and similarly for the mapping $\mathcal{D}_2^k \to \mathcal{V}_2^k$, $k = 1, 2$. The associated determinants result:

$$\left\{ \begin{aligned} \det(\mathcal{J}_1^1) &= 0.5 \\ \det(\mathcal{J}_1^2) &= 1 \\ \det(\mathcal{J}_1^3) &= 1.5 \\ \det(\mathcal{J}_1^4) &= x-1 \\ \det(\mathcal{J}_1^5) &= 1 \\ \det(\mathcal{J}_1^6) &= 3-x \end{aligned} \right. \quad (66)$$

The Jacobian of the inverse transformations also accepts a fully separated representation because as just proved the
determinant of the Jacobian only involves functions of x. For example, when considering $(\mathcal{J}_1^6)^{-1}$

$$ (\mathcal{J}_1^6)^{-1} = \frac{1}{\det(\mathcal{J}_1^6)} \begin{pmatrix} 3-x & y-3 \\ 0 & 1 \end{pmatrix} = \begin{pmatrix} 1 & \frac{y-3}{3-x} \\ 0 & \frac{1}{3-x} \end{pmatrix} \quad (67) $$

that proves its fully separability.

In the most general case in which the determinant involves a two-terms sum where each term consists of the product of a function of x times a function of y, the components of $(\mathcal{J}_1^j)$ have not a direct separated representation [4]. Approximated separated representations can be obtained by invoking the SVD (or the HOSVD) as indicated before.

After being calculated the Jacobian of the direct and inverse transformations of each quadrilateral we can compact the notation by introducing the characteristic function of each quadrilateral $\mathcal{V}_1^j$ and $\mathcal{V}_2^k$:

$$ \chi_1^j(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \mathcal{V}_1^j = (x_1^{j-}, x_1^{j+}) \times (y_1^{j-}, y_1^{j+}) \\ 0 & \text{elsewhere} \end{cases} \quad (68) $$

and

$$ \chi_2^k(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \mathcal{V}_2^k = (x_2^{k-}, x_2^{k+}) \times (y_2^{k-}, y_2^{k+}) \\ 0 & \text{elsewhere} \end{cases} \quad (69) $$

from which we can write the Jacobian of the mappings $\Omega_1 \to \mathcal{T}_1$ and $\Omega_2 \to \mathcal{T}_2$, respectively $\mathcal{J}_1$ and $\mathcal{J}_2$:

$$ \mathcal{J}_1(\mathbf{x} \in \mathcal{T}_1) = \sum_{j=1}^{6} \chi_1^j(\mathbf{x}) \cdot \mathcal{J}_1^j(\mathbf{x}) \quad (70) $$

and

$$ \mathcal{J}_2(\mathbf{x} \in \mathcal{T}_2) = \sum_{k=1}^{2} \chi_2^k(\mathbf{x}) \cdot \mathcal{J}_2^k(\mathbf{x}) \quad (71) $$

and similarly for the inverse transformations

$$ (\mathcal{J}_1)^{-1}(\mathbf{x} \in \mathcal{T}_1) = \sum_{j=1}^{6} \chi_1^j(\mathbf{x}) \cdot (\mathcal{J}_1^j(\mathbf{x}))^{-1} \quad (72) $$

and

$$ (\mathcal{J}_2)^{-1}(\mathbf{x} \in \mathcal{T}_2) = \sum_{k=1}^{2} \chi_2^k(\mathbf{x}) \cdot (\mathcal{J}_2^k(\mathbf{x}))^{-1} \quad (73) $$

The mapping determinants can be also expressed in the compact form

$$ \det(\mathcal{J}_1(\mathbf{x} \in \mathcal{T}_1)) = \sum_{j=1}^{6} \chi_1^j(\mathbf{x}) \cdot \det(\mathcal{J}_1^j(\mathbf{x})) \quad (74) $$

and

$$ \det(\mathcal{J}_2(\mathbf{x} \in \mathcal{T}_2)) = \sum_{k=1}^{2} \chi_2^k(\mathbf{x}) \cdot \det(\mathcal{J}_2^k(\mathbf{x})) \quad (75) $$

## References

1. Ammar A, Mokdad B, Chinesta F, Keunings R (2006) A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. J Non Newton Fluid Mech 139:153–176

2. Ammar A, Mokdad B, Chinesta K, Keunings R (2007) A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: Transient simulation using space-time separated representation. J Non Newton Fluid Mech 144:98–121

3. Ammar A, Chinesta F, Diez P, Huerta A (2010) An error estimator for separated representations of highly multidimensional models. Comput Methods Appl Mech Eng 199:1872–1880

4. Ammar A, Huerta A, Chinesta F, Cueto E, Leygue A (2014) Parametric solutions involving geometry: a step towards efficient shape optimization. Comput Methods Appl Mech Eng 68(1):178–193

5. Ben Dhia H (1998) Multiscale mechanical problems: the Arlequin method. CR Acad Sci 326:899–904

6. Ben Dhia H, Rateau G (2005) The Arlequin method as a flexible engineering design tool. Int J Numer Methods Eng 62:1442–1462

7. Ben Dhia H (2008) Further insights by theoretical investigations of the multiscale Arlequin method. Int J Multiscale Comput Eng 6(3):215–232

8. Bernardi C, Maday Y, Patera AT (1994) A new non conforming approach to domain decomposition: The mortar element method. In H. Brezzi et al (eds) Nonlinear partial differential equations and their applications. Paris, pp 13–51

9. Bognet B, Leygue A, Chinesta F, Poitou A, Bordeu F (2012) Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity. Comput Methods Appl Mech Eng 201:1–12

10. Bognet B, Leygue A, Chinesta F (2012) On the fully 3D simulation of thermoelastic models defined in plate geometries. Eur J Comput Mech 21(1–2):40–51

11. Bognet B, Leygue A, Chinesta F (2014) Separated representations of 3D elastic solutions in shell geometries. Adv Model Simul Eng Sci 1(1):1–34

12. Cancès E, Defranceschi M, Kutzelnigg W, Le Bris C, Maday Y (2003) Computational Quantum Chemistry: a primer. Handbook of Numerical Analysis, vol X. Elsevier, pp 3–270

13. Chinesta F, Leygue A, Bognet B, Ghnatios Ch, Poulhaon F, Bordeu F, Barasinski A, Poitou A, Chatel S, Maison-Le-Poec S (2014) First steps towards an advanced simulation of composites manufacturing by Automated Tape Placement. Int J Mater Form 7(1):81–92

14. Chinesta F, Leygue A, Bordeu F, Aguado JV, Cueto E, Gonzalez D, Alfaro I, Ammar A, Huerta A (2013) PGD-based computational vademecum for efficient design, optimization and control. Arch Comput Methods Eng 20:31–59
15. Chinesta F, Keunings R, Leygue A (2013) The Proper Generalized Decomposition for advanced numerical simulations. A primer, Springerbriefs, Springer

16. Farhat C, Roux FX (1991) A method of finite element tearing and interconnecting and its parallel solution algorithm. Int J Numer Meth Eng 32:1205-1227

17. Farhat C, Mandel J (1998) The two-level FETI method for static and dynamic plate problems - part I: an optimal iterative solver for biharmonic systems. Comput Methods Appl Mech Eng 155:129-152

18. Ghnatios CH, Ammar A, Cimetiere A, Hamdouni A, Leygue A, Chinesta F (2012) First steps in the space separated representation of models defined in complex domains. ASME 2012: 11th Biennial conference on engineering systems design and analysis -ESDA 2012-, Nantes

19. Gonzalez D, Ammar A, Chinesta F, Cueto E (2010) Recent advances in the use of separated representations. Int J Numer Methods Eng 81(5):637-659

20. Ladeveze P, Lorong Ph (1992) A large time increment approach with domain decomposition technique for mechanical non linear problems. Comput Methods Appl Sci Eng INRIA 569-578

21. Ladeveze P, Dureisseix D (1998) A 2-level and mixed domain decomposition approach for structural analysis. Contemp Math 218:246-253

22. Ladeveze P, Nouy A (2003) On a multiscale computational strategy with time and space homogenization for structural mechanics. Comput Methods Appl Mech Eng 192(28-30):3061-3087

23. Ladeveze P, Néron D, Gosselet P (2007) On a mixed and multiscale domain decomposition method. Comput Methods Appl Mech Eng 96:1526-1540

24. Ladeveze P, Passieux J-C, Néron D (2010) The latin multiscale computational method and the proper generalized decomposition. Comput Methods Appl Mech Eng 199(21-22):1287-1296

25. Ladeveze P, Chamoin L (2011) On the verification of model reduction methods based on the proper generalized decomposition. Comput Methods Appl Mech Eng 200:2032-2047

26. Leygue A, Chinesta F, Beringhier M, Nguyen TL, Grandidier JC, Pasavento F, Schrefler B (2013) Towards a framework for non-linear thermal models in shell domains. Int J Numer Methods Heat Fluid Flow 23(1):55-73

27. Maday Y, Mavriplis C, Patera AT (1989) Nonconforming mortar element methods: application to spectral discretizations. Domain decomposition methods, Los Angeles, CA, 1988. SIAM, Philadelphia, pp 392-418

28. Schwarz HA (1890) Gesammelte mathematische abhandlungen, vol 2. Springer, Berlin