Enhanced parametric shape descriptions in PGD-based space separated representations

Abstract

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 $\varvec{x}$ and time t, as well as on a series of parameters, grouped in vector ${\varvec{\mu }}$, related to the considered physics, the loading or the domain in which the problem is defined. Thus, the problem reads

$$\begin{aligned} {\mathcal {L}}(u(\varvec{x},t;{\varvec{\mu }})) = f(\varvec{x},t;{\varvec{\mu }}), \end{aligned}$$

(1)

where $\varvec{x} \in \Omega _{\varvec{x}}({\varvec{\mu }}) \subset {\mathbb {R}}^3$, $t \in \Omega _t \subset {\mathbb {R}}$ and the parameteres ${\varvec{\mu }}\in \Omega _\mu \subset {\mathbb {R}}^\mathtt 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 transient 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 $\phi _i(\varvec{x})$ for representing those solutions.

A generalization consists of expressing $u(\varvec{x},t)$ from the finite sum

$$\begin{aligned} u(\varvec{x},t) \approx \sum \limits _{i=1}^{\mathtt M} T_i(t) X_i(\varvec{x}), \end{aligned}$$

(2)

with both functions $T_i(t)$ and $X_i(\varvec{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 $\mathtt M$ problems in space for calculating functions $X_i(\varvec{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 $\Omega _{\varvec{x}}$ can be fully or partially decomposed, i.e. $\Omega _{\varvec{x}} = \Omega _x \times \Omega _y \times \Omega _z$ or $\Omega _{\varvec{x}} = \Omega _{xy} \times \Omega _z$ respectively, the separated representation based approximations read respectively

$$\begin{aligned} u(x,y,z) \approx \sum \limits _{i=1}^{\mathtt M} X_i(x) Y_i(y) Z_i(z), \end{aligned}$$

(3)

that allows calculating the 3D solution from a sequence of 1D problems, and

$$\begin{aligned} u(x,y,z) \approx \sum \limits _{i=1}^{\mathtt M} X_i(x,y) Z_i(z). \end{aligned}$$

(4)

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

$$\begin{aligned} u(\varvec{x},t,\mu _1, \dots , \mu _\mathtt P) \approx \sum \limits _{i=1}^{\mathtt M} X_i(\varvec{x}) T_i(t) \prod \limits _{j=1}^{\mathtt P} M_i^j(\mu _j), \end{aligned}$$

(5)

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 extra-coordinates, 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.

Geometrical description enabling fully separated representations

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 $\Omega _{\varvec{x}}$ reads

$$\begin{aligned} \frac{\partial }{\partial x}\left( K(T)\frac{\partial T}{\partial x}\right) +\frac{\partial }{\partial y}\left( K(T)\frac{\partial T}{\partial y}\right) +\frac{\partial }{\partial z}\left( K(T)\frac{\partial T}{\partial z}\right) +f=0 , \end{aligned}$$

(6)

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 $\Gamma \equiv \partial \Omega _{\varvec{x}}$

$$\begin{aligned} T={\overline{T}}\ \ \ \ \mathrm {on} \ \ \ \ \ \ {{\varGamma }}_D \end{aligned}$$

(7)

and

$$\begin{aligned} q={\overline{q}}\ \ \ \ \mathrm {on}\ \ \ \ \ \ {{\varGamma }}_N \end{aligned}$$

(8)

with the heat flux $q = - K \nabla T \cdot \varvec{n}$ (with $\nabla$ the gradient operator and $\varvec{n}$ the unit outwards vector defined on the domain boundary), and ${\overline{T}}$ and ${\overline{q}}$ the prescribed temperature and heat flux on the Dirichlet ${{\Gamma }}_D$ and Neumann ${{\Gamma }}_N$ boundaries respectively.

The Galerkin weak formulation of the heat equation reads

$$\begin{aligned}&\int _{{\varOmega _{\varvec{x}}}}{\left[ \frac{\partial T^*}{\partial x}\left( K\frac{\partial T}{\partial x}\right) +\frac{\partial T^*}{\partial y}\left( K\frac{\partial T}{\partial y}\right) +\frac{\partial T^*}{\partial z}\left( K\frac{\partial T}{\partial z}\right) \right] d\varvec{x}}\nonumber \\&\quad = \int _{{\varOmega _{\varvec{x}}}}{T^*fd\varvec{x}}+\int _{{{\varGamma }}_N}{T^*{\overline{q}}d\varvec{x}} , \end{aligned}$$

(9)

where $T^*$ represents the test function, assumed regular enough and vanishing in $\Gamma _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

$$\begin{aligned} {\varvec{x}}= & {} \sum ^n_{i=1}{R^p_i\left( \varvec{\xi }\right) {{\varvec{P}}}_i},\ \ \ \ \ \ \ \ \ \ {\varvec{x}}\in C, \end{aligned}$$

(10)

$$\begin{aligned} {\varvec{x}}= & {} \sum ^n_{i=1}{\sum ^m_{j=1}{R^{pq}_{ij}\left( \varvec{\xi }\right) {{\varvec{P}}}_{ij}}},\ \ \ \ \ \ \ \ \ \ {\varvec{x}}\in S, \end{aligned}$$

(11)

and

$$\begin{aligned} {\varvec{x}}=\sum ^n_{i=1}{\sum ^m_{j=1}{\sum ^l_{k=1}{R^{pqr}_{ijk}\left( \varvec{\xi }\right) {{\varvec{P}}}_{ijk}}}},\ \ \ \ \ \ \ \ \ \ {\varvec{x}}\in V , \end{aligned}$$

(12)

respectively, where, ${\varvec{P}}$ refers to the control points (the vertices of the so-called control net) in the physical domain (curve, surface or volume).

Equations (10)–(12) represent the application that maps any point ${\varvec{\xi }}$ into the computational domain $\Omega _{{\varvec{\xi }}}$ ([0, 1] in 1D, $[0,1]^2$ in 2D and $[0,1]^3$ in 3D) into a point $\varvec{x}$ in physical domain $\Omega _{\varvec{x}}$. The form of the NURBS-based approximation functions involved in these expressions is given in Appendix A.

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({\varvec{x}})$ defined in the physical domain $\Omega _{\varvec{x}}$ and $g({\varvec{\xi }})$, its counterpart in the computational domain $\Omega _{{\varvec{\xi }}}$. A separated approximate representation of $g(\varvec{\xi })$ in the computational domain reads

$$\begin{aligned} g(\varvec{\xi })\approx \sum ^{NG}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j({\xi }_j){{\varvec{G}}}_{ji}}} , \end{aligned}$$

(13)

where ${\varvec{M}}_j({\xi }_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 (${\xi }_1,{\xi }_2,{\xi }_3 \equiv \xi ,\eta ,\zeta$) and ${\varvec{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$(\varvec{\xi })$.

The procedure for calculating the unknown coefficients $\varvec{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

$$\begin{aligned} x(\varvec{\xi })= & {} \sum ^{NX}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{X}}}_{ji}}} , \end{aligned}$$

(14)

$$\begin{aligned} y(\varvec{\xi })= & {} \sum ^{NY}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{Y}}}_{ji}}} \end{aligned}$$

(15)

and

$$\begin{aligned} z(\varvec{\xi })=\sum ^{NZ}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{Z}}}_{ji}}} \end{aligned}$$

(16)

where, NX, NY and NZ are the number of modes required to reach the desired level of accuracy in the SAR of $x(\varvec{\xi })$, $y(\varvec{\xi })$ and $z\left( \varvec{\xi }\right)$ respectively. The associated Jacobian expression is derived in Appendix .

Proper generalized decomposition constructor

The separated representation of the temperature field $\mathrm {T}\left( \varvec{{\upxi }}\right)$ in the computational domain $\Omega _{{\varvec{\xi }}}$ is constructed by progressive enrichment. Thus, with (n-1) modes already calculated, the n term approximate reads

$$\begin{aligned} T\left( \varvec{\xi }\right) \approx \sum ^{n-1}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{T}}}_{ji}}}+\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{T}}}_{jn}} . \end{aligned}$$

(17)

Finding the coefficient vectors ${{\varvec{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 ${{\varvec{T}}}_{jn}$, is considered unknown and the others are considered known.

Therefore, the test field involved in the weighted residual formulation, reads

$$\begin{aligned} T^{*d}\left( \varvec{\xi }\right) ={{\varvec{M}}}_d^T \varvec{T}^{*d}\prod ^{ND}_{\small {\begin{array}{c} j=1 \\ j\ne d \end{array} }}{{{\varvec{M}}}^T_j{{\varvec{T}}}_{jn}} \end{aligned}$$

(18)

where d represents the considered direction, for calculating $\varvec{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

$$\begin{aligned} K\left( T\left( \varvec{\xi }\right) \right) =K\left( \varvec{\xi }\right) =\sum ^{NK}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{K}}}_{ji}}} . \end{aligned}$$

(19)

When considering three dimensions, i.e. $ND=3$, the problem weak form reads

$$\begin{aligned}&\int _{{\varOmega _{\varvec{x}}}}{\frac{\partial T^{*d}}{\partial x}\frac{\partial T}{\partial x}Kd\varvec{x}}+\int _{{\varOmega _{\varvec{x}}}}{\frac{\partial T^{*d}}{\partial y}\frac{\partial T}{\partial y}Kd\varvec{x}}+\int _{{\varOmega _{\varvec{x}}}}{\frac{\partial T^{*d}}{\partial z}\frac{\partial T}{\partial z}Kd\varvec{x}}\nonumber \\&\quad = \int _{{\varOmega _{\varvec{x}}}}{T^{*d}fd\varvec{x}}+\int _{{{\varGamma }}_N}{T^{*d}{\overline{q}}d\varvec{x}} . \end{aligned}$$

(20)

The weak form is then mapped into the computational domain $\Omega _{{\varvec{\xi }}}$ 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

$$\begin{aligned}&\int _{\Omega _{{\varvec{\xi }}} }{\left( h_{11}\frac{\partial T^{*d}}{\partial \xi }+h_{12}\frac{\partial T^{*d}}{\partial \eta }+h_{13}\frac{\partial T^{*d}}{\partial \zeta }\right) }\nonumber \\&\quad {\left( h_{11}\frac{\partial T}{\partial \xi }+h_{12}\frac{\partial T}{\partial \eta }+h_{13}\frac{\partial T}{\partial \zeta }\right) K({\varvec{\xi }}) \left| {{\varvec{J}}}\right| d{\varvec{\xi }}} , \end{aligned}$$

(21)

with $K({\varvec{\xi }})$ 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 $\left| {{\varvec{J}}}\right|$ 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

$$\begin{aligned} \int _{\Omega _{{\varvec{\xi }}} }{\frac{\partial T^{*d}}{\partial \xi }\frac{\partial T}{\partial \xi }h_{11}h_{11}K({\varvec{\xi }})\left| {{\varvec{J}}}\right| d{\varvec{\xi }}} . \end{aligned}$$

(22)

The 27 resulting terms involved in the left-hand side of Eq. (20) have the same structure, expressed by

$$\begin{aligned} \int _{\Omega _{{\varvec{\xi }}} }{\frac{\partial T^{*d}}{\partial {\xi }_a}\frac{\partial T}{\partial {\xi }_b}h_{ca}h_{cb}K({\varvec{\xi }})\left| {{\varvec{J}}}\right| d{\varvec{\xi }}} , \end{aligned}$$

(23)

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

$$\begin{aligned} e= \left( \frac{\int _{\Omega _{\varvec{x}}} (T^{ref}-T^{PGD})^2 \ d \varvec{x}}{ \int _{\Omega _{\varvec{x}}} ( T^{ref} )^2 \ d \varvec{x} } \right) ^{1/2} , \end{aligned}$$

(24)

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 $\xi$ and $q=1$ in $\eta$, associated to the knots (the interested reader can find the main concepts and notation related to the NURBS in Appendix A)

$$\begin{aligned} \left\{ \begin{array}{l} {\varvec{\kappa }}_{\xi }={[0,0,0,0.5,0.5,1,1,1]}^T \\ {\varvec{\kappa }}_{\eta }={[0,0,1,1]}^T \end{array} \right. . \end{aligned}$$

(25)

Fig. 1

Domain geometry (a) and control net for the construction of the NURBS representation (b)

A repeated knot is added in the knot vector ${\varvec{\kappa }}_{\xi }$ 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]

$$\begin{aligned} w={[1,cos(\pi /8),1,cos(\pi /8),1,1,1,1,1,1]}^T . \end{aligned}$$

(26)

A four-modes SAR of the domain geometry is shown in Fig. 2.

Fig. 2

Four-modes ($i=1,2,3,4$) SAR of the NURBS-based mapping: a $\varvec{M}^T_1(\xi ) \varvec{X}_{1i}$; b $\varvec{M}^T_2(\eta ) \varvec{X}_{2i}$; c $\varvec{M}^T_1(\xi ) \varvec{Y}_{1i}$; and d $\varvec{M}^T_2(\eta ) \varvec{Y}_{2i}$

Then, by considering the boundary conditions on the different boundary segments reported in Fig. 1a

$$\begin{aligned} T(x,y)=\left\{ \begin{array}{ll} 0 &{} {\ \mathrm{on} \ AB} \\ (4-2x)(x^2-1) &{} {\ \mathrm{on} \ BC} \\ 0 &{} {\ \mathrm{on} \ CDE} \\ (4-2y)(y^2-1) &{} { \mathrm{on} \ AE} \end{array} \right. , \end{aligned}$$

(27)

and the temperature dependency of the thermal conductivity expressed by

$$\begin{aligned} K(T)=1+T , \end{aligned}$$

(28)

the exact solution reads

$$\begin{aligned} T\left( x,y\right) =\left( x-2\right) \left( y-2\right) (x^2+y^2-1) . \end{aligned}$$

(29)

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.

Fig. 3

Six-modes separated representation of the temperature field: Modes in $\xi$ (a) and in $\eta$ (b)

Fig. 4

Temperature field in the computational domain (a) and in the physical domain (b)

Fig. 5

Relative error norm versus the number of modes involved in the temperature separated representation for different meshes

Fig. 6

NURBS surface in 2D and its control net (a), single mode SAR (b); and two modes SAR (c)

Fig. 7

Modes involved in the solution: Modes in $\xi$ (a) and modes in $\eta$ (b)

Fig. 8

Temperature field: In the computational domain (a); In physical domain (b); and reference solution obtained by using a commercial FEM solver

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.

Fig. 9

Domain geometry: Dimensions of the generative surface (a); the resulting volume (b); and the control net for the NURBS representation (c)

Fig. 10

SAR of the NURBS-based mapping: a $\varvec{M}^T_1(\xi ) \varvec{X}_{1i}$; b $\varvec{M}^T_2(\eta ) \varvec{X}_{2i}$; c $\varvec{M}^T_3(\zeta ) \varvec{X}_{3i}$; d $\varvec{M}^T_1(\xi ) \varvec{Y}_{1i}$; e $\varvec{M}^T_2(\eta ) \varvec{Y}_{2i}$; f $\varvec{M}^T_3(\zeta ) \varvec{Y}_{3i}$; g $\varvec{M}^T_1(\xi ) \varvec{Z}_{1i}$; h $\varvec{M}^T_2(\eta ) \varvec{Z}_{2i}$; i $\varvec{M}^T_3(\zeta ) \varvec{Z}_{3i}$

Fig. 11

Six modes involved in the solution: a Modes in $\xi$; b modes in $\eta$; and c modes in $\zeta$

Fig. 12

Temperature field: a Computational domain; b Physical domain

Fig. 13

Relative error norm versus the number of modes for different meshes

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 $\xi$ ($p=2$) and one along direction $\eta$ ($q=1$), and using the knots

$$\begin{aligned} {\varvec{\kappa }}_{\xi }={[0,0,0,0.25,0.4,0.5,0.6,0.75,1,1,1]}^T , \end{aligned}$$

(30)

and

$$\begin{aligned} {\varvec{\kappa }}_{\eta }={[0,0,1,1]}^T . \end{aligned}$$

(31)

the SAR of coordinate $x(\varvec{\xi })$ and $y(\varvec{\xi })$ 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 $\xi$ and 31 in $\eta$ where considered. Even if a standard mesh-based discretization would imply $61 \times 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, $\xi$ and $\eta$ 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.

Fig. 14

Modes involved in the separated representation: a modes involving the coordinate $\xi$; b Modes involving the coordinate $\eta$; and c Modes involving the thermal conductivity K

Fig. 15

Temperature distribution for different conductivities in the computational domain

Fig. 16

Temperature distribution for different conductivities in the physical domain

Fig. 17

Domain and function contours (a) and 3D representation of the considered function (b)

Fig. 18

Different control nets to represent the same geometry but with different stretching to capture accurately steep gradients (a, c, e) and (b, d, f) corresponding NURBS geometry

Fig. 19

Convergence analysis

Fig. 20

Function distribution along the coordinate $\eta$

Fig. 21

Solution f(r), with r the radius

Fig. 22

Geometry parameterization of a complex Y shaped domain: extreme positions of the moving control point

Fig. 23

Temperature field in the computational domains particularized for different positions S of the control point indicated in Fig. 22, that is, for different shapes of the physical domain

Fig. 24

Solutions in the physical domain related to those shown in Fig. 23

Fig. 25

Modes involved in the separated representation: a modes involving the coordinate $\xi$; b Modes involving the coordinate $\eta$; and c Modes involving the control point position S

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 $\xi$, $q=1$ in $\eta$ and $r=1$ in $\zeta$ with the knot vectors given by

$$\begin{aligned} \left\{ \begin{array}{l} {\varvec{\kappa }}_{\xi }={[0,0,0,1,1,1]}^T \\ {\varvec{\kappa }}_{\eta }={[0,0,1,1]}^T \\ {\varvec{\kappa }}_{\zeta }={[0,0,1,1]}^T \end{array} \right. . \end{aligned}$$

(32)

The weighting coefficients corresponding to control points 2, 5, 8 and 11 are $\mathrm {cos}\mathrm {}(\frac{\pi }{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 $\xi$, $\eta$ and $\zeta$ 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

$$\begin{aligned} T\left( x,y,z\right) =\frac{1}{50}xyz(z-2)(x^2+y^2-1)\left( x^2+y^2-{ \left( 3-\frac{1}{2}z \right) }^2\right) , \end{aligned}$$

(33)

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 \in [1,3]$.

Figure 14 represents the different modes involved in the solution separated representation, whereas Figs. 15 and 16 show the solution in both, the computational and the physical domains respectively.

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)

$$\begin{aligned} f(x,y) = \frac{2}{\pi } \mathrm {atan} \left( 3 (x^2+y^2 -4)\right) , \end{aligned}$$

(34)

depicted in 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 $\eta$ 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.

For that purpose, we consider the Y-shaped domain already introduced in section , where the location of one control point is considered as a parameter, that is, as a model extra-coordinate within the PGD rationale. The extreme positions of that control point and the corresponding geometries, are shown in Fig. 22. The temperature distribution in the computational and physical domains are given respectively in Figs. 23 and 24 on different slices corresponding to different values of the coordinate S (representing the control point position). Finally, the first 5 modes of the separated representation, concerning the different coordinates, are depicted in Fig. 25.

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.

Appendices

Appendix A: NURBS fundamentals

NURBS enable exactly representing the geometry of most curves, surfaces and volumes encountered in engineering applications, being nowadays a usual technology employed in CAD.

For the sake of completeness this section revisits the main concepts involved in the NURBS construction. For additional details the interested reader can refer to [30].

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 $\xi \in [0,1]$, the knot vector ${\varvec{\kappa }}_{\xi }$ is defined as follows

$$\begin{aligned} \varvec{\kappa }_{\xi }={[{\xi }_1,\ {\xi }_2,\ \dots ,{\xi }_{n+p+1}]}^T. \end{aligned}$$

(35)

Denote $N_{ap}(\xi )$ as the B-spline basis function of order p in the a-th knot span $\xi \in [{\ }{\xi }_a,{\ }{\xi }_{a+1})$. The following recursive equations can be used to compute the univariate B-spline basis function $N_{ap}(\xi )$ [30]

$$\begin{aligned} N_{a0}(\xi )=\left\{ \begin{array}{ll} 1 &{} {\xi }_a\le \xi <{\xi }_{a+1} \\ 0 &{} {\mathrm {otherwise} } \end{array} \right. , \end{aligned}$$

(36)

for $p=0$, and for $p>0$

$$\begin{aligned} N_{ap}\left( \xi \right) =\frac{\xi -{\xi }_a}{{\xi }_{a+p}-{\xi }_a}N_{a\left( p-1\right) }+\frac{{\xi }_{a+p+1}-\xi }{{\xi }_{a+p+1}-{\xi }_{a+1}}N_{(a+1)\left( p-1\right) } \end{aligned}$$

(37)

For illustrative purposes a set of basis functions of order 2 using a non-uniform knot vector ${\varvec{\kappa }}_{\xi }\mathrm {=}{\mathrm {[0,0,0,0.3,0.4,0.5,0.6,0.7,0.7,1,1,1]}}^T$ is shown in Fig. 26.

Fig. 26

The B-spline basis functions of order 2 for non-uniform knot vector $\varvec{\kappa }={[\mathrm {0,0,0,0.3,0.4,0.5,0.6,0.7,0.7,1,1,1}]}^T$

A proper choice of the knot vectors allows obtaining rich behavior of the basis functions and enough flexibility to describe complex geometries. For example, the basis functions given in Fig. 26 define an interpolation at the ends of the interval [0, 1] and also at the position of $\xi =0.7$ due to the repeated knot at these locations.

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(\xi )$ reads [30]

$$\begin{aligned} R^p_i\left( \xi \right) =\frac{N_{ip}\left( \xi \right) w_i}{\sum ^n_{\alpha =1}{N_{\alpha p}\left( \xi \right) w_{\alpha }}} . \end{aligned}$$

(38)

The bivariate NURBS basis functions $R^{pq}_{ij}\left( \xi ,\eta \right)$ in 2D domains and trivariate NURBS basis functions $R^{pqr}_{ijk}\left( \xi ,\eta ,\zeta \right)$ in 3D domains can be obtained using tensor product of univariate basis functions

$$\begin{aligned} R^{pq}_{ij}\left( \varvec{\xi }\right) =\frac{N_{ip}\left( \xi \right) N_{jq}\left( \eta \right) w_{ij}}{\sum ^n_{\alpha =1}{\sum ^m_{\beta =1}{N_{\alpha p}\left( \xi \right) N_{\beta q}\left( \eta \right) w_{\alpha \beta }}}} , \end{aligned}$$

(39)

and

$$\begin{aligned} R^{pqr}_{ijk}\left( \varvec{\xi }\right) =\frac{N_{ip}\left( \xi \right) N_{jq}\left( \eta \right) N_{kr}\left( \zeta \right) w_{ijk}}{\sum ^n_{\alpha =1}{\sum ^m_{\beta =1}{\sum ^l_{\gamma =1}{N_{\alpha p}\left( \xi \right) N_{\beta q}\left( \eta \right) N_{\gamma r}\left( \zeta \right) w_{\alpha \beta \gamma }}}}} , \end{aligned}$$

(40)

respectively.

In the previous equations p, q, and r are the order of B-splines in directions $\xi$, $\eta$, and $\zeta$, respectively and $\varvec{\xi }\varvec{=}\left( \xi ,\eta ,\zeta \right)$ 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 ${\varvec{\kappa }}_{\eta }$ and ${\varvec{\kappa }}_{\zeta }$ applying in the coordinates $\eta$ and $\zeta$ read

$$\begin{aligned} {\varvec{\kappa }}_{\eta }={[{\eta }_1,\ {\eta }_2,\ \dots ,{\eta }_{m+q+1}]}^T\ \ \ \ \ \ \ \ \ \eta \in [0,1] , \end{aligned}$$

(41)

and

$$\begin{aligned} {\varvec{\kappa }}_{\zeta }={[{\zeta }_1,\ {\zeta }_2,\ \dots ,{\zeta }_{l+r+1}]}^T\ \ \ \ \ \ \ \ \ \zeta \in [0,1] , \end{aligned}$$

(42)

respectively, being n, m and l the number of basis functions in directions $\xi$, $\eta$ and $\zeta$ respectively.

Using the NURBS basis functions just described, the single variate curve C, bivariate surface S and trivariate volume V can be described from

$$\begin{aligned} {\varvec{x}}= & {} \sum ^n_{i=1}{R^p_i\left( \varvec{\xi }\right) {{\varvec{P}}}_i}\ \ \ \ \ \ \ \ \ \ {\varvec{x}}\in C, \end{aligned}$$

(43)

$$\begin{aligned} {\varvec{x}}= & {} \sum ^n_{i=1}{\sum ^m_{j=1}{R^{pq}_{ij}\left( \varvec{\xi }\right) {{\varvec{P}}}_{ij}}}\ \ \ \ \ \ \ \ \ \ {\varvec{x}}\in S, \end{aligned}$$

(44)

and

$$\begin{aligned} {\varvec{x}}=\sum ^n_{i=1}{\sum ^m_{j=1}{\sum ^l_{k=1}{R^{pqr}_{ijk}\left( \varvec{\xi }\right) {{\varvec{P}}}_{ijk}}}}\ \ \ \ \ \ \ \ \ \ {\varvec{x}}\in V , \end{aligned}$$

(45)

respective, where, ${\varvec{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 ${\varvec{\xi }}$ in the computational domain $\Omega _{{\varvec{\xi }}}$ to a point $\varvec{x}$ in physical domain $\Omega _{\varvec{x}}$.

Figure 27 illustrates the just described construction, where the control net and control points are shown in the left, the physical domain $\Omega _{\varvec{x}}$ described using the NURBS-based mapping (45), in the central image, and the corresponding computational domain $\Omega _{{\varvec{\xi }}}$, the unit cube shown in the right image.

Fig. 27

The control net for generating a NURBS volume (a), NURBS-based description of the physical domain (b) and computational domain (c)

Now, using the notation $\varvec{x}^T = (x,y,z) \equiv (x_1,x_2, x_3)$ and ${\varvec{\xi }}^T = (\xi ,\eta ,\zeta ) \equiv (\xi _1,\xi _2, \xi _3)$, for facilitating the mapping description compactness, the terms involved in the transformation of the differential operator from $\Omega _{\varvec{x}}$ to $\Omega _{{\varvec{\xi }}}$ read

$$\begin{aligned} \frac{\partial x_a}{\partial {\xi }_b}=\sum ^n_{i=1}{\sum ^m_{j=1}{\sum ^l_{k=1}{\frac{\partial \ R^{pqr}_{ijk}\left( \varvec{\xi }\right) }{\partial \ {\xi }_b}\ \ P^a_{ijk}}}},\ \ \ \ \ \ \ a,b=1,2,3 , \end{aligned}$$

(46)

where $P^a_{ijk}$ refers to the a-component of $\varvec{P}_{ijk}$.

The previous equation involves the NURBS derivatives. Thus, for example the derivative of $R_i^p(\xi )$ reads

$$\begin{aligned} \frac{d}{d\xi } R^p_i (\xi ) = w_i \frac{ N_{ip}' (\xi ) W (\xi ) - N_{ip} (\xi ) W'(\xi ) }{ { W(\xi ) }^2 } \end{aligned}$$

(47)

where, $W\left( \xi \right)$ 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(\varvec{\xi })$. It might be given analytically or it may be possible to numerically evaluate it at any point $\varvec{\xi }$. We are interested in obtaining a SAR of $g(\varvec{\xi })$ 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 ${\varvec{M}}$ to find the coefficients in ${\varvec{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

$$\begin{aligned} r^{n-1}\left( \varvec{\xi }\right) =g\left( \varvec{\xi }\right) \varvec{-}\sum ^{n-1}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j(\xi _j) {{\varvec{G}}}_{ji}}} . \end{aligned}$$

(48)

Therefore, the n-th mode must be obtained in such a way that makes an approximation for the residual, i.e.

$$\begin{aligned} \prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{G}}}_{jn}}\approx r^{n-1}\left( \varvec{\xi }\right) . \end{aligned}$$

(49)

To find the coefficients ${{\varvec{G}}}_{jn}$ that provides the best approximation of the residual $r^{n-1}\left( \varvec{\xi }\right)$, the error

$$\begin{aligned} \left\| e\right\| =\int _{\Omega _{{\varvec{\xi }}} }{{\left( \prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{G}}}_{jn}}-r^{n-1}\left( \varvec{\xi }\right) \right) }^2d {\varvec{\xi }}} , \end{aligned}$$

(50)

is enforced to be minimum by prescribing null derivatives

$$ \begin{aligned} \frac{\partial \left\| e\right\| }{\partial G^{\alpha }_{1n}}=0\ \ ,\ \ \frac{\partial \left\| e\right\| }{\partial G^{\beta }_{2n}}=0\ \ \& \ \ \frac{\partial \left\| e\right\| }{\partial G^{\gamma }_{3n}}=0, \end{aligned}$$

(51)

where, $G^{\alpha }_{1n}$ , $G^{\beta }_{2n}$ and $G^{\gamma }_{3n}$ refer to the $\alpha$-th, $\beta$-th and $\gamma$-th component of corresponding vectors.

Consider, for instance

$$\begin{aligned} \frac{\partial \left\| e\right\| }{\partial G^{\alpha }_{1n}}=0, \end{aligned}$$

(52)

it results

$$\begin{aligned}&\int _{\Omega _{{\varvec{\xi }}} }{M^{\alpha }_1\left( {{\varvec{M}}}^T_1{{\varvec{G}}}_{1n}\right) {\left( {{\varvec{M}}}^T_2{{\varvec{G}}}_{2n}\right) }^2{\left( {{\varvec{M}}}^T_3{{\varvec{G}}}_{3n}\right) }^2d{\varvec{\xi }}} \nonumber \\&\quad = \int _{\Omega _{{\varvec{\xi }}} }{M^{\alpha }_1\left( {{\varvec{M}}}^T_2{{\varvec{G}}}_{2n}\right) \left( {{\varvec{M}}}^T_3{{\varvec{G}}}_{3n}\right) r^{n-1}(\varvec{\xi })d{\varvec{\xi }}} , \end{aligned}$$

(53)

that reduces to

$$\begin{aligned}&{\left( \int ^1_0{M^{\alpha }_1{{\varvec{M}}}_1d\xi }\right) }^T{{\varvec{G}}}_{1n} \nonumber \\&\quad = \frac{\int ^1_0{\int ^1_0{\int ^1_0{M^{\alpha }_1\left( {{\varvec{M}}}^T_2{{\varvec{G}}}_{2n}\right) \left( {{\varvec{M}}}^T_3{{\varvec{G}}}_{3n}\right) \ r^{n-1}\left( \varvec{\xi }\right) \ d\zeta d\eta d\xi }}}}{\int ^1_0{{\left( {{\varvec{M}}}^T_2{{\varvec{G}}}_{2n}\right) }^2d\eta }\int ^1_0{{\left( {{\varvec{M}}}^T_3{{\varvec{G}}}_{3n}\right) }^2d\zeta }} , \ \forall \alpha \end{aligned}$$

(54)

By considering the other two relations in Eq. (51), each one applying $\forall \beta$ and $\forall \gamma$ respectively, a nonlinear discrete system is obtained whose solution result in the vectors $\varvec{G}_{1n}$, $\varvec{G}_{2n}$ and $\varvec{G}_{3n}$.

Appendix C: Separated geometry mapping

As before stated in Section , the most appealing features of the PGD is reducing of a high dimensional problem to a series of problems defined in lower dimensional spaces.

To apply this method in non-separable geometries one needs first to map the problem from the physical domain $\Omega _{\varvec{x}}$ into the computational domain $\Omega _{{\varvec{\xi }}}$. The main issue, as previously discussed, is the non-separable expression of the terms related to the Jacobian of the coordinate transformation.

The derivatives transform according to

$$\begin{aligned} \left[ \begin{array}{c} {\partial (\cdot )}/{\partial x} \\ {\partial (\cdot )}/{\partial y} \\ {\partial (\cdot )}/{\partial z} \end{array} \right] ={\varvec{h}}\left[ \begin{array}{c} {\partial (\cdot )}/{\partial \xi } \\ {\partial (\cdot )}/{\partial \eta } \\ {\partial (\cdot )}/{\partial \zeta } \end{array} \right] , \end{aligned}$$

(55)

where, ${\varvec{h}}$ is the inverse of the Jacobian, ${\varvec{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

$$\begin{aligned} x(\varvec{\xi })= & {} \sum ^{NX}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{X}}}_{ji}}} , \end{aligned}$$

(56)

$$\begin{aligned} y(\varvec{\xi })= & {} \sum ^{NY}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{Y}}}_{ji}}} \end{aligned}$$

(57)

and

$$\begin{aligned} z(\varvec{\xi })=\sum ^{NZ}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{Z}}}_{ji}}} \end{aligned}$$

(58)

where, NX, NY and NZ are the number of modes required to reach the desired level of accuracy in the SAR of $x(\varvec{\xi })$, $y(\varvec{\xi })$ and $z\left( \varvec{\xi }\right)$ respectively.

Using the same rationale, we can derive the SAR of Jacobian determinant and the transformation derivatives, respectivaly

$$\begin{aligned} \left| {\varvec{J}} \right| =\sum ^{NJ}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{ {\varvec{J}}}_{ji} }} , \end{aligned}$$

(59)

and

$$\begin{aligned} h_{ab}\left( \varvec{\xi }\right) =\sum ^{NH^{ab}}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{H}}}^{ab}_{ji}}},\ \ \ \ a,b=1,2,3 . \end{aligned}$$

(60)

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 B-splines along the direction $\xi$ ($p=2$) and one along direction $\eta$ ($q=1$), and using the knots

$$\begin{aligned} {\varvec{\kappa }}_{\xi }={[0,0,0,0.25,0.4,0.5,0.6,0.75,1,1,1]}^T , \end{aligned}$$

(61)

and

$$\begin{aligned} {\varvec{\kappa }}_{\eta }={[0,0,1,1]}^T \end{aligned}$$

(62)

the SAR of coordinate $x(\varvec{\xi })$ and $y(\varvec{\xi })$ 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.

Fig. 28

Two-modes ($i=1,2$) SAR of the NURBS-based mapping: a $\varvec{M}^T_1(\xi ) \varvec{X}_{1i}$; b $\varvec{M}^T_2(\eta ) \varvec{X}_{2i}$; c $\varvec{M}^T_1(\xi ) \varvec{Y}_{1i}$; and d $\varvec{M}^T_2(\eta ) \varvec{Y}_{2i}$

Table 1 Maximum error in the separated approximate representation of coordinates x and y

Appendix D: PGD constructor

First, we define the operator

$$\begin{aligned} \overbrace{{{\varvec{M}}}_j}^{a}=\left\{ \begin{array}{cc} \frac{\partial {{\varvec{M}}}_a}{\partial {\xi }_a} &{} \mathrm {if}\ \ a=j \\ {{\varvec{M}}}_j &{} \mathrm {else} \end{array} \right. . \end{aligned}$$

(63)

By considering the above operator the derivatives of T and $T^{*d}$ read

$$\begin{aligned}&\frac{\partial T}{\partial {\xi }_b}=\frac{\partial }{\partial {\xi }_b}\left( \sum ^{n-1}_{i=1}{\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{T}}}_{ji}}}+\prod ^{ND}_{j=1}{{{\varvec{M}}}^T_j{{\varvec{T}}}_{jn}}\right) \nonumber \\&\quad = \sum ^{n-1}_{i=1}{\prod ^{ND}_{j=1}{\overbrace{{{\varvec{M}}}^T_j}^{b}{{\varvec{T}}}_{ji}}}+\prod ^{ND}_{j=1}{\overbrace{{{\varvec{M}}}^T_j}^{b}{{\varvec{T}}}_{jn}} \end{aligned}$$

(64)

and

$$\begin{aligned} \frac{\partial T^{*d}}{\partial {\xi }_a}=\frac{\partial }{\partial {\xi }_a}\left( { \varvec{T}^{*d T}{\varvec{M}}}_d\prod ^{ND}_{\small { \begin{array}{c} j=1 \\ j\ne d \end{array}} }{{{\varvec{M}}}^T_j{{\varvec{T}}}_{jn}}\right) =\varvec{T}^{*d T}\overbrace{{{\varvec{M}}}_d}^{a}\prod ^{ND}_{\small { \begin{array}{c} j=1 \\ j\ne d \end{array}} }{\overbrace{{{\varvec{M}}}^T_j}^{a}{{\varvec{T}}}_{jn}}, \end{aligned}$$

(65)

respectively.