Domain decomposition involving subdomain separable space representations for solving parametric problems in complex geometries

A domain decomposition technique combined with an enhanced geometry mapping based on the use of NURBS is considered for solving parametrized models in complex geometries (non simply connected) within the so-called proper generalized decomposition (PGD) framework, enabling the expression of the solution in each subdomain in a fully separated form, involving both the space and the model parameters. NURBS allow a compact and powerful domain mapping into a fully separated reference geometry, while the PGD allows recovering an affine structure of the problem in the reference domain, facilitating the use of the standard PGD solver for computing the parametric solution in each subdomain first, and then by enforcing the interface transmission conditions, in the whole domain.


Introduction
Proper generalized decomposition (PGD) is based on the use of separated representations, fact that alleviates the computational complexity when solving multidimensional problems by reducing them to a sequence of lower dimensional problems. PGD was successfully applied for addressing four problem typologies: • Transient models can be decomposed into a series of space and time problems, leading to a non-incremental integrator where the solution history is computed simultaneously instead of incrementally as most of time integrators perform. The interested reader can refer to [7] and the numerous references therein; • Parametric models can be decomposed into a series of space, time and parameter problems, enabling the calculation of the whole solution history for any value of the parameters involved in the model and here assumed as model extra-coordinates, as described in [8,10] and the numerous references therein; • The solution of models defined in a degenerated domain (like plate, shells, beam, ... geometries) where at least one of its characteristic dimensions is much smaller than the others, can be successfully accomplished within the separated representation framework, where each dimension is treated sequentially, weakly coupled with the others, reducing the 3D problem to a sequence of one-dimensional ones (or a series of two and one dimensional problems) which alleviates the meshing issues usually encountered when performing 3D analyses in this kind of domains [7]. • Problems defined in non-degenerated but non-separable domains that could become separable in a domain decomposition framework.
This work concerns the fourth problem typology, that is, the one related to spaceseparation with domain decomposition.
In the past it was claimed that this kind of decomposition works well when the domain can be expressed from the cartesian product of one-dimensional domains, or a twodimensional and a one-dimensional ones as it is the case when addressing plates and shell geometries, laminates or extruded profiles. More complex geometries, non-separable in a direct way, can be separated by using appropriate mappings [2,13,14] or penalty techniques [13], and thus, fully space-separated representation can be applied to reduce the computational cost of solving 3D problems to the one characteristic of the solution of 2D or 1D problems.
However, the proposed mappings generality remains quite limited when the domain shape complexity increases. In [1] we proposed a two levels discretization, the coarser to parametrize the shape and the finest for approximating the problem solution, being the main remaining challenges: (i) the accurate geometry parametrization, looking for its most compact representation when that geometry is being considered as a model parameter (a very valuable route in shape optimization); and (ii) the affine expression of the problem defined in the reference (fully separable) domain to ensure the efficiency of the separated representation constructor (PGD solver).
For alleviating the just referred issues domain decomposition techniques could be envisaged. Domain decomposition techniques have been widely employed in a diversity of applications. Some techniques are based on the subdomains overlapping, other consider non-overlaping coupling strategies. Among the overlapping-based strategies we can mention the classical Schwartz [27] or the more recent and versatile Arlequin technique and its variants [3][4][5]. Among the non-overlapping ones we can mention the mortar-based techniques [6,23], the FETI strategy and its variants [11,12] and the LATIN-multiscale and LATIN-multidomain approaches widely developed by Ladeveze and coauthors [18][19][20][21] [22]. These works considered multi-domain 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. Authors proposed different techniques based on domain decomposition, within the Arlequin coupling paradigm [25] and within a more standard decomposition (without overlapping) technique, assembling parametric patches [15].
Recently, in [16] we proposed a powerful mapping able to represent complex CAD geometries (employing NURBS) and able to efficiently parametrize complex mappings while recover efficiently an affine decomposition of the differential operators involved in the model. However, the presented technology only applied in simply connected domains, that is, for example in 2D domains without holes.
In the present work the technique proposed in [16] is extended for addressing non simply connected domains. For that purpose, we consider the mapping proposed in [16] in tandem with a domain decomposition technique, where the required continuity is ensured at the patches interface level. The proposed technique will be applied in some numerical experiments to prove its potential.

Geometrical mapping and domain decomposition
To use the full potential of the PGD technique for solving parametric and high dimensional problems defined in non-regular geometries (non rectangular or hexahedral) it is necessary to map the problem from the non-regular physical domain, x , into a regular (in the sense of separable) computational domain, ξ . A schematic geometric mapping is shown in Fig. 1. It can be noticed that the geometric mapping should ensure a separable nature of the problem. Recently, the authors have proposed a general and easily separable mapping in conjunction with the NURBS basis functions [16] for the geometry description. This technique makes possible using the PGD for solving parametric problems in 2D and 3D non-separable domains.
This technique can be generalized for solving problems in non simply connected domains as the one illustrated in Fig. 2 containing a hole, by transforming it into a simplyconnected domain. For that purpose, one could cut the domain creating a new boundary connecting the inner and outer boundaries and then applying in it the continuity constraints when solving the problem, for guaranteeing the solution periodicity across this   artificial boundary. More details regarding the enforcement of the continuity constraints will be presented later in subsequent sections.
Although the above mentioned procedure could be used to convert any domain with one internal hole into a simply connected domain, the mapping could become very complex in some cases affecting the solution efficiency. It also must be noted that in the case of domains with many holes this technique needs more than a simple single cut, making difficult the proposed procedure.
To overcome this difficulty, it is proposed in the current work to use more than one computational domain while applying the PGD in each computational domain (subdomain). In this approach, the physical domain would be decomposed into some relatively simple subdomains or patches and then each patch mapped into a separable computational domain.
To preserve the continuity of the solution on the subdomain interfaces, appropriate continuity constraints must be applied between adjacent patches expressed in their computational domains.
A schematic representation of this multi-patch technique is shown in Figs. 3 and 4 for domains containing one and two internal holes, respectively. Such domain decomposition is possible for other topologies addressed later. Details concerning the continuity enforcement between adjacent patches will be presented in the following sections.

Problem formulation
For the sake of simplicity and without loss of generality in the present work we evaluate the proposed multi-domain technique by considering the nonlinear heat equation with a temperature dependent thermal conductivity, expressed from: where u is the unknown temperature field defined in physical domain x , k(u) is the temperature dependent thermal conductivity and f is the distributed heat source. The nonlinear models will be addressed by using an adequate linearization. In the context of the PGD addressed later, the simplest route consists in evaluating the nonlinear terms, the conductivity in the present case, from the temperature field computed at the previous iteration.
The physical domain, x , results from the union of non-overlapping subdomains or patches s x .
where N p is the total number of patches. The weighted residual form of Eq. (1) reads where w and u are the test and trial functions respectively, and B and L the associated bi-linear and linear forms. The problem can be expressed in the variational form The unknown field u(x) must be obtained in such a way that minimize the functional I in Eq. (4). In the proposed multi-patch technique, a penalty method is used to enforce the Dirichlet boundary conditions, by considering the extended variational form whereū is the prescribed value of the field on the Dirichlet boundary D and α is the associated penalty parameter. By minimizing the functional it results wherek represents the conductivity linearization andq is the prescribed heat flux on the Neumann boundary N . In Eq. (6), u * is the first variation of the field variable u.

Solution procedure
The use of NURBS enables the efficient description of general curves, surfaces and volumes. According to the notation introduced in Appendix A, the mapping from the physical to the computational domain in the 2D case reads where are R pq ij (ξ) are the NURBS-based shape functions, P s ij refers to the control points (the vertices of the so-called control net) related to the s-th patch in the physical domain. Eq. (7) expresses the geometrical mapping, with ξ in the computational domain s ξ to a point x in physical domain s x . The previous equation allows not only transforming the domain, but also transforming the differential operators and computing the transformation Jacobian (see Appendix A).
Even if the computational domain in which the physical one is mapped, is fully separable by construction, the differential operators are not, because of the fact that the just referred shape functions involved in the geometry mapping are not separated. Thus, to obtain a separated representation of a generic function g(x) into the computational domain, it suffices looking for the separated form with M j (ξ j ) the vector of approximation functions involving the j-th coordinate and G ji is the vector of coefficients of the i-th mode in the j-th coordinate direction. N D is the total number of problem dimensions (in the present case N D = 2 if the problem is 2D or N D > 2 when some problem parameters are considered as model extra-coordinates) and N G is the number of modes needed for reaching the desired accuracy (as described in Appendix B).
Thus, the so-called Separated Approximate Representation -SAR-proceeds by enforcing in a weak form the equality g(ξ) = g(x(ξ)) in s ξ , with g(ξ) expressed in the separated form (8), following the procedure deeply addressed in our former works [9,16], to finally obtain the coefficients G ji .
This procedure, the just described SAR, is then applied for expressing both the Jacobian and the differential operator transformation (see Appendix C), in an affine form, making possible the application of the standard PGD procedure.
Thus, the separated representation of the temperature field u(ξ) in the computational subdomain s ξ reads at the n-enrichment (the dependence on the subdomain considered, referred by the upper-script • s , is omitted for the sake of notational simplicity) where U jn , j = 1, 2, . . . , N D must be calculated from the problem weak form by using the standard fixed-point alternated directions fixed point algorithm [9]. When updating the function concerned by the d-coordinate (in the present case d = 1 or d = 2), the test function reads: that allows updating coefficients U dn . Nonlinear models can be addressed by computing nonlinear terms at the current solution (for instance u n−1 (ξ)) whose affine description can be recovered by applying again the SAR. For additional details the interested reader can refer to [16].

Enforcing continuity constraints
As stated before, continuity constraints must be enforced between adjacent patches to enforce the continuity of the field. The continuity constraints may be defined between opposite sides of a single patch or between sides of adjacent patches. In the former case, the standard fixed point iterative method involved in the PGD constructor, can be applied for solving the resulting system of nonlinear algebraic equations.
Application of the continuity constraints between the adjacent edges of a multi-patch topology is more complex and the standard fixed point iterative approach is not longer applicable. In this case, a multi-action fixed point iterative approach is proposed in the present work to solve the system of nonlinear algebraic equations. These techniques are discussed in the following subsections.

Single patch PGD with continuity constraints
As previously discussed a valuable approach to solve the problem in a domain containing an internal hole consists of applying a cut line and transform the non simply connected region into a simply connected one as shown schematically in Fig. 2.
A continuity constraint then should be enforce to ensure the continuity of the field along that cut-line. In this approach, the problem can be solved using a single patch (one computational domain). Consider, for instance, that the continuity constraint is defined between the edges parallel to ξ axis, that according to Eq. (9) implies some constraints to be applied in the modes concerning the direction η.
In other words, when operating with M T 2 U 2n , the constraint reads: When using standard piecewise linear functions for expressing M 2 , verifying the Kroenecker delta property, the previous constraint implies that the first and last degree of freedom coincides. That is, in N 2 nodes are used for approximating the modes involving the η coordinate, the previous continuity constraint simply reads that can be directly enforced when looking for the coefficients U 2n .

Multi-patch PGD with continuity constraints
When more than one computational domain is considered, the continuity constraint and fixed point iteration schema becomes more complex. To explain the basic idea and the procedure, three sample topologies for computational domains are considered. More complex cases are considered in subsequent sections.

First topology
For the first sample topology, consider the two patches shown in Fig. 5. There is a single continuity constraint and is shown in that figure. Consider in the n-th enrichment step, we would like to find U s dn , where s is the patch number, and d is the space direction (d = 1 for ξ and d = 2 for η).
The continuity constraint (see Fig. 5) indicates that the unknown field at the edge BC of patch 1 must be equal to the field at the edge BC of patch 2. Therefore, the modes involving the η direction in patch 1 must be the same as the modes in this same direction in patch 2. In addition, the value of the modes involving the ξ direction at the common boundary must match. To make it more clear, refer to Fig. 6. This figure shows schematic representation of modes in ξ and η directions for two patches 1 and 2. To satisfy the continuity between two patches, Fig. 6a, b show that the modes of ξ should have same value at the patches interface. On the other hand, Fig. 6c, d show that the modes of η should be same function for both patches. Mathematically, the continuity constraint for The PGD system of algebraic equations regarding direction d of patch s at the n-th enrichment step reads (see [16]): To solve the equations and obtain the n-th enrichment modes, a modified fixed point iteration approach is used here considering the continuity constraint given in Eq. (13). This modified fixed point schema consists of the three steps: 1. Regarding the continuity constraint given in the first row of Eq. (13), the system to be solved reads 2. Regarding the continuity constraint given in the second row of Eq. (13), the system to be solved reads The Eqs. (15)- (17) are located in the inner loop of fixed point algorithm. It is clear that the first row of Eq. (13) must be considered when solving Eq. (15) and it consists of simple matrix operations. Each system of Eqs. in (15)-(17) is solved by a direct system solver.

Second topology
The second sample topology consists of three patches shown in Fig. 7. The continuity constraints are also shown in that figure. Consider in the n-th enrichment step, we would like to find U s dn , where s is the patch number, and d is the space direction (d = 1 for ξ and d = 2 for η).
The continuity constraint 1 (in reference to Fig. 7) indicates that the unknown field at the edge BC of patch 1 must be equal to the field at the edge BC of patch 2. Therefore, the modes involving the η direction in patch 1 must be the same as the ones in this same direction in patch 2. In addition, the value of the modes involving the ξ direction at the common boundary must match, i.e.: Cont. Cons. 1 : The same rationale applies for enforcing the continuity constraint 2 (see Fig. 7). To make it more clear, the modes involving ξ direction in patch 2 must be the same as the modes in ξ direction in patch 3. And, the modes involving η direction must match at the common interface of patches 2 and 3. That reads: The PGD system of algebraic equations regarding direction d of patch s at the n-th enrichment step reads The modified fixed point schema consists of the three steps: 1. Regarding the continuity constraint 1 given in Eq. (18), the system to be solved reads 2. Regarding the continuity constraint 2 given in Eq. (19), the system to be solved reads 3. If other parameters are involved as extra-coordinates, that can be represented by a generic d, d > 2, the associated degrees of freedom results from The systems of Eqs.  (21)-(23) is solved by a direct system solver but the whole process is iterative.

Third topology
For the third sample topology, we consider a computational domain composed of the four patches shown in Fig. 8 that shows the four continuity constraints.
By using the notation that was used previously, the continuity of the filed along all adjacent edges, the continuity constraints 1 to 4 are given in Table 1. Note that the operator Flip(·) that appears in some terms of Table 1 acts on a vector, by reversing the elements of the vector. It also acts in a similar way on the associated matrices. For example concerning the second constraint, η increases in patch 2 while the associated coordinate, ξ in patch 3 decreases.
Taking into account the usual PGD update given by Eq. (20) and the continuity constraints given in Table 1, the procedure reads: 1. Compute the less-constrained modes where only the constraints in the last column of Table 1 are enforce in a direct manner, ⎡ 2. The most constrained modes are then calculated, in fact only one because the first column in Table 1 implies that it propagates along all the connected patches. This mode reads 3. If other parameters are involved as extra-coordinates, that can be represented by a generic d, d > 2, the associated degrees of freedom results from Similar to the previous topology, the systems of Eqs.

Numerical examples
In this section some numerical examples are presented to evaluate the applicability and potential of the proposed multi-patch separable PGD technique.

First case study
The first numerical example concerns the Poisson problem in a circular domain containing an internal hole. The problem domain, its dimensions and the control net that is used for representing the physical domain from NURBS are shown in Fig. 9. The domain is transformed into a simply-connected domain by using a cut-line, the segment AC, shown in Fig. 10(left).
The computational domain and the continuity constraint are shown in Fig. 10(right). To construct the NURBS basis, the order of the B-Spline basis functions is of two concerning the dimension ξ and one concerning η. The knot vectors κ ξ and κ η (introduced in Appendix A) are given by and respectively. Homogeneous Dirichlet boundary conditions are enforced on both boundaries, the inner and outer, that is, on the boundaries ABA and CDC of the computational domain depicted in Fig. 10.
The heat source f (x, y) in Eq. (1) reads: This particular choice leads to the exact solution The separated representation of the physical coordinates x(ξ , η) and y(ξ , η) (according to Appendix C) involve a single mode for the former and two modes for the last, represented in Figs. 11 and 12, for x and y, respectively.
The modes involved in the separated representation of the solution involved 41 nodes for the modes affecting the ξ coordinate and 21 the ones concerning η.   The first five modes involved in the separated representation of u(ξ , η) are shown in Fig.  13. The contours of the reconstructed field in the physical and computational domains are shown in Fig. 14. The computed PGD solution is compared with exact one in Fig.  15 along the line x = 0 for 1 < y < 4, proving the excellent accuracy of the computed solution.

Second case study
The second case study concerns a domaine containing two holes. The problem domain is divided into 8 subdomains or patches. The dimensions of the physical domain and the control nets related to the NURB representation of the geometry are shown in Fig. 16.  The patches in the physical and computational domains and the continuity constraints are shown in Figs. 17 and 18, respectively. In the present numerical study all the patches have the same structure, therefore, the same knot vectors are used for all patches, in particular A nonzero Dirichlet boundary condition is enforced on the boundary of both internal holes, with u = 100 and u = 10 on the left and right holes respectively, that will be enforced by using the penalty technique discussed before. In the present case, a temperature dependent thermal conductivity is considered, according to k(u) = 1 + 0.04u. No thermal source is considered, i.e., f (x, y) = 0.
The 9 continuity constraints which are shown in Fig. 18 are expressed in detail in Table  2. These continuity constraints will be considered as previously discussed.
The three steps procedure previously described is here applied again, where the third step does not apply because there are not extra-coordinates to be addressed:    Table 2 Continuity constraints to be enforced in the second use case, with N 1 = N 2 ≡ N for all the subdomains Continuity constraint Equations 1. The modes that concerns each subdomain, with as only constraints the ones reported in the last column of Table 2, leads to the system ⎡ 2. The second step concerns the calculation of the modes that propagate throughout different patches for ensuring the continuity (first column in Table 2). It is easy to see that in the present case these modes verify the system In the present case study the 1D discretizations of modes involving the ξ and η coordinates is done by using 21 nodes in each direction. The contours of the reconstructed solution, in the physical u(x, y) and computational u(ξ , η) domains, are shown in Figs. 19 and 20, respectively.
The first six modes involved in the solution u(ξ , η), are shown for the eight patches in Fig. 21. To evaluate the PGD results, the problem is also solved by using the finite element method with linear quadrilateral elements. The solution on the upper edge of the domain,

Third case study
In the last case study we are including parameters as model extra-coordinates. In particular the domain geometry is parameterized here and the field variable is not only a function of coordinates ξ and η but also it is a function of few extra-parameters describing the domain shape.
We solve the Laplace equation, in the domain shown in Fig. 23, where the internal hole has a shape represented by the four digits NACA profile that involves three parameters M, P and T [24]. Thus, the problem involves 5 coordinates, the space ones, ξ and η, and the three extra-coordinates M, P and T , that within the PGD rationale, its solution will reduce to a sequence of five one-dimensional problems.
The physical domain consists of the internal boundary (where a null normal gradient is enforced) and an outer boundary (ellipse) on which the solution in enforced (Dirichlet boundary condition), u(x, y) = x. A single computational domain (single patch) is used here to map the physical domain. The geometry mapping and continuity constraint are shown in Fig. 23. Both internal and external boundaries are described by using NURBS, each involving 31 control points.  The continuity constraint is very simple as previously discussed, it reduces to The ranges of the parameters describing the shape are: M ∈ [0, 5]%, P ∈ [35, 50]% and T ∈ [15, 30]%. The 1D discretization is done by using 11 nodes in each directions ξ , M, P and T and 31 nodes in the η direction. The contours of the reconstructed solution, u(ξ , η, M, P, T ), in the physical domain are shown in Fig. 24 for four different choices of the parameters defining the internal shape. The first six modes involved in the solution separated representation are shown in Fig. 25.

Conclusions
The present paper proposed a valuable methodology for addressing complex geometries non simply connected, while enabling a fully separated representation with respect to the   space coordinates as well as on some problem parameters, considered as model extracoordinates within the PGD rationale.
The proposed technique combines domain decomposition, a geometrical description based on the use of NURBS and the use of the PGD for recovering an affine representation after the mapping that transforms complex geometries into fully separated ones.
Model extra-coordinates can be easily included in the formalism, that only needs make special care on the continuity enforcement across the subdomains interfaces.

Appendix A: NURBS-based geometry description
NURBS enables accurately representing the geometry of curves, surfaces and volumes, being nowadays a usual technology employed in CAD. For the sake of completeness this appendix revisits the main concepts involved in the NURBS construction. For additional details refer to [26].
NURBS result from weighted combinations of B-spline functions. To define a set of n B-spline functions of order p in a univariate parametric space ξ ∈ [0, 1], the knot vector κ ξ is defined as follows.
Denote N ap (ξ ) as the B-spline basis function of order p in the a-th knot span ξ ∈ [ξ a , ξ a+1 ]. The following recursive equations can be used to compute the univariate Bspline basis function N ap (ξ ) [26] N a0 (ξ ) = 1 ξ a ≤ ξ < ξ a+1 0 otherwise , for p = 0 (36) N ap (ξ ) = ξ − ξ a ξ a+p − ξ a N a(p−1) + ξ a+p+1 − ξ ξ a+p+1 − ξ a+1 N (a+1)(p−1) , for p > 0 (37) A proper choice of the knot vector allows obtaining rich behavior of the basis functions and enough flexibility to describe complex geometries. The NURBS basis functions can be obtained using a rational weighted sum of the B-splines basis functions. With w i the weight, the univariate NURBS basis function, R p i (ξ ) reads [26] The bivariate NURBS basis functions R pq ij (ξ , η) in 2D domains can be obtained using tensor product of univariate basis functions.
In the above equation p and q are the order of B-splines in directions ξ and η, respectively and ξ = (ξ , η) is the coordinates vector in the computational domain. In this equation, w ij denotes the geometry related weight parameter. The knot vector κ η applying in the coordinates η read being n and m the number of basis functions in directions ξ and η respectively. Using the NURBS basis functions just described, the bivariate physical domain x can be described from where, P refers to the control points (the vertices of the so-called control net) in the physical domain. Equation (41) represent a one to one relation that maps any point ξ in the computational domain ξ to a point x in physical domain x . Now, using the notation (x, y) ≡ (x 1 , x 2 ) and (ξ , η) ≡ (ξ 1 , ξ 2 ), for facilitating the mapping description compactness, the terms involved in the transformation of the differential operator from x to ξ read ∂x ∂ξ a = n i=1 m j=1 ∂R pq ij (ξ) ∂ξ a P ij , a,b= 1, 2.
More details and derivations regarding basic concepts and computations of NURBS basis functions and their derivatives could be found in [26]. Note that the geometric mapping given in Eq. (41) and its derivatives in Eq. (42) leads to a non-separable mapping in the sense of the PGD technique, compromising the effectiveness of the PGD solver. To overcome this difficulty, we apply a Separated Approximate Representation (SAR) of the Jacobian of the transformation [16], revisited in Appendix B. separated representation of the coordinates mapping, x(ξ), as follows: and where, N X and N Y are the number of modes required to reach the desired level of accuracy in the SAR of x(ξ) and y(ξ), respectively. Using the same rationale, we can derive the SAR of Jacobian determinant, |J (ξ)|, and the all 4 elements of the transformation derivatives tensor, h(ξ), as follows: More details regarding separated geometry mapping and its performance in parametric solution of field problems are given in [16].