Enhanced parametric shape descriptions in PGD-based space separated representations

Space separation within the Proper Generalized Decomposition—PGD—rationale allows solving high dimensional problems as a sequence of lower dimensional ones. In our former works, different geometrical transformations were proposed for addressing complex shapes and spatially non-separable domains. Efficient implementation of separated representations needs expressing the domain as a product of characteristic functions involving the different space coordinates. In the case of complex shapes, more sophisticated geometrical transformations are needed to map the complex physical domain into a regular one where computations are performed. This paper aims at proposing a very efficient route for accomplishing such space separation. A NURBS-based geometry representation, usual in computer aided design—CAD—, is retained and combined with a fully separated representation for allying efficiency (ensured by the fully separated representations) and generality (by addressing complex geometries). Some numerical examples are considered to prove the potential of the proposed methodology.


Introduction
A generic problem in physics consists of a differential operator acting on the so-called unknown field. In what follows for the sake of clarity we will assume a scalar unknown field depending on space x and time t, as well as on a series of parameters, grouped in vector μ, related to the considered physics, the loading or the domain in which the problem is defined. Thus, the problem reads where x ∈ x (μ) ⊂ R 3 , t ∈ t ⊂ R and the parameteres μ ∈ μ ⊂ R P . For approximating the unknown function, and in absence of a priori knowledge, any multi-purpose polynomial approximation basis can be used. The Finite Element Method-FEM-uses simple approximation bases with interpolative properties and compact support. To approximate very general functions the approximation basis must be rich enough, with the computational efficiency impact, in particular when addressing nonlinear tran-sient models, where the computational complexity scales with the space and time mesh resolution.
For alleviating such complexity, the first route consists in reducing the number of degrees of freedom. Thus, the so-called Proper Orthogonal Decomposition-POD- [10], extract from some offline computations the most significant modes φ i (x) for representing those solutions.
A generalization consists of expressing u(x, t) from the finite sum with both functions T i (t) and X i (x) assumed unknown and calculated on the fly by using an appropriate linearization [25]. The interested reader can refer to [3,9,26,29,32,33] and the references therein for additional references. Equation (2) expresses that the solution of M problems in space for calculating functions X i (x) and the same number of one-dimensional problems defined in time for computing function T i (t) suffices for computing it, in fact few more because of the nonlinearity of the solution procedure.
PGD has also been used with the attempt of alleviating the complexity of the problems involving the space coordinates. For that purpose, when the spatial domain x can be fully or partially decomposed, i.e. x = x × y × z or x = xy × z respectively, the separated representation based approximations read respectively that allows calculating the 3D solution from a sequence of 1D problems, and The former is specially suitable when the problem is defined in hexahedral domains, whereas the later applies in degenerated domains like plates [5], shells [6] or extruded domains (e.g. beams) [7,27]. In both cases the computing time savings could be impressive, with a complexity scaling with standard 1D discretization for the former and 2D for the last, while computing a fully 3D solution [9]. The same approach was extensively considered in structural plate and shell models in [31,[36][37][38], squeeze flows [13] and electromagnetism problems [35] defined in laminates.
When the domain is not intrinsically separable, fully separated representations can be performed by immersing of the non-separable domain onto a fully separable one and use penalty formulations [15,18]. Another valuable route consists in using appropriate geometrical mappings [2,[15][16][17].
A step forward in the use of separated representations consisted in assuming model parameters as extra-coordinates leading to the definition of multi-parametric solutions of parametrized models, that once calculated offline, enable simulation, optimization, inverse analysis, uncertainty propagation and simulation-based control, all them under the real-time constraint [8,10].
An appealing use of those separated representations concern their use in shape optimization as reported in [1]. However, several difficulties were noticed, among them the necessity of easily parametrizing complex geometries in a simple and compact form when addressing CAD geometries. There are techniques in CAD enabling accurate geometry descriptions based on the employ of NURBS. Recently such geometrical descriptors were combined with a discretization operating on the same framework, leading to the so-called Isogeometric Analysis (IGA) [12,19]. Thus, the IGA unifies the NURBS-based geometry description and a NURBS-based field approximation. Some advantages of the IGA are [19]: elimination of mesh generation process, exact (or very accurate) geometric representation, easy formulation of high order element and superior numerical accuracy than standard FEM. Several works have shown superior accuracy and robustness of IGA and excellent behavior has been proven in a number of publications [4,11,14,[19][20][21][22]28]. The NURBS can be also used in conjunction with usual discretization. Thus, they were also applied for enhancing standard finite element discretization [34]. Beside the above-mentioned methods, hyper reduction [24] and reduced basis isogeometric analysis [23] were also considered for geometry parametrization within the model order reduction setting.
The present paper aims at exploring the use of NURBS-based geometrical descriptions combined with the PGD-based fully 3D space separation. Moreover, considering the control points as model extra-coordinates within the PGD rationale we will address parametric geometries, that could be enriched by considering other parameters as extracoordinates, like material parameters.
All these points will be addressed in the next sections, before illustrating the solution capabilities from the solution and analysis of few selected problems. For the sake of simplicity, and without loss of accuracy, the different capabilities will be illustrated on the heat equation.

Problem statement
As stated before, the main goal of the present work is to develop a procedure for solving boundary value problems (BVP) in non-separable geometries using the PGD. Therefore, to describe the proposed method and also to evaluate it, we are considering the nonlinear heat equation with temperature dependent thermal conductivity. The heat equation defined in the domain x reads where T is the temperature field, K (T ) is the thermal conductivity assumed depending on the temperature and here assumed isotropic, and f is the volumetric heat source. Appropriate boundary conditions must be enforced on the domain boundary ≡ ∂ x and with the heat flux q = −K ∇T · n (with ∇ the gradient operator and n the unit outwards vector defined on the domain boundary), and T and q the prescribed temperature and heat flux on the Dirichlet D and Neumann N boundaries respectively.
The Galerkin weak formulation of the heat equation reads where T * represents the test function, assumed regular enough and vanishing in D .

NURBS-based geometrical approximation
Using the NURBS basis functions (introduced in Appendix A), a single variate curve C, a bivariate surface S and a trivariate volume V can be described according to and respectively, where, P refers to the control points (the vertices of the so-called control net) in the physical domain (curve, surface or volume). Equations These expressions allow computing the different partial derivatives involved in the differential operators as well as in the transformation Jacobian. The main drawback of such a mapping is that those derivatives have not an affine structure. For that purpose we apply a separated approximate representation -SAR-described in the next section.

Separated approximate representation (SAR)
Consider the generic function g(x) defined in the physical domain x and g(ξ), its counterpart in the computational domain ξ . A separated approximate representation of g(ξ) in the computational domain reads where M j (ξ j ) is the vector of approximation functions (such as usual piecewise linear function like the ones employed in FEM) in term of j-th coordinate direction (ξ 1 , ξ 2 , ξ 3 ≡ ξ , η, ζ ) and G is the vector of coefficients. ND is the number of coordinates (2 in 2D and 3 in 3D) and NG is the number of modes used to express function g(ξ). The procedure for calculating the unknown coefficients G ji was addressed in [9], here revisited in Appendix B.
The SAR is used for performing the affine decomposition of the terms involved in the geometrical mapping and where, NX, NY and NZ are the number of modes required to reach the desired level of accuracy in the SAR of x(ξ), y(ξ) and z (ξ) respectively. The associated Jacobian expression is derived in Appendix .

Proper generalized decomposition constructor
The separated representation of the temperature field T (ξ) in the computational domain ξ is constructed by progressive enrichment. Thus, with (n-1) modes already calculated, the n term approximate reads Finding the coefficient vectors T jn , consists of using a fixed-point iteration [9]. At each iteration of the fixed-point algorithm, one vector among the ND unknown vectors T jn , is considered unknown and the others are considered known.
Therefore, the test field involved in the weighted residual formulation, reads where d represents the considered direction, for calculating T dn . The other source of nonlinearity is the thermal conductivity that is assumed depending on the temperature. The simplest approach to overcome this issue, among many other alternatives, consists in evaluating it at the previous temperature approximate, and then apply the SAR according to When considering three dimensions, i.e. ND = 3, the problem weak form reads The weak form is then mapped into the computational domain ξ where the separated representation of the temperature field will apply, as well as the one of the conductivity.
For this purpose we will use the SAR of the transformation derivatives and the Jacobian determinant.
Thus, for example the first term in Eq. (20) leads to with K (ξ) given by the SAR (19), and where h 11 , h 12 and h 13 are the coefficients relating the derivatives in the physical and computational domains, with |J | the determinant of the Jacobian of the transformation (the expression of all them are given in Appendix ). Equation (21) results at its turn in nine terms, where the first reads The 27 resulting terms involved in the left-hand side of Eq. (20) have the same structure, expressed by where indices a, b, c take values in the set {1, 2, 3}. The solution procedure is detailed in Appendix D.

Numerical results
In this section, to evaluate the applicability of the proposed method, some numerical examples are solved and the results are compared with exact solutions or the ones obtained by using the FEM that will be considered as reference.
To evaluate the accuracy the following relative error norm is used where T ref refers the reference solution and T PGD to the approximate numerical solution obtained using the proposed method.

2D geometry
The domain geometry and control net associated to a NURBS surface are shown in Fig. 1.
The order of the B-splines considered is p = 2 in ξ and q = 1 in η, associated to the knots (the interested reader can find the main concepts and notation related to the NURBS in Appendix A) A repeated knot is added in the knot vector κ ξ to enforce the sharp corner at the control point 8 in Fig. 1b. To exactly represent the circular part of the boundary we consider the weighting coefficients [30] w = [1, cos(π/8), 1, cos(π/8), 1, 1, 1, 1, 1, 1] T .
A four-modes SAR of the domain geometry is shown in Fig. 2. Then, by considering the boundary conditions on the different boundary segments reported in Fig. 1a T (x, y) = and the temperature dependency of the thermal conductivity expressed by the exact solution reads The forcing function f (x, y) can be obtained by substitution of the exact solution given by Eq. (29) into the problem model expressed by Eq. (6). The problem is then solved by using the proposed methodology whose first 6 modes are shown in Fig. 3.
The temperature distribution in the computational and physical domains is shown in Fig. 4. To evaluate the accuracy of the method and the effect the number of modes the problem is solved for different number of nodes and in each case the error norm with respect to the exact solution is computed using Eq. (24) and it is reported in Fig. 5.
The results show that the method converges very rapidly and the solution only requires few modes. As expected the solution accuracy increases by considering finer meshes, with no impact on the computational performances because of the use of separated representation that only solves one-dimensional problems.

2D complex geometry
In the second example, we consider the geometry given in Fig. 6, whose left image represents a 2D NURBS surface and its control net. When considering second order B-splines along the direction ξ (p = 2) and one along direction η (q = 1), and using the knots and the SAR of coordinate x(ξ) and y(ξ) related to a single mode, i.e. NX = NY = 1, results in the reconstructed geometry depicted in Fig. 6(centre). As it can be noticed in Fig. 6(right) two SAR modes, NX = NY = 2, suffice for generating an almost perfect representation of the mapping. The nonlinear thermal problem is solved in the domain shown in Fig. 6 for proving the ability to perform separated representations in geometrically complex domains. Homogeneous Dirichlet boundary conditions are enforced on the domain boundary and the temperature dependent conductivity is taken the same that in the previous numerical test. A volumetric thermal source applies in the domain, f = 300. A mesh containing 61 nodes in ξ and 31 in η where considered. Even if a standard mesh-based discretization  would imply 61 × 31 nodes (even more if operating directly in the physical domain for accurately adapting to the domain geometry), the separated representation only needs solving a sequence of 1D problems with 61 degrees of freedom or 31 depending on the direction considered, ξ and η respectively. Again 6 modes (depicted in Fig. 7) were enough for representing accurately the solution, when compared with a standard FEM solution on a mesh consisting of quadratic elements and involving 4767 nodes. The converged temperature field is shown in Fig. 8 in the computational and physical domains, with the FEM reference solution. The maximum relative error concerning the maximum temperature is of 0.24% that can be considered excellent, and it could be reduced even more by considering either finer meshes or richer separated approximations.

3D geometry
In the present case, the thermal problem is solved in the 3D domain generated by rotating the 2D surface depicted in Fig. 9a around the z axis. The generated 3D volume is shown in Fig. 9b. That volume can be easily represented by using NURBS with 12 control points. The control net is shown in Fig. 9c. The order of the B-splines are p = 2 in ξ , q = 1 in η and r = 1 in ζ with the knot vectors given by The weighting coefficients corresponding to control points 2, 5, 8 and 11 are cos( π 4 ) and the weighting coefficients of the rest of control points are unit values. This selection of weights ensures the exact representation of the conic surface [30].
The modes of the separated representation of the physical coordinates x, y and z in terms of the computational coordinates ξ , η and ζ are shown in Fig. 10. To represent the domain geometry, only two modes for x and y and one mode for z suffice.
Homogeneous Dirichlet boundary conditions are applied on the domain boundary. The thermal conductivity is again the same than the one considered in the previous examples. The exact temperature distribution was enforced to be needing for a volumetric thermal source that results from introducing the temperature field (33) into the problem strong form (6). The solution involves the six modes shown in Fig. 11. The solution is depicted on some domain cross-sections of the computational domain in Fig. 12a whereas the corresponding solution in the physical domain is shown in Fig. 12b. The problem is solved for different number of modes and nodes, with the relative error norm, obtained by using Eq. (24), reported in Fig. 13.

Parametric solutions
The present section addresses the case in which the PGD operates first, as previously discussed, for separating the space coordinates in the computational domain, and then for including the material conductivity as an extra-coordinate. For that purpose we consider the problem considered in Section with now the conductivity as extra-coordinate. Thus, we are interested in computing the temperature field at each point and for any value of the thermal conductivity (linear, isotropic and uniformly distributed in the physical domain), taking values in the interval K ∈ [1, 3].

Describing a solution exhibiting a steep gradient
To emphasize the performances of the proposed approach for describing at very low cost complex geometries in a very accurate way, in this section we consider the domain depicted in Fig. 17a, where the solution in each point is described by function f (x, y) f (x, y) = 2 π atan 3(x 2 + y 2 − 4) ,  Fig. 17b. As it can be noticed the solution exhibits a high gradient when the radius approaches r = 2. Figure 18 depicts different geometrical descriptions enabling finer meshes in the region where the solution exhibits higher gradients (r = 2). Figure 19 depicts the convergence in the solution representation when decreasing the mesh size in the three considered control nets. Finally Figs. 20 and 21 compare the solutions along the η coordinate (in the computational domain) and the radial one respectively, proving an excellent accuracy.

Parametrized shapes
In this section instead of using as extra-coordinate a parameter related to the physics itself (the conductivity in Section ), here we consider the control point positions as parameters, and then, within the PGD rationale, as extra-coordinates, for in this manner computing from one shot the solution for any possible geometry (shape) according to the considered parameterization.

Conclusions
PGD makes possible advanced parametric analyses, by considering parameters as model extra-coordinates, while circumventing the resulting curse of dimensionality by using separated representations. When applying separated representation to the space coordinates, the solution of 3D problems becomes a sequence of one-dimensional problems as soon as the 3D domain can be expressed in a separated form (hexahedra). This fact limits significantly the domain of applicability of space separated representations.
The present work enlarges the applicability of the PGD rationale to general non-regular geometries, expressed by using NURBS-based representations, very usual in CAD technologies. A separated representation is then proposed for the geometrical transformation, enabling the transformation of a 3D problem in a complex geometry into a sequence of one-dimensional problems solution. Many examples served to prove the capabilities and potential of the proposed methodology.
More complex geometries need multiple NURBS for approximating them. In particular, these complex geometries can be decomposed in simpler subdomains, each one approximated as described in the present paper, with some constraints applied on the NURBS, for taking into account the different transmission conditions on the interfaces. This topic is addressed in the second part of the present paper, that is being finalized.
The bivariate NURBS basis functions R pq ij (ξ , η) in 2D domains and trivariate NURBS basis functions R pqr ijk (ξ , η, ζ ) in 3D domains can be obtained using tensor product of univariate basis functions and respectively.
In the previous equations p, q, and r are the order of B-splines in directions ξ , η, and ζ , respectively and ξ = (ξ , η, ζ ) is the coordinates vector in the computational domain, and w ijk denote the weighting parameters. All the univariate basis functions are defined based on the B-splines defined in their respective coordinates. The knot vector κ η and κ ζ applying in the coordinates η and ζ read and respectively, being n, m and l the number of basis functions in directions ξ , η and ζ respectively. Using the NURBS basis functions just described, the single variate curve C, bivariate surface S and trivariate volume V can be described from and respective, where, P refers to the control points (the vertices of the so-called control net) in the physical domain. Equations (43)-(45) represent the application that maps any point ξ in the computational domain ξ to a point x in physical domain x . Figure 27 illustrates the just described construction, where the control net and control points are shown in the left, the physical domain x described using the NURBS-based mapping (45), in the central image, and the corresponding computational domain ξ , the unit cube shown in the right image. Now, using the notation x T = (x, y, z) ≡ (x 1 , x 2 , x 3 ) and ξ T = (ξ , η, ζ ) ≡ (ξ 1 , ξ 2 , ξ 3 ), for facilitating the mapping description compactness, the terms involved in the transformation of the differential operator from x to ξ read where P a ijk refers to the a-component of P ijk . The previous equation involves the NURBS derivatives. Thus, for example the derivative of where, W (ξ ) is the denominator of Eq. (38). The derivatives of bivariate and trivariate basis functions can also be obtained following the same rationale, and can be found in [30].
As it can be noticed, the main drawback of the proposed mapping is that expressions like (47) have not a separate structure, compromising the effectiveness of the PGD solver. For that purpose we apply a separated approximate representation described in B.

Appendix B: SAR constructor
Consider function g(ξ). It might be given analytically or it may be possible to numerically evaluate it at any point ξ. We are interested in obtaining a SAR of g(ξ) in terms of a set of approximation functions in different directions and a set of coefficients as given in Eq. (13). Any approximation basis could be considered in M to find the coefficients in G related to a finite sum involving NG terms providing an approximation accurate enough.
We proceed in an iterative manner, and we consider that we have already calculated the first (n − 1) modes and we would like to find the next mode, n.
The residual at the present iteration reads Therefore, the n-th mode must be obtained in such a way that makes an approximation for the residual, i.e.
To find the coefficients G jn that provides the best approximation of the residual r n−1 (ξ), the error is enforced to be minimum by prescribing null derivatives where, G α 1n , G β 2n and G γ 3n refer to the α-th, β-th and γ -th component of corresponding vectors.
Consider, for instance previously discussed, is the non-separable expression of the terms related to the Jacobian of the coordinate transformation. The derivatives transform according to ⎡ ⎢ ⎣ ∂(·)/∂x ∂(·)/∂y ∂(·)/∂z where, h is the inverse of the Jacobian, J . As previously indicated the coordinate transformation (46) leads to non separable expressions that compromise the PGD effectiveness. An alternative procedure consists in contructing a separated representation of the coordinates mapping and where, NX, NY and NZ are the number of modes required to reach the desired level of accuracy in the SAR of x(ξ), y(ξ) and z (ξ) respectively. Using the same rationale, we can derive the SAR of Jacobian determinant and the transformation derivatives, respectivaly To illustrate the SAR of NURBS we consider the case given in Fig. 6, whose left image represents a 2D NURBS surface and its control net. When considering second order Bsplines along the direction ξ (p = 2) and one along direction η (q = 1), and using the knots and κ η = [0, 0, 1, 1] T (62) the SAR of coordinate x(ξ) and y(ξ) are constructed by using the procedure just described, where a single mode representation, i.e. NX = NY = 1, results in the reconstructed geometry depicted in Fig. 6(centre). As it can be noticed in Fig. 6(right) two SAR modes NX = NY = 2 suffice for generating an almost perfect representation of the mapping. The different modes involved in the approximation are depicted in Fig. 28. The maximum error of SAR of x and y are given in Table 1 proving that only two modes largely suffices for capturing the geometric details.