Parametric solutions involving geometry integrated with computer-aided design

TL;DR

This paper introduces a general framework for computing parametric solutions in geometry using NURBS within CAD environments, enabling efficient sensitivity analysis and integration with engineering workflows.

Contribution

It presents a novel, general approach for off-line parametric solutions involving NURBS geometries, applicable to 2D and 3D, with integration into CAD and efficient PGD-based computations.

Findings

Valid for 2D and 3D geometries

Achieves optimal convergence rates for flow problems

Circumvents curse of dimensionality with PGD

Abstract

The main objective of this work is to describe a general and original approach for computing an off-line solution for a set of parameters describing the geometry of the domain. That is, a solution able to include information for different geometrical parameter values and also allowing to compute readily the sensitivities. Instead of problem dependent approaches, a general framework is presented for standard engineering environments where the geometry is defined by means of NURBS. The parameters controlling the geometry are now the control points characterising the NURBS curves or surfaces. The approach proposed here, valid for 2D and 3D scenarios, allows a seamless integration with CAD preprocessors. The proper generalised decomposition (PGD), which is applied here to compute explicit geometrically parametrised solutions, circumvents the curse of dimensionality. Moreover, optimal convergence rates are shown for PGD approximations of incompressible flows.

Figures

Equations

$$\{\begin{array}{rlrl}-\text{boldmath∇unboldmath}\cdot (\mathbf{K}\text{boldmath∇unboldmath}{u}^{\text{boldmathμunboldmath}})& =s& & \text{in Ωboldmathμunboldmath,}\\ u^{\text{boldmathμunboldmath}}& =u_D& & \text{on ΓDboldmathμunboldmath,}\\ \text{boldmathnunboldmath}\cdot (\mathbf{K}\text{boldmath∇unboldmath}{u}^{\text{boldmathμunboldmath}})& =t& & \text{on ΓNboldmathμunboldmath,}\end{array}$$

$$a(\text{boldmathμunboldmath};u^{\text{boldmathμunboldmath}},v)=\ell (\text{boldmathμunboldmath};v),$$

$$a(\text{\boldmath$\mu$\unboldmath};u,v):=\bigl{(}\text{\boldmath$\nabla$\unboldmath}v,\mathbf{K}\text{\boldmath$\nabla$\unboldmath}u\bigr{)}_{\!\Omega^{\text{\boldmath$\mu$\unboldmath}}}\text{ and }\;\ell(\text{\boldmath$\mu$\unboldmath};v):=\bigl{(}v,s\bigr{)}_{\!\Omega^{\text{\boldmath$\mu$\unboldmath}}}+\bigl{\langle}v,t\bigr{\rangle}_{\Gamma^{\text{\boldmath$\mu$\unboldmath}}_{N}},$$

$$\bigl{(}u,v\bigr{)}_{\!\Omega}=\int_{\Omega}u\,v\,d\Omega\text{, }\;\quad\bigl{(}\text{\boldmath$u$\unboldmath},\text{\boldmath$v$\unboldmath}\bigr{)}_{\!\Omega}=\int_{\Omega}\text{\boldmath$u$\unboldmath}\cdot\text{\boldmath$v$\unboldmath}\,d\Omega\text{, and }\;\quad\bigl{\langle}u,v\bigr{\rangle}_{\Gamma}=\int_{\Gamma}{u\,v}\,d\Gamma,$$

$$A(u,v)=L(v),\forall v\in {\mathcal{H}}_{{\Gamma }_D^{\text{boldmathμunboldmath}}}^1\otimes {\mathcal{L}}_2(I_1)\otimes {\mathcal{L}}_2(I_2)\otimes \cdots \otimes {\mathcal{L}}_2(I_{{\text{n}}_{\text{pa}}}),$$

$$\begin{array}{rl}A(u,v)& :={\int }_{{\mathcal{I}}_1}{\int }_{{\mathcal{I}}_2}\cdots {\int }_{{\mathcal{I}}_{{\text{n}}_{\text{pa}}}}a(\text{boldmathμunboldmath};u,v)d{\mu }_{{\text{n}}_{\text{pa}}}\cdots d{\mu }_{2}d{\mu }_1\text{ and }\\ L(v)& :={\int }_{{\mathcal{I}}_1}{\int }_{{\mathcal{I}}_2}\cdots {\int }_{{\mathcal{I}}_{{\text{n}}_{\text{pa}}}}\ell (\text{boldmathμunboldmath};v)d{\mu }_{{\text{n}}_{\text{pa}}}\cdots d{\mu }_{2}d{\mu }_1.\end{array}$$

$$a(\text{\boldmath$\mu$\unboldmath};u,v)=\sum_{q=1}^{Q}\biggl{(}\prod_{i=1}^{\texttt{n}_{\texttt{pa}}}\Theta^{q}_{i}(\mu_{i})\biggr{)}a^{q}(u,v).$$

$$u\approx u_{\text{PGD}}^n=\sum _{m=1}^n{\psi }_m^u(\text{boldmathμunboldmath})w_m(\text{boldmathxunboldmath}).$$

$${\text{boldmathMunboldmath}}_{\text{boldmathμunboldmath}}:\Omega \times \text{boldmathIunboldmath}$$

$$(\text{boldmathXunboldmath},\text{boldmathμunboldmath})$$

$$\begin{split}a(\text{\boldmath$\mu$\unboldmath};u,v)&=\bigl{(}\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}}^{-1}\text{\boldmath$\nabla$\unboldmath}\!\!_{\text{\boldmath$X$\unboldmath}}v,\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\mathbf{K}\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}}^{-1}\text{\boldmath$\nabla$\unboldmath}\!\!_{\text{\boldmath$X$\unboldmath}}u\bigr{)}_{\!\Omega}=\bigl{(}\text{\boldmath$\nabla$\unboldmath}\!\!_{\text{\boldmath$X$\unboldmath}}v,\mathbf{H}_{\text{\boldmath$\mu$\unboldmath}}\,\text{\boldmath$\nabla$\unboldmath}\!\!_{\text{\boldmath$X$\unboldmath}}u\bigr{)}_{\!\Omega}\text{ and }\\ \ell(\text{\boldmath$\mu$\unboldmath};v)&=\bigl{(}v,\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\,s\bigr{)}_{\!\Omega}+\bigl{\langle}v,\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\,t\bigr{\rangle}_{\Gamma_{N}},\end{split}$$

$${\mathbf{H}}_{\text{boldmathμunboldmath}}:=\frac{\mathrm{adj}({\mathbf{J}}_{\text{boldmathμunboldmath}}^T)\mathbf{K}\mathrm{adj}({\mathbf{J}}_{\text{boldmathμunboldmath}})}{\det ({\mathbf{J}}_{\text{boldmathμunboldmath}})},$$

$$a(\text{\boldmath$\mu$\unboldmath};\text{\boldmath$u$\unboldmath},\text{\boldmath$v$\unboldmath})=\bigl{(}\text{\boldmath$\nabla$\unboldmath}\text{\boldmath$v$\unboldmath},\underline{\underline{\mathbf{C}}}:\text{\boldmath$\nabla$\unboldmath}\text{\boldmath$u$\unboldmath}\bigr{)}_{\!\Omega^{\text{\boldmath$\mu$\unboldmath}}}=\bigl{(}\text{\boldmath$\nabla$\unboldmath}\!\!_{\text{\boldmath$X$\unboldmath}}\text{\boldmath$v$\unboldmath},\underline{\underline{\mathbf{\hat{C}}}}^{\text{\boldmath$\mu$\unboldmath}}:\text{\boldmath$\nabla$\unboldmath}\!\!_{\text{\boldmath$X$\unboldmath}}\text{\boldmath$u$\unboldmath}\bigr{)}_{\!\Omega},$$

$${\text{boldmathMunboldmath}}_{\text{boldmathμunboldmath}}(\text{boldmathXunboldmath},\text{boldmathμunboldmath})=\text{boldmathXunboldmath}+\sum _{m=1}^N{\psi }_m(\text{boldmathμunboldmath}){\text{boldmathdunboldmath}}^m(\text{boldmathXunboldmath}),$$

$${\mathbf{J}}_{\text{boldmathμunboldmath}}(\text{boldmathXunboldmath},\text{boldmathμunboldmath})=\frac{\partial \text{boldmathxunboldmath}}{\partial \text{boldmathXunboldmath}}(\text{boldmathXunboldmath},\text{boldmathμunboldmath})={\mathbf{I}}_{{\text{n}}_{\text{sd}}}+\sum _{m=1}^N{\psi }_m{\mathbf{A}}_m,$$

$$\mathrm{adj}(\mathbf{A})=\sum _{s=0}^{n-1}{\mathbf{A}}^s\sum _{k_1,k_2,\dots ,k_{n-1}}\prod _{l=1}^{n-1}\frac{(-1{)}^{k_l+1}}{l^{k_l}{k}_l!}\mathrm{tr}({\mathbf{A}}^l{)}^{k_l}$$

$$\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})=\Bigl{(}1+\sum_{m=1}^{N}\psi_{m}(\text{\boldmath$\mu$\unboldmath})\bigl{[}\mathbf{A}_{m}\bigr{]}_{11}\Bigr{)}\Bigl{(}1+\sum_{m=1}^{N}\psi_{m}(\text{\boldmath$\mu$\unboldmath})\bigl{[}\mathbf{A}_{m}\bigr{]}_{22}\Bigr{)}\\ -\Bigl{(}\sum_{m=1}^{N}\psi_{m}(\text{\boldmath$\mu$\unboldmath})\bigl{[}\mathbf{A}_{m}\bigr{]}_{21}\Bigr{)}\Bigl{(}\sum_{m=1}^{N}\psi_{m}(\text{\boldmath$\mu$\unboldmath})\bigl{[}\mathbf{A}_{m}\bigr{]}_{12}\Bigr{)},$$

$$\mathrm{adj}(\mathbf{A})=\mathrm{tr}(\mathbf{A}){\mathbf{I}}_2-\mathbf{A},$$

$$\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})=\mathbf{I}_{2}+\sum_{m=1}^{N}\psi_{m}\,\operatorname{adj}\bigl{(}[\partial\text{\boldmath$d$\unboldmath}^{m}/\partial\text{\boldmath$X$\unboldmath}]\bigr{)}=\mathbf{I}_{2}+\sum_{m=1}^{N}\psi_{m}\,\operatorname{adj}(\mathbf{A}_{m}).$$

$$\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\;\mathbf{H}_{\text{\boldmath$\mu$\unboldmath}}=\mathbf{K}+\sum_{m=1}^{N}\psi_{m}\,\mathbf{K}\,\operatorname{adj}(\mathbf{A}_{m})+\sum_{m=1}^{N}\psi_{m}\bigl{[}\mathbf{K}\,\operatorname{adj}(\mathbf{A}_{m})\bigr{]}^{T}\\ +\sum_{m=1}^{N}\sum_{l=1}^{N}\psi_{m}\psi_{l}\,\operatorname{adj}(\mathbf{A}_{m})^{T}\,\mathbf{K}\,\operatorname{adj}(\mathbf{A}_{l}).$$

$$\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})=\sum_{\sigma\in\mathcal{S}_{3}}\operatorname{sgn}(\sigma)\prod_{i=1}^{3}\Bigl{(}1{+}\sum_{m=1}^{N}\psi_{m}(\text{\boldmath$\mu$\unboldmath})\bigl{[}\mathbf{A}_{m}\bigr{]}_{i\sigma(i)}\Bigr{)},$$

$$\operatorname{adj}(\mathbf{A})=\tfrac{1}{2}\bigl{[}\bigl{(}\operatorname{tr}(\mathbf{A})\bigr{)}^{2}-\operatorname{tr}(\mathbf{A}^{2})\bigr{]}\mathbf{I}_{3}-\operatorname{tr}(\mathbf{A})\mathbf{A}+\mathbf{A}^{2},$$

$$\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})=\mathbf{I}_{3}+\sum_{m=1}^{N}\psi_{m}\bigl{[}\operatorname{tr}(\mathbf{A}_{m})\mathbf{I}_{3}-\mathbf{A}_{m}\bigr{]}\\ +\sum_{m=1}^{N}\sum_{l=1}^{N}\psi_{m}\psi_{l}\bigl{[}\frac{1}{2}\bigl{(}\operatorname{tr}(\mathbf{A}_{m})\operatorname{tr}(\mathbf{A}_{l})-\operatorname{tr}(\mathbf{A}_{m}\mathbf{A}_{l})\bigr{)}\mathbf{I}_{3}-\operatorname{tr}(\mathbf{A}_{m})\mathbf{A}_{l}+\mathbf{A}_{m}\mathbf{A}_{l}\bigr{]}.$$

$${\mathbf{H}}_{\text{boldmathμunboldmath}}\approx {\mathbf{H}}_{\text{boldmathμunboldmath}}^{\text{sep}}(\text{boldmathxunboldmath},\text{boldmathμunboldmath})=\sum _{m=1}^{N_H}{\psi }_m^{\mathbf{H}}(\text{boldmathμunboldmath}){\mathbf{A}}_m^{\mathbf{H}}(\text{boldmathXunboldmath}),$$

$$\bigl{(}\text{\boldmath$v$\unboldmath},\mathbf{H}^{\texttt{sep}}_{\text{\boldmath$\mu$\unboldmath}}\bigr{)}_{\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath}}=\bigl{(}\text{\boldmath$v$\unboldmath},\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}}^{T})\,\mathbf{K}\,\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})/\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\bigr{)}_{\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath}}$$

$$\bigl{(}\text{\boldmath$v$\unboldmath},\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\,\psi^{\mathbf{H}}_{n}\,\mathbf{A}^{\mathbf{H}}_{n}\bigr{)}_{\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath}}=\bigl{(}\text{\boldmath$v$\unboldmath},\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}}^{T})\,\mathbf{K}\,\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\bigr{)}_{\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath}}\\ -\bigl{(}\text{\boldmath$v$\unboldmath},\det(\mathbf{J}_{\!\text{\boldmath$\mu$\unboldmath}})\,\sum_{m=1}^{n-1}\psi^{\mathbf{H}}_{m}\,\mathbf{A}^{\mathbf{H}}_{m}\bigr{)}_{\Omega\times\text{\boldmath$\mathcal{I}$\unboldmath}},$$

$$\partial {\Omega }^{\text{boldmathμunboldmath}}=\bigcup _{j=1}^M{\text{boldmathCunboldmath}}_j^{\text{boldmathμunboldmath}}([0,1]).$$

$${\text{boldmathCunboldmath}}_{}(\lambda )=\sum _{i=0}^{{\text{n}}_{\text{cp}}}{\text{boldmathBunboldmath}}_i{R}_i(\lambda )\lambda \in [0,1]$$

$$R_{i}(\lambda)=\nu_{i}\,C_{i}^{q}(\lambda)\biggm{/}\left(\sum_{i=0}^{\texttt{n}_{\texttt{cp}}}\nu_{i}\,C_{i}^{q}(\lambda)\right).$$

$$C_i^0(\lambda )$$

$$C_i^k(\lambda )$$

$$\Lambda =\{\underset{q+1}{\underbrace{0,\dots ,0}},{\lambda }_{q+1},\dots ,{\lambda }_{{\text{n}}_{\text{k}}-q-1},\underset{q+1}{\underbrace{1,\dots ,1}}\},$$

Full text

Keywords: Geometry parametrisation; Reduced order model; Computer-aided design (CAD); Proper generalised decomposition (PGD).

1 Introduction

The current role of computational simulations in modern engineering design is limited by the complexity of the simulations that are required, particularly during the final stages of a design. The main problem is motivated by the number of configurations that need to be tested (e.g. loads, boundary conditions, material parameters and geometric configurations).

One alternative to decrease the computational complexity in this scenario is to introduce a reduced order model [27]. The main idea involves projecting the governing equations describing the full model onto a space with lower dimension that is described using a reduced order basis. Well known methods to produce reduced order basis are Krylov-based methods [16], the reduced basis method [29] and the proper orthogonal decomposition (POD) [4, 20]. More recently, the proper generalised decomposition (PGD) [8, 9] has gained popularity due to its ability to build reduced basis with no prior knowledge of the solution. The PGD starts by considering the solution not only as a function of the standard coordinates (i.e., space and time) but also of any parameter of interest (e.g. boundary conditions, external loads, material parameters). The problem involving a range of all the desired parameters can be solved at the cost of several problems of the same size as the original problem for a particular choice of the parameters. This expensive calculation, usually referred as the off-line stage, is performed only once, usually making use of high performance computing resources, to build the reduced order basis that is written as an explicit function of the coordinates (space and time) and the parameters (i.e. a computational vademecum, see [10]). Then, the on-line stage consists of just a particularisation of the solution by using, in the simplest case, an interpolation of the already computed results.

The PGD has been successfully applied to numerous multi-dimensional problems involving boundary conditions, material parameters and external loads as extra coordinates, to name a few, see [8, 9] and references therein. The application to problems involving geometrically parametrised domains is generally more challenging and the existing work is usually limited to simple geometries where the extra coordinates are scalings: the length of an interval in one-dimensional problems [11], the thickness of extruded geometries [19, 5] or problem specific parameters [17]. More recently, an approach based on an initial subdivision of the computational domain in macro-elements was proposed in [2]. This idea was extended to domains with interfaces in [34] and has been also applied to an engineering design process in [12].

In this paper, a new approach to incorporate the geometric parameters as extra coordinates in a PGD framework is proposed. The objective is to produce a general methodology that enables to obtain the solution of a particular problem when the geometry of the domain is parametrised using the NURBS boundary representation of the domain, as usually done in a computer-aided design (CAD) environment. The control points of the NURBS entities defining the boundary are considered as extra coordinates and a mapping between a reference domain and the current configuration is proposed by using a solid mechanics analogy. An explicit and separated representation of the mapping is derived and the application of the proposed methodology to Stokes flow problems is presented using examples of increasing difficulty in two and three dimensions. Contrary to the approach in [2, 34], the methodology proposed in this paper considers geometric parameters that are independent on the spatial discretisation, i.e. the control points of the NURBS entities defining the boundary representation of the domain. In addition, the technique presented here enables the solution of the multi-dimensional problem to be computed using high-order finite elements whereas the technique introduced in [2, 34] requires an affine mapping between a reference element and the macro-elements that are used to parametrise the geometry.

The structure of the remainder of the paper is as follows. Section 2 presents the problem statement using the Poisson equation on a geometrically parametrised domain and summarises the application of the PGD. Section 3 describes in detail the proposed technique to build a generalised solution assuming that a mapping between a reference configuration and the current one can be written in separated form. In Section 4, a methodology to integrate this approach in a CAD environment is presented. Considering the geometric parameters as the control points of the NURBS entities describing the boundary of the domain, the methodology to build a mapping that can be explicitly written in separated form is detailed. Finally, Section 5 presents a series of numerical examples of increasing difficulty involving the solution of Stokes flow problems in two and three dimensions.

2 Problem statement and geometrically parametrised solutions

2.1 The Poisson equation on a parametrised domain

The methodology proposed here can be directly extended to second order linear problems such as the Stokes problem that is studied in the examples. However, in order to simplify the presentation, the heat problem (Poisson) is presented in detail for a parametrised domain \Omega^{\text{\boldmath\mu\unboldmath}}\subset\mathbb{R}^{\texttt{n}_{\texttt{sd}}} (with $\texttt{n}_{\texttt{sd}}$ number of spatial dimensions), whose boundary \partial\Omega^{\text{\boldmath\mu\unboldmath}} is characterised by a set of geometric parameters \text{\boldmath\mu\unboldmath}\in\text{\boldmath\mathcal{I}\unboldmath}\subset\mathbb{R}^{\texttt{n}_{\texttt{pa}}} (with $\texttt{n}_{\texttt{pa}}$ number of parameters characterising the geometry) and is partitioned into Dirichlet, \Gamma^{\text{\boldmath\mu\unboldmath}}_{D}, and Neumann, \Gamma^{\text{\boldmath\mu\unboldmath}}_{N}, frontiers such that \overline{\partial\Omega}^{\text{\boldmath\mu\unboldmath}}=\overline{\Gamma}^{\text{\boldmath\mu\unboldmath}}_{D}\cup\overline{\Gamma}^{\text{\boldmath\mu\unboldmath}}_{N} and \Gamma^{\text{\boldmath\mu\unboldmath}}_{D}\cap\Gamma^{\text{\boldmath\mu\unboldmath}}_{N}=\emptyset. Note that the set \text{\boldmath\mathcal{I}\unboldmath}\subset\mathbb{R}^{\texttt{n}_{\texttt{pa}}}, which characterises the admissible range for parameters $\mu$, can be defined as the Cartesian combination of the range for each parameter, namely, \text{\boldmath\mathcal{I}\unboldmath}:=\mathcal{I}_{1}\times\mathcal{I}_{2}\times\dotsb\times\mathcal{I}_{\texttt{n}_{\texttt{pa}}} with $\mu_{i}\in\mathcal{I}_{i}$ for $i=1,\dotsc,\texttt{n}_{\texttt{pa}}$.

For each set of parameter $\mu$, the objective is to determine the parametric solution u^{\text{\boldmath\mu\unboldmath}}(\text{\boldmathx\unboldmath}), with \text{\boldmathx\unboldmath}\in\Omega^{\text{\boldmath\mu\unboldmath}}, of the boundary value problem

[TABLE]

where $\mathbf{K}$ is the thermal conductivity (symmetric and positive definite) matrix, $s$ is a source term, $u_{D}$ is the imposed temperature, $t$ is the imposed heat normal flux and $n$ is the outward unit normal vector. The standard variational form of the previous problem reads: find u^{\text{\boldmath\mu\unboldmath}}\in\mathcal{V}^{\text{\boldmath\mu\unboldmath}} for all v\in\mathcal{H}^{1}_{\Gamma^{\text{\boldmath\mu\unboldmath}}_{D}} such that

[TABLE]

where

[TABLE]

denote, respectively, the $\mathcal{L}_{2}$ product of scalar/vector functions in $\Omega$ and its traces over $\Gamma$.

2.2 The multi-dimensional parametric problem

There are different alternatives to obtain, for any given set of parameters $\mu$, an approximation to the solution u^{\text{\boldmath\mu\unboldmath}}(\text{\boldmathx\unboldmath}) of problem (1). The obvious option of solving a new problem for every instance of $\mu$ is feasible but too costly. A standard strategy to reduce the cost is to pre-compute off-line some representative samples of the parametric family of solutions (e.g. snapshots for reduced basis methods, principal components for POD). Then, any other instance is computed on-line with a small computational overhead. Here, the PGD is preferred because the off-line phase provides an explicit description of the parametric solution, i.e. a computational vademecum see [10]. Thus, in spite of an off-line phase more involved, the on-line phase is a simple functional evaluation with a negligible computational overhead.

In practice, this can be interpreted as taking $\mu$ as additional independent variables (or parametric coordinates) instead of problem parameters. Hence, the unknown temperature field is not interpreted any more as a parametric solution, denoted as u^{\text{\boldmath\mu\unboldmath}}, but it is seen now as a function in a larger dimensional space and it is written as u(\text{\boldmathx\unboldmath},\text{\boldmath\mu\unboldmath}) with (\text{\boldmathx\unboldmath},\text{\boldmath\mu\unboldmath})\in\Omega^{\text{\boldmath\mu\unboldmath}}\times\text{\boldmath\mathcal{I}\unboldmath}. Consequently, formally $u$ lies in a tensor product space, namely, u\in\mathcal{V}^{\text{\boldmath\mu\unboldmath}}\otimes\mathcal{L}_{2}(\mathcal{I}_{1})\otimes\mathcal{L}_{2}(\mathcal{I}_{2})\otimes\dotsb\otimes\mathcal{L}_{2}(\mathcal{I}_{\texttt{n}_{\texttt{pa}}}). A standard weighted residuals approach, with integrals in \Omega^{\text{\boldmath\mu\unboldmath}}\times\text{\boldmath\mathcal{I}\unboldmath} and the usual integration by parts only in \Omega^{\text{\boldmath\mu\unboldmath}} produces a weak form in this multi-dimensional setup. Namely, find u\in\mathcal{V}^{\text{\boldmath\mu\unboldmath}}\otimes\mathcal{L}_{2}(\mathcal{I}_{1})\otimes\mathcal{L}_{2}(\mathcal{I}_{2})\otimes\dotsb\otimes\mathcal{L}_{2}(\mathcal{I}_{\texttt{n}_{\texttt{pa}}}) such that

[TABLE]

Obviously, the number of dimensions of the solution domain increases with the number of parameters. To circumvent the curse of dimensionality, the PGD approach [3, 11, 8, 9] is employed here. This approach assumes a separable structure in the function that approximates $u$. Note that the tensor product space \mathcal{V}^{\text{\boldmath\mu\unboldmath}}\otimes\mathcal{L}_{2}(\mathcal{I}_{1})\otimes\mathcal{L}_{2}(\mathcal{I}_{2})\otimes\dotsb\otimes\mathcal{L}_{2}(\mathcal{I}_{\texttt{n}_{\texttt{pa}}}) inherits the multi-dimensional complexity of the problem and, in principle, does not assume separability of the functions.

Moreover, the solution of (3a) requires an affine parameter dependence of the different forms. This is standard in reduced order methods and it is very well discussed in [23, 28]. More precisely, it is required that the different forms are expressed (or at least well approximated) by the sum of products of parameter-dependent functions and parameter-independent operators, for instance

[TABLE]

Note that the forms $a^{q}(u,v)$ do not depend on the parameters (in particular, they are integrated over domains parameter independent). In fact, finding the affine parameter dependence of (2b) is a major concern in subsequent sections and will enable to obtain a separated approximation of the solution, namely

[TABLE]

3 Separated spatial mapping to determine generalised solutions

If the affine parameter dependence must be enforced, it is necessary to integrate in space over domains not depending on the parameters. Note that forms a(\text{\boldmath\mu\unboldmath};u,v) and \ell(\text{\boldmath\mu\unboldmath};v), see (2b), are integrated on spacial parametrised domains (domains depending on parameters $\mu$). As suggested in [2, 34], a mapping \text{\boldmath\mathcal{M}\unboldmath}_{\text{\boldmath\mu\unboldmath}} is necessary (not sufficient) in order to have an affine parameter dependence. This mapping transforms a reference domain $\Omega$ into the geometrically parametrised (“deformed”) domain \Omega^{\text{\boldmath\mu\unboldmath}}, namely

[TABLE]

As classically in computational mechanics, the reference configuration $\Omega$ is associated to a reference coordinate system denoted by $X$, whereas the distorted domain \Omega^{\text{\boldmath\mu\unboldmath}} will be associated to the spatial description $x$. Following this analogy, a displacement field $d$ is used to relate both coordinate systems.

The mapping \text{\boldmath\mathcal{M}\unboldmath}_{\text{\boldmath\mu\unboldmath}} can be defined in an ad hoc manner for each problem [2, 12] or by a more general strategy [2, 34], but always has to induce an affine parameter dependence of forms (2b). The introduction of this mapping allows to rewrite (2b) as integrals over the reference computational domain $\Omega$ and its corresponding boundary $\partial\Omega$, partitioned into the Dirichlet, $\Gamma_{D}$, and Neumann, $\Gamma_{N}$, boundaries, all independent of parameter $\mu$. Namely,

[TABLE]

where [\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}]_{ij}=[\partial x_{j}/\partial X_{i}] is the Jacobian matrix of the mapping \text{\boldmath\mathcal{M}\unboldmath}_{\text{\boldmath\mu\unboldmath}} and, more importantly, it is the only element in the above equations that depends upon parameters $\mu$. Moreover, in order to compact the notation a new matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} is introduced,

[TABLE]

where the definition of the adjoint operator, \operatorname{adj}(\text{\boldmathA\unboldmath})=\det(\text{\boldmathA\unboldmath})\,\text{\boldmathA\unboldmath}^{-1}, has been used.

Note that even for ad hoc [2] or a piecewise linear [2, 34] mappings, the affine parameter dependence of the different forms in (5) must be now determined.

Remark 1** (Mapping for the Stokes problem).**

Despite being more cumbersome, the same rationale can be applied for the Stokes problem without any extra conceptual challenge. The bilinear viscosity form is reproduced here for illustration purposes

[TABLE]

where $\text{C}_{ijkl}=\nu\,\delta_{ik}\,\delta_{jl}$, \hat{\text{C}}_{ijkl}^{\text{\boldmath\mu\unboldmath}}=\nu\det(\mathbf{J_{\text{\boldmath\mu\unboldmath}}})\,\delta_{ik}\,[\mathbf{J_{\text{\boldmath\mu\unboldmath}}^{-T}}]_{js}\,[\mathbf{J_{\text{\boldmath\mu\unboldmath}}^{-1}}]_{sl} and $\nu>0$ is the kinematic viscosity. Note that $[\mathbf{A}]_{ij}$ denotes the component $ij$ of a matrix $\mathbf{A}$.

3.1 Separated displacements

The methodology proposed here allows to obtain an affine parameter dependence (namely, a separable expression for \mathbf{H}_{\text{\boldmath\mu\unboldmath}}) quasi-analytically if the mapping can be written (or at least, well approximated) with a separated representation, that is as a sum of separated terms, namely

[TABLE]

which induces a separated Jacobian matrix

[TABLE]

where \mathbf{A}_{m}:=[\partial\text{\boldmathd\unboldmath}^{m}/\partial\text{\boldmathX\unboldmath}] for $m=1,\dotsc,N$.

3.2 Affine parameter dependence

The affine parameter dependence of a(\text{\boldmath\mu\unboldmath};u,v) as defined in (5), is in practice determined by obtaining a separated expression of \mathbf{H}_{\text{\boldmath\mu\unboldmath}}. In order to obtain a separated expression for matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, both \operatorname{adj}(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}) and \det(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}) are analytically separated. Then, the Higher-Order PGD-Projection proposed in [22] is used to obtain a compact separation for \mathbf{H}_{\text{\boldmath\mu\unboldmath}}.

The separated representation for the determinant, \det(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}), can be obtained using Leibniz formula from equation (8), whereas for \operatorname{adj}(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}), the Leverrier’s algorithm [15] is employed. This method is a consequence of the Cayley-Hamilton theorem [6, 18] and the Newton’s identities [21]. It enables to express the adjoint of a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ in terms of its trace and its powers, namely

[TABLE]

where $k_{l}\in\mathbb{N}_{0}$ and $s+\sum_{l=1}^{n-1}lk_{l}=n-1$. In practice, the two cases of interest are 2D ($n=2$) and 3D ($n=3$) problems, which are detailed next.

3.3 Two-dimensional approach

In 2D the Jacobian \mathbf{J}_{\!\text{\boldmath\mu\unboldmath}} is a $2\times 2$ matrix. From (8), the separated expression for its determinant is

[TABLE]

again $[\mathbf{A}]_{ij}$ denotes the component $ij$ of the matrix $\mathbf{A}$.

To obtain the expression for the adjoint requires to particularise (9), namely

[TABLE]

which is a linear mapping because the trace is also linear. Thus, the Jacobian in separated form as presented in (8) can be written as

[TABLE]

Consequently, recalling that $\operatorname{adj}(\mathbf{A}^{T})=\operatorname{adj}(\mathbf{A})^{T}$, the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} in (6) can be rewritten as

[TABLE]

The separated expression of \mathbf{H}_{\text{\boldmath\mu\unboldmath}} is efficiently obtained with a numerical Higher-Order PGD-Projection [22].

3.4 Three-dimensional approach

Following the previous rationale, in 3D the Jacobian \mathbf{J}_{\!\text{\boldmath\mu\unboldmath}} is a $3\times 3$ matrix. The separated expression for the determinant, \det(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}), is obtained using Leibniz formula as

[TABLE]

where $\mathcal{S}_{3}$ is the set of the six permutations of the integers $\{1,2,3\}$, where the element in position $i$ after the reordering $\sigma$ is denoted $\sigma(i)$, and $\operatorname{sgn}(\sigma)$ denotes the signature of $\sigma$ (i.e. $+1$ for even $\sigma$ and $-1$ for odd $\sigma$). Note that such a separation will induce sums of order $N^{3}$.

The adjoint of the Jacobian is also separated by means of particularising (9) to 3D,

[TABLE]

and given the separation of the Jacobian in (8) implies

[TABLE]

As done in 2D, once the adjoint is separated, from (6), a separated expression for \det(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}})\mathbf{H}_{\text{\boldmath\mu\unboldmath}}=\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}^{T})\,\mathbf{K}\,\operatorname{adj}(\mathbf{J}_{\!\text{\boldmath\mu\unboldmath}}) can be computed. In this case, it consists of $N^{4}$ terms. As previously noted for 2D, the separated expression of \mathbf{H}_{\text{\boldmath\mu\unboldmath}} is obtained with the Higher-Order PGD-Projection described next.

3.5 The Higher-Order PGD-Projection

Finally, to obtain a separated approximation of \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, namely

[TABLE]

equation (6) has to be solved. This is done by means of the Higher-Order PGD-Projection [22], which in practice obtains the separated approximation with an $\mathcal{L}_{2}$ projection, say

[TABLE]

for all $v$ in a suitable space. In practice, the PGD strategy is implemented. Thus, a greedy approach is used. For computational efficiency, it is critical to use the exact separated expressions for the determinant, see (10) and (12), and the adjoint of the Jacobian, see (11) and (13). In practice, each mode of \mathbf{H}_{\text{\boldmath\mu\unboldmath}}^{\texttt{sep}} is obtained by solving

[TABLE]

where now, as usual in PGD, $v$ lives in the tangent space of the mode. Note that the Higher-Order PGD-Projection allows many parameters without any prior knowledge, is computationally efficient, and the precision of this approximation (i.e. the total number of terms $N_{H}$) can be controlled by the user.

4 Integration within a CAD environment

In this section a procedure to integrate the methodology described in the previous section within a CAD environment is proposed. To simplify the presentation, the two dimensional case is presented here and the details for the extension to three dimensional domains are given in A.

The boundary of the parametrised domain, \partial\Omega^{\text{\boldmath\mu\unboldmath}}, is assumed to be described by a set of NURBS curves \{\text{\boldmathC\unboldmath}_{\!j}^{\text{\boldmath\mu\unboldmath}}\}_{j=1,\dotsc,M}, being $M$ the total number of curves, namely

[TABLE]

Next, the necessary concepts about NURBS curves are briefly recalled and the proposed strategy to build a geometric mapping \text{\boldmath\mathcal{M}\unboldmath}_{\text{\boldmath\mu\unboldmath}} that can be written in the separated form (7) is presented in detail.

4.1 NURBS curves

A $q$th-degree non-uniform rational B-spline (NURBS) curve is a piecewise rational function defined in parametric form as

[TABLE]

where \{\text{\boldmathB\unboldmath}_{i}\} are the coordinates of the $\texttt{n}_{\texttt{cp}}+1$ control points (forming the control polygon) and $\{R_{i}(\lambda)\}$ are rational basis functions defined as

[TABLE]

In the above expression $\{\nu_{i}\}$ are the control weights associated to the control points and $\{C_{i}^{q}(\lambda)\}$ are the normalized B-spline basis functions of degree $q$, which are defined recursively by

[TABLE]

for $k=1,\dotsc,q$ and where $\lambda_{i}$ (for $i=0,\dotsc,\texttt{n}_{\texttt{k}}$) are the knots or breakpoints, which are assumed ordered $0\leq\lambda_{i}\leq\lambda_{i+1}\leq 1$. They form the so-called knot vector,

[TABLE]

which uniquely describes the B-spline basis functions. The number of control points, $\texttt{n}_{\texttt{cp}}+1$, and knots, $\texttt{n}_{\texttt{k}}+1$, are related to the degree of the parametrisation, $q$, by the relation $\texttt{n}_{\texttt{k}}=\texttt{n}_{\texttt{cp}}+q+1$, see [25] for more details.

Figure 1 shows an example of a two dimensional domain \Omega^{\text{\boldmath\mu\unboldmath}} where the boundary is described by five NURBS curves.

The curve \text{\boldmathC\unboldmath}_{\!3} in Figure 1 (b) is represented in Figure 2 with the corresponding control polygon, formed by six control points and the breakpoints. The knot vector for this curve is given by

[TABLE]

4.2 Geometric parameters

The geometric parameters \text{\boldmath\mu\unboldmath}\in\text{\boldmath\mathcal{I}\unboldmath}\subset\mathbb{R}^{\texttt{n}_{\texttt{pa}}} are defined as the variations of the original coordinates of the control points of the NURBS curves describing the boundary. More precisely, for each NURBS curve \text{\boldmathC\unboldmath}_{\!j}, with $j=1,\dotsc,M$, having $\texttt{n}_{\texttt{cp}}^{j}+1$ control points, the undisturbed boundary is characterised by the coordinates of the control points: \text{\boldmathB\unboldmath}_{i}^{j}, for $i=0,\dotsc,\texttt{n}_{\texttt{cp}}^{j}$. This configuration will be used as the reference one in $\Omega$ and will be associated to a reference coordinate system denoted by $X$. The distorted domain, \Omega^{\text{\boldmath\mu\unboldmath}}, will be associated to the spatial description $x$. The boundary in the spatial domain, \partial\Omega^{\text{\boldmath\mu\unboldmath}}, is defined by the position of the displaced control points, namely \text{\boldmathB\unboldmath}_{i}^{j}+\delta\!\text{\boldmathB\unboldmath}_{i}^{j}. The displacement range for each control point is characterised by

[TABLE]

In fact, each displacement of a control point $i$ on the $j$-th NURBS, \delta\!\text{\boldmathB\unboldmath}_{i}^{j}, might depend upon the parameters and can be written as

[TABLE]

where \text{\boldmathe\unboldmath}_{i}, for $i=1,2$, are the unit coordinate vectors. Then \text{\boldmath\mu\unboldmath}\in\text{\boldmath\mathcal{I}\unboldmath}:=\text{\boldmath\mathcal{I}\unboldmath}^{1}\times\text{\boldmath\mathcal{I}\unboldmath}^{2}\times\dotsb\times\text{\boldmath\mathcal{I}\unboldmath}^{M}, where \text{\boldmath\mathcal{I}\unboldmath}^{j}=\text{\boldmath\mathcal{I}\unboldmath}^{j}_{1}\times\text{\boldmath\mathcal{I}\unboldmath}^{j}_{2}\times\cdots\text{\boldmath\mathcal{I}\unboldmath}^{j}_{\texttt{n}_{\texttt{cp}}^{j}+1} is the range of variation of the coordinates of the control points of the curve \text{\boldmathC\unboldmath}_{\!j}.

Consequently, the maximum number of geometric parameters is

[TABLE]

but in practical problems, the number of geometric parameters $\texttt{n}_{\texttt{pa}}$ is drastically lower than $m$ because not all the control points of all the curves are to be modified during the design stage.

Remark 2**.**

In a practical setting it is common to introduce some restrictions on the motion of the control points (viz. pure translations, rotations, expansions…), meaning that the motion of a set of control points can be expressed with a significantly low number of parameters.

4.3 Separated representation of the boundary displacement

The variation of a control point \text{\boldmathB\unboldmath}_{i}^{j} of a NURBS curve \text{\boldmathC\unboldmath}_{\!j}, namely \delta\!\text{\boldmathB\unboldmath}_{i}^{j}, changes the definition of the original curve only in the support of the basis function $R_{i}^{j}$, given by the subspace of the parametric space $[\lambda_{i},\lambda_{i+q^{j}+1}]$. The modified NURBS curve is parametrised by

[TABLE]

Figure 3 illustrates the effect of modifying the coordinates of one control point of a NURBS. The curve in red is the result of modifying the coordinates of the control point \text{\boldmathB\unboldmath}_{4}^{3} of the original curve in black, also depicted in Figure 2. It can be observed that both curves are identical in the interval [0,1/3] of the parametric space whereas they differ in the interval [1/3,1], which is the support of the basis function $R_{4}^{3}$ associated to the control point \text{\boldmathB\unboldmath}_{4}^{3}.

Given a computational mesh for the reference configuration, $\Omega$, the boundary mesh nodes affected by the motion of a control point can be easily identified. The procedure starts by finding the NURBS curve to which each boundary mesh node belongs and its associated parametric coordinate by using a standard NURBS point projection algorithm [25]. For each boundary mesh node \text{\boldmathX\unboldmath}_{k}, the index $j\in\{1,\dotsc,M\}$ and parametric coordinate $\lambda_{k}$ such that \text{\boldmathC\unboldmath}_{\!j}(\lambda_{k})=\text{\boldmathX\unboldmath}_{k} are computed. Then, for a deformed configuration, induced by the variation of a control point \text{\boldmathB\unboldmath}_{i}^{j} of a NURBS curve \text{\boldmathC\unboldmath}_{\!j}, namely \delta\!\text{\boldmathB\unboldmath}_{i}^{j}, the new position of each boundary node is computed as \text{\boldmathx\unboldmath}_{k}^{\text{\boldmath\mu\unboldmath}}=\text{\boldmathC\unboldmath}_{\!j}^{\text{\boldmath\mu\unboldmath}}(\lambda_{k}).

Therefore, the variation of the control points of a boundary curve \text{\boldmathC\unboldmath}_{\!j} induces a displacement of the boundary mesh nodes, namely

[TABLE]

where \mathcal{S}=\left\{l\in\left\{1,2,\ldots,\texttt{n}_{\texttt{mn}}\right\}:\text{\boldmathX\unboldmath}_{l}\in\partial\Omega\right\} is the index set of mesh nodes on the boundary of the computational domain and $\texttt{n}_{\texttt{mn}}$ is the total number of mesh nodes. Using the expression of the original and modified NURBS boundary curves, Equations (15) and (18) respectively, the displacement of the boundary mesh node \text{\boldmathX\unboldmath}_{k}=\text{\boldmathC\unboldmath}_{\!j}(\lambda_{k}) that belongs to the NURBS curve \text{\boldmathC\unboldmath}_{\!j} can be written in separated form as

[TABLE]

where the dependence of the displacements of the control points in terms of the parameters described in (16) has been used. Moreover, since the NURBS parameter $\lambda_{k}$ is only dependent on the spatial coordinates \text{\boldmathX\unboldmath}_{k}, and not on the geometric parameters $\mu$, the previous equation can be written as,

[TABLE]

which characterises the displacement of the boundary nodes and has a separated expression of the form

[TABLE]

Note that this expression is compatible with the desired structure of the displacements described in (7). In fact, it can be further compacted by means of the Higher-Order PGD-Projection [22] to obtain a more compact separation in the form of

[TABLE]

which now coincides with (7) but only for the boundary nodes. Precisely, the next section describes the extension to any point in the domain.

It is also important to note that the total number of terms in the separated representation (19) will be, in the majority of cases, much lower than $m$, defined in (17). In practical applications, a large percentage of mesh boundary nodes will not be affected by a variation of the control points of a particular NURBS curve describing the boundary, so the displacement will be zero for a large number of boundary nodes. More precisely, the set of mesh nodes affected by the variation of a control point \text{\boldmathB\unboldmath}_{i}^{j} of a NURBS curve \text{\boldmathC\unboldmath}_{\!j} can be defined as

[TABLE]

4.4 Separated representation of the geometric mapping

The proposed strategy to build a mapping \text{\boldmath\mathcal{M}\unboldmath}_{\text{\boldmath\mu\unboldmath}} between the reference configuration, $\Omega$, and the current configuration, \Omega^{\text{\boldmath\mu\unboldmath}}, consists on solving a solid mechanics problem. The reference configuration is assumed to be a linear elastic medium and the displacement of mesh boundary nodes, induced by the variation of NURBS control points, is interpreted as a Dirichlet boundary condition. The following problem governing the static deformation of $\Omega$ is considered

[TABLE]

where $f$ is an external force defined by the user. The stress tensor $\sigma$ is given by

[TABLE]

where $E$ and $\nu$ denote the Young’s Modulus and the Poisson’s ratio of the elastic medium and the deformation tensor is defined as

[TABLE]

with $d$ being the unknown displacement field.

This strategy has been successfully applied in the context of high-order curved mesh generation, see [24, 33, 26] for further details. The key aspect here is to write the approximated solution, \text{\boldmathd\unboldmath}_{h}, in separated form by using the separated representation of the imposed boundary displacement derived in Section 4.3.

The discretisation of the weak formulation associated to the strong form of the solid mechanics problem (21) leads to a system of linear equations that can be written as

[TABLE]

where $\mathbf{d}$ and $\mathbf{\delta d}$ are vectors containing the nodal values of the approximated displacement \text{\boldmathd\unboldmath}_{h} and the imposed displacement of boundary nodes respectively.

The matrix $\mathbf{A}_{11}$ in (22) is symmetric and positive definite. Thus, solving for $\mathbf{d}$ induces a linear application applied on $\mathbf{\delta d}$, see B for more details. Since $\mathbf{\delta d}$ has a separable expression, see (19), the solution of the previous system will induce a separated expression for all nodal values and, consequently, the piecewise interpolation, standard in finite elements, produces a separated representation at any point of the domain, namely

[TABLE]

which is exactly the desired structure of the displacements.

5 Numerical examples

This section presents three numerical examples that show the optimal approximation properties of the proposed PGD approach and its potential for two and three dimensional problems involving geometric parameters as extra coordinates. The examples involve the simulation of Stokes flows using a separable expression that employs the same parametric function for both velocity and pressure. This alternative was shown to be superior to other approaches in [13]. For instance, to satisfy the so-called Ladyzhenskaya-Babuška-Brezzi (LBB) condition [14]. Namely, in all the examples $k$ denotes the degree of approximation used for the velocity field and the parametric functions, whereas a degree of approximation $k-1$ is used for the pressure field. This ensures satisfaction of the LBB.

5.1 Rotating Couette flow

The first example considers the Couette flow around two infinite coaxial circular cylinders centred at the origin and with radius $R_{\text{in}}$ and $R_{\text{out}}$ respectively, with $R_{\text{in}}<R_{\text{out}}$, as represented in Figure 4. The boundary conditions correspond to known angular velocities, $\Omega_{\text{in}}$ and $\Omega_{\text{out}}$, at $R_{\text{in}}$ and $R_{\text{out}}$, respectively.

It is worth noting that the pressure must be specified at a point to remove its indeterminacy, as only Dirichlet boundary conditions for the velocity are considered. Here the pressure is imposed at one point of the outer boundary.

The analytical solution for this problem is known [7]. The azimuthal component of the velocity is given by

[TABLE]

where r=\|\text{\boldmathx\unboldmath}\|_{2}. The magnitude of the velocity of the analytical solution is also depicted in Figure 4 for $\Omega_{\text{in}}=0$ and $\Omega_{\text{out}}=1$.

The inner radius $R_{\text{in}}$ is considered an extra parameter within the proposed PGD framework and the objective is to find, in the off-line stage, the generalised velocity and pressure fields for $R_{\text{in}}\in[1,2.5]$. It is worth noting that the variation of the inner radius induces the variation of the eight control points (i.e. 16 parameters in two dimensions) of the NURBS curve describing the inner circle, as represented in Figure 4. However, as mentioned in Remark 2, the motion of these control points can be controlled by a single parameter $\mu$ representing the variation of the radius of the inner circle.

In this example the reference configuration corresponds to $R_{\text{in}}=1$ and $R_{\text{out}}=5$ and $\mu\in\mathcal{I}_{1}=[0,1.5]$. Three unstructured triangular meshes of the reference domain, with 251, 1,023 and 4,256 elements respectively, are represented in Figure 5. These meshes are generated using the technique proposed in [30] to guarantee that elements without an edge on a curved boundary can be mapped to a reference triangle using an affine mapping.

The quality of the coarsest mesh of the reference domain, measured using the scaled Jacobian [33], is represented in Figure 6.

The minimum quality observed is, as expected near the inner boundary where the elements show the maximum distortion.

Figure 7 shows the deformed configuration for three values of the parameter $\mu$, namely $\mu=0.5$, $\mu=1$ and $\mu=1.5$. These plots also represent the quality of the deformed mesh obtained by using the elastic analogy described in Section 4.4.

In all cases, the minimum quality is higher than 0.65.

To further illustrate the robustness of the mesh deformation technique employed within the proposed PGD framework, Figure 8 shows the evolution of the quality of the deformed meshes as a function of the parameter $\mu$ using four different meshes and three different degrees of approximation.

As expected, the worst case scenario corresponds to the maximum deformation induced by a parameter $\mu=1.5$, but in all cases the minimum quality is always higher than 0.6. Furthermore, it is worth noting that the influence of the parameter $\mu$ on the quality of the deformed meshes is less important for finer meshes.

A crucial aspect of the proposed PGD strategy is the separation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} defined in Equation (6). As discussed previously, \mathbf{H}_{\text{\boldmath\mu\unboldmath}} does not generally admit an exact separable expression and therefore a separable approximation is computed here via the higher-order PGD-projection [22].

Figure 9 shows the first eight normalised spatial modes of the component [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{11}.

The results suggest that the first two modes capture the main global features of [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{11}, whereas the rest of modes capture the local variations near the inner circle. The results for the second diagonal component, [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{22}, not displayed for brevity, show the same behaviour but with the expected rotation of 90 degrees due to the symmetry of the domain and the displacement field. Similarly, Figure 10 shows the first eight normalised spatial modes of the component [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{12}.

The results show, again how the first modes capture the global behaviour of the component [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{12}, whereas the last modes show relevant spatial variations in the close vicinity of the inner circle.

The first eight normalised parametric modes of \mathbf{H}_{\text{\boldmath\mu\unboldmath}} are represented in Figure 11.

It is worth recalling that the same parametric modes are associated to all the components of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}.

The amplitude, $\alpha_{m}$, corresponding to the mode $m$ of the separation of \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, is computed as the product of the Euclidean norms of the spatial and parametric functions. Figure 12 shows the amplitudes of [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{11} using three different meshes and three different degrees of approximation.

The results show that, for this example, the number of terms required to obtain a separable approximation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} using the higher-order PGD-projection is completely independent on the spatial discretisation. In all cases, 12 modes provide a decrease in the amplitude $\alpha_{m}$ of exactly 13 orders of magnitude and the amplitude is the same in all meshes and for all degrees of approximation.

Using the separation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, the rotating Couette flow problem is solved to obtain the generalised solution of the Stokes problem. The first eight normalised spatial modes of the magnitude of the velocity field are shown in Figure 13. The simulation was performed using the mesh of Figure 5 (c), 800 equally-spaced elements in $\mathcal{I}_{1}=[0,1.5]$ and a degree of approximation $k=4$.

The solution of the coaxial Couette flow problem requires imposing Dirichlet boundary conditions for the velocity at the inner and outer boundaries of the spatial domain. This is implemented within the PGD framework by adding one initial mode that fulfils the non-homogeneous Dirichlet boundary conditions. Its spatial part is computed as the solution of the Stokes problem at the reference mesh and its parametric component is taken as a constant function equal to one. The next PGD modes are computed by imposing homogeneous Dirichlet boundary conditions, ensuring that the separated PGD solution satisfies the required Dirichlet boundary conditions.

The first eight normalised parametric modes associated to the spatial modes of Figure 13 are represented in Figure 14.

It is worth recalling that the same parametric mode is associated to all the components of the velocity and the pressure fields. The results reveal that the parametric modes of the velocity field show a similar behaviour compared to the parametric modes of the separation of \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, shown in Figure 11. The first two modes are smooth whereas the next modes, that contribute less to the global solution, show a more oscillatory character.

To illustrate the gain in accuracy as the number of modes increases, Figure 15 shows the absolute value of the error in the magnitude of the velocity computed with the proposed PGD approach using $N$ modes for different values of the parameter $\mu$.

The computation has been performed on the mesh of Figure 5 (c), with 800 elements in the parametric dimension and a degree of approximation $k=4$. In all cases, the results illustrate that the computation with four modes is almost three orders of magnitude more accurate than using a single mode.

To further analyse the accuracy of the proposed PGD approach when the number of modes is increased, the relative error in the $\mathcal{L}_{2}(\Omega\times\mathcal{I}_{1})$ norm, defined as

[TABLE]

is studied. Figure 16 depicts the evolution of $\varepsilon_{\texttt{PGD}}$ as a function of the number of PGD modes for three different degrees of approximation and four different meshes.

The discontinuous lines in Figure 16 show the relative error of the standard FEM solution measured in the $\mathcal{L}_{2}(\Omega\times\mathcal{I}_{1})$ norm. The evaluation of the error of the FEM solution in the $\mathcal{L}_{2}(\Omega\times\mathcal{I}_{1})$ norm requires the solution of as many FEM problems as the number of elements in the parametric dimension used by the PGD approach, multiplied by the number of integration points within each element.

The results clearly illustrate the increased level of accuracy as the number of modes is increased. In addition, the increased accuracy induced by the use of finer meshes or higher degrees of approximation can be observed. It worth noting that the PGD solution achieves its maximum accuracy (i.e. the accuracy of a standard FEM solution) using a low number of modes. For coarse meshes, the maximum accuracy of the PGD solution is achieved with three or four modes. For finer meshes, the PGD approach provides accurate results with only three or four modes and requires between six and 11 modes to reach the same level of error as the standard FEM solution.

Next, the accuracy of the PGD solution for different values of the parameter $\mu$ is studied. Figure 17 shows the difference between the PGD solution and the FEM solution, in logarithmic scale, as a function of the parameter $\mu$ and the number of PGD modes using different degrees of approximation.

The results show that the accuracy of the PGD solution depends only weakly on the value of the parameter selected. For values of $\mu=0$, corresponding to the undeformed configuration, the error is lower for a moderate number of PGD modes. This is expected because for $\mu=0$ the displacement function $d$ is null and the Jacobian \mathbf{J}_{\!\text{\boldmath\mu\unboldmath}} is equal to the identity matrix. But for a large enough number of modes, the difference between the PGD solution and the FEM solution shows very little dependence on the value of $\mu$, illustrating the robustness of the proposed PGD approach.

The last numerical study for the rotation Couette flow considers a mesh convergence analysis for the proposed PGD approach. The PGD solution is computed on the three meshes shown in Figure 5 using three different orders of approximation. Figure 18 shows the evolution of the relative error in the $\mathcal{L}_{2}(\Omega\times\mathcal{I}_{1})$ norm, $\varepsilon_{\texttt{PGD}}$, as a function of the characteristic element size $h$.

For each simulation the minimum number of PGD modes required to achieve the maximum accuracy is considered, as described earlier when presenting the results of Figure 16. The results show that, when enough modes are considered within the PGD framework, the error $\varepsilon_{\texttt{PGD}}$ converges with near an optimal rate $h^{k+1}$.

5.2 Stokes flow around two circular cylinders

The second example is inspired by the studies of [31, 1] on the analysis of micro-swimmers. The problem involves the computation of the Stokes flow in a rectangular channel of dimension $40\times 14$ with two circular cylinders of radius $R_{1}$ and $R_{2}$, where the distance between their centres is $D$, as represented in Figure 19. Slip boundary conditions are considered on the surface of the cylinders, a free slip boundary condition on the top and bottom boundaries, a imposed horizontal velocity of unit magnitude on the left boundary and a homogeneous Neumann boundary condition on the right part of the boundary.

The geometric parameters considered are the radius of the cylinders and the distance between their centres. In addition, as done in the context of microswimmers, it is assumed that the variation of the two radius is linked so that $R_{1}^{2}+R_{2}^{2}$ remains constant, that is the area enclosed by the two cylinders is constant. It is worth noting that the variation of these geometric parameters involve the variation of all the control points (i.e. 32 parameters) of the two NURBS describing the circles represented in Figure 19. However, as done in the previous example, it is possible to re-parametrise the motion in terms of only two parameters, $\mu_{1}\in\mathcal{I}_{1}=[-1,1]$ which controls the radius of the cylinders and $\mu_{2}\in\mathcal{I}_{2}=[-1,1]$ controlling the distance between the cylinders.

The reference configuration corresponds to $R_{1}=R_{2}=0.8$ ($\mu_{1}=0$) and $D=14$ ($\mu_{2}=0$). The minimum and maximum values for the radius of the cylinders are 0.3578 ($\mu_{1}=-1$) and 1.0733 ($\mu_{1}=1$) which, given the link between $R_{1}$ and $R_{2}$, correspond to the cases where the area of one of the circles is 90% and 10% respectively of the total area occupied by both circles. For the distance, the maximum and minimum values are 15.5 ($\mu_{2}=-1$) and 12.5 ($\mu_{2}=1$) respectively.

A computational mesh with 2,338 triangular elements is generated for the reference configuration. For a cubic degree of approximation (for geometry and velocity), Figure 20 shows the quality of the mesh for the reference configuration and for two deformed configurations corresponding to the extreme cases with $\mu_{1}=\mu_{2}=-1$ and $\mu_{1}=\mu_{2}=1$.

It can be observed that, even for the large deformations corresponding to the extreme cases of Figures 20 (b) and (c), the minimum quality (measured as the scaled Jacobian) is near 0.7, which is similar to the quality observed in the previous example for a simpler problem.

Next, the separation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, defined in Equation (6), is considered. As discussed previously \mathbf{H}_{\text{\boldmath\mu\unboldmath}} does not generally admit an exact separable expression and therefore a separable approximation is computed here via the higher-order PGD-projection described in [22].

Figure 21 shows the first six normalised spatial modes of the component [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{11}.

Similarly, Figures 22 and 23 show the first six normalised spatial modes of the component [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{12} and [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{22} respectively.

It is worth noting that some of the modes in Figure 21 resemble the modes obtained in the previous example and represented in Figure 9, whereas other modes display a completely different spatial variation. This indicates that the similar modes are the ones that carry information about the change of the radius of the cylinders whereas the different ones are related to the variation of the distance between the cylinders.

Contrary to the previous example, it is apparent that more modes are necessary to describe the global behaviour of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}. For instance, in the previous example only two modes were needed to describe the global information of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} whereas now the fourth mode of all the components of \mathbf{H}_{\text{\boldmath\mu\unboldmath}} and the fifth mode of [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{22} contain global information. Also, contrary to the previous example, it is interesting to observe the different behaviour of the modes associated to the diagonal terms of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}. This is due to the more complex motion induced by the geometric parameters, compared to the previous example.

The first eight normalised parametric modes of \mathbf{H}_{\text{\boldmath\mu\unboldmath}} are represented in Figure 24.

It is worth recalling that the same parametric modes are associated to all the components of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}. It is interesting to observe the similar qualitative behaviour of the parametric modes of $\mu_{1}$, represented in Figure 24 (a), and the parametric modes obtained in the previous example, represented in Figure 11. This is expected as in both cases these modes are related to the variation of the radius of a circle. The parametric modes associated to the distance between the cylinders, represented in Figure 24 (b), display a less oscillatory character than the modes associated to the variation of the radius of a circle.

To illustrate the increased complexity due to the introduction of two geometric parameters, Figure 25 shows the amplitude, $\alpha_{m}$, corresponding to the mode $m$ of the separation of \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, computed as the product of the Euclidean norms of the spatial and parametric functions.

The amplitude, $\alpha_{m}$, corresponding to the mode $m$ of the separation of \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, is computed as the product of the Euclidean norms of the spatial and parametric functions. Figure 12 shows the amplitudes of [\mathbf{H}_{\text{\boldmath\mu\unboldmath}}]_{11} using four different meshes and three different degrees of approximation.

The results show that, for this example, the number of terms required to obtain a separable approximation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} using the higher-order PGD-projection is significantly higher than in the previous example. To provide a decrease in the amplitude $\alpha_{m}$ of six orders of magnitude the number of modes vary from 50 and 63 for an order of approximation ranging from $k=2$ to $k=4$, showing again that the number of modes is not highly dependent upon the spatial discretisation.

Using the separation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, the Stokes flow around the two cylinders is computed with the proposed approach to obtain the generalised solution. The first six normalised spatial modes of the magnitude of the velocity and the pressure fields are shown in Figure 26.

The simulation was performed using the mesh of Figure 20 (a), 25 equally-spaced elements in $\mathcal{I}_{1}=[-1,1]$, 45 equally-spaced elements in $\mathcal{I}_{2}=[-1,1]$ and a degree of approximation $k=4$.

The first eight normalised parametric modes associated to the spatial modes of Figure 26 are represented in Figure 27.

It is worth recalling that the same parametric mode is associated to all the components of the velocity and the pressure fields.

Contrary to the previous example, it can be observed that the spatial modes for the velocity and the pressure do not resemble the spatial modes of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, illustrating the increased complexity of the current problem. In addition, it is worth noting that the parametric modes associated to the first parameter $\mu_{1}$ shows a more oscillatory character than $\mu_{2}$ when the separation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} is studied (as shown in Figure 24), whereas the the second parameter $\mu_{2}$ shows a more oscillatory character than $\mu_{1}$ when the separation of the velocity and pressure is considered (as shown in Figure 27). This indicates that the flow around the two cylinders changes slowly when the radius of the cylinders is varied whereas the flow changes more rapidly when the distance between the cylinders is varied.

Next, the solutions obtained with the proposed PGD framework are represented for different values of the geometric parameters. Figure 28 shows the magnitude of the velocity and the pressure fields for the three configurations shown in Figure 20.

To illustrate the accuracy of the proposed approach, Figure 29 shows the difference between the standard FE solution and the PGD solution for the magnitude of the velocity and the pressure corresponding to the three configurations shown in Figure 20.

For all the configurations, the results obtained with the proposed PGD framework are in excellent agreement with the results obtained using a standard FE solver. For the first configuration displayed in Figure 29, corresponding to $\mu_{1}=\mu_{2}=0$, the difference between the magnitude of the velocity field using FEM and the PGD in the \mathcal{L}_{2}(\Omega^{\text{\boldmath\mu\unboldmath}}) norm is 0.0143 and the difference between the pressure field using FEM and the PGD in the \mathcal{L}_{2}(\Omega^{\text{\boldmath\mu\unboldmath}}) norm is 0.0277. For the second configuration of Figure 29, corresponding to $\mu_{1}=\mu_{2}=-1$, the difference between the magnitude of the velocity and pressure fields using FEM and the PGD in the \mathcal{L}_{2}(\Omega^{\text{\boldmath\mu\unboldmath}}) norm are 0.0267 and 0.0592 respectively. Finally, for the last configuration of Figure 29, corresponding to $\mu_{1}=\mu_{2}=1$, the difference between the magnitude of the velocity and pressure fields using FEM and the PGD in the \mathcal{L}_{2}(\Omega^{\text{\boldmath\mu\unboldmath}}) norm are 0.0218 and 0.0787 respectively. It is worth emphasising that the accurate results obtained, with differences between PGD and standard FE simulations lower than 8% for the whole range of geometric configurations, have been obtained with very coarse meshes in the parametric spaces, only 25 elements in $\mathcal{I}_{1}$ and 45 elements in $\mathcal{I}_{2}$.

5.3 Stokes flow around an arbitrarily shaped vesicle

The last example is inspired by the study of vesicles suspended in a viscous flow presented in [32]. The characterisation of such flows is of interest in many biomechanical applications and the simulations often require the computation to be performed for a large variety of geometric configurations. This example is used to demonstrate the applicability of the proposed technique in three dimensions by using three geometric parameters that lead to substantial variations in the geometric model.

The problem considered here involves the simulation of the Stokes flow around an arbitrarily shaped vesicle in a channel of dimension $10\times 10\times 20$. Slip boundary conditions are considered on the surface of the vesicle, an imposed vertical velocity of unit magnitude in the inflow part of the boundary ($z=-10$), a homogeneous Neumann boundary condition on the outflow part of the boundary ($z=10$) and free slip boundary conditions on the remaining parts of the boundary.

The generic vesicle considered is modelled using a single degenerate cubic NURBS with 20 control points and four patches, as illustrated in Figure 30.

The control points are arranged in three groups, according to their $z$ coordinate being positive, zero or negative. The geometric parametrisation considers the motion of all the three groups independently by using three parameters, $\mu_{1}$, $\mu_{2}$ and $\mu_{3}$ respectively. Each parameter is used to vary the radial position of the control points whilst maintaining its $z$ coordinate. For the control points with positive $z$ coordinate the radial coordinate is given by $r+3\mu_{1}/4$ with $\mu_{1}\in[-1,1]$. Similarly for the control points with zero or positive $z$ coordinate the radial coordinate is $r+3\mu_{2}/4$ and $r+3\mu_{3}/4$ respectively, with $\mu_{2}\in[-1,1]$ and $\mu_{3}\in[-1,1]$.

The reference configuration, corresponding to \text{\boldmath\mu\unboldmath}=(\mu_{1},\mu_{2},\mu_{3})=(0,0,0), is shown in Figure 30.

The generalised PGD solution is computed using a tetrahedral mesh with 6,712 elements, 10 equally-spaced elements in each parametric interval $\mathcal{I}_{1}=\mathcal{I}_{2}=\mathcal{I}_{3}=[-1,1]$ and a degree of approximation $k=4$.

To illustrate the variation in the geometry induced by the selected geometric parameters, Figure 31 shows six different geometric configurations using different values for the geometric parameters.

It can be observed that the variations of the control points considered lead to substantial changes in the geometric model. These variations induce a deformation of the mesh generated for the reference configuration, computed in a separated form by employing the strategy described in Section 4.4. Figure 32 shows six deformed meshes near the vesicle that correspond to the six variations of the geometric parameters depicted in Figure 31.

The colours represent the quality of the elements, measured as the scaled Jacobian.

Following the strategy described in this work, the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}} is separated first using the higher-order PGD-projection. To provide a decrease in the amplitude $\alpha_{m}$ of five orders of magnitude, 50 modes are required in this example. It is worth noting that this is almost the same number of modes required for the previous two dimensional example to provide the same degree in the amplitude of the modes of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}.

Using the separation of the matrix \mathbf{H}_{\text{\boldmath\mu\unboldmath}}, the Stokes flow around the vesicle is computed with the proposed approach to obtain the generalised solution. A total of 160 modes are computed, being the amplitude of the last mode three orders of magnitude lower than the amplitude of the first mode.

Figures 33 and 34 show the velocity and pressure fields for the six geometric configurations shown in Figure 31.

These configurations are obtained in the on-line phase as particularisations of the genearlised, six dimensional, PGD solution. Both the velocity and pressure fields illustrate the ability of the proposed PGD approach to capture significant changes in the flow field induced by geometric variations of the CAD model.

To quantify the accuracy of the proposed technique in a more complex problem in three dimensions and with three geometric parameters, the particularised solutions for the velocity and pressure are compared to standard finite element computations performed on the deformed configurations for the six cases displayed in Figures 33 and 34. For the velocity field, the difference between the PGD and the finite element solution, measured in the $\mathcal{L}_{2}(\Omega)$ norm, is 0.7%, 0.8%, 0.8%, 0.4%, 0.6% and 0.7% for the six cases shown in Figure 33 respectively. Similarly, for the pressure field, the difference between the PGD and the finite element solution, measured in the $\mathcal{L}_{2}(\Omega)$ norm, is 7%, 7%, 7%, 3%, 5% and 6% for the six cases shown in Figure 34 respectively. It is worth emphasising that the higher accuracy observed in the velocity is due to the use of a higher interpolation degree for the velocity field, compared to the pressure field, in order to satisfy the LBB conditions.

6 Concluding remarks

A computational framework for the computation of off-line solutions for a set of parameters describing the geometry of a domain has been presented. The proposed approach considers as parameters of the generalised PGD solution the control points of the CAD boundary representation of the computational domain. A mapping between a reference configuration and the current configuration is proposed by interpreting the geometric changes as a displacement field derived from a linear elastic problem. The key aspect of the proposed approach is that the displacement field is explicitly written in a separable form. This approach enables the incorporation of the PGD rationale in a high-order finite element context.

The potential of the proposed approach is shown for a variety of problems involving the solution of the Stokes equation in geometrically parametrised domains, both in two and three dimensions. The problems, of increasing difficulty show the optimal approximation properties of the method and its ability to accurately capture the flow features for significant changes of the geometric model. For the most complex problem considered in this work, the generalised solution computed with three geometric parameters in three dimensions show good agreement when compared to a standard finite element computation, with errors in the velocity field lower than 1% and errors in the pressure field below 8%.

Appendix A Integration within a CAD environment in 3D

In three dimensions, the boundary of the parameterized domain, \partial\Omega^{\text{\boldmath\mu\unboldmath}}, is assumed to be described by using a set NURBS surfaces \{\text{\boldmathS\unboldmath}_{\!j}^{\text{\boldmath\mu\unboldmath}}\}_{j=1,\dotsc,M}, being $M$ the total number of surfaces, namely

[TABLE]

Next, the necessary concepts about NURBS surfaces are briefly recalled and the minimum changes that are required to extend the technique presented in Section 4 to three dimensional domains are detailed.

A.1 NURBS surfaces

A NURBS surface of degree $q$ in $\lambda$ and degree $r$ in $\kappa$, is a piecewise rational function defined in parametric form as

[TABLE]

where \{\text{\boldmathB\unboldmath}_{il}\} are the coordinates of the $(\texttt{n}_{\texttt{cp}}+1)(\texttt{m}_{\texttt{cp}}+1)$ control points (defining the control net) and $\{R_{il}\}$ are rational basis functions defined as

[TABLE]

In the above expression $\{\nu_{il}\}$ are the control weights associated to the control points and $\{S_{i,l}^{q,r}(\lambda,\kappa)\}$ are the 2D B-spline basis functions of degree $q$ in $\lambda$ and $r$ in $\kappa$. Each 2D B-Spline basis function is defined as a tensor product of 1D basis functions, that is

[TABLE]

A.2 Geometric parameters

The geometric parameters \text{\boldmath\mu\unboldmath}\in\text{\boldmath\mathcal{I}\unboldmath}\subset\mathbb{R}^{\texttt{n}_{\texttt{pa}}} are defined as the variations of the original coordinates of the control points of the NURBS surfaces describing the boundary. More precisely, for each NURBS curve \text{\boldmathS\unboldmath}_{\!j}, with $j=1,\dotsc,M$, having $(\texttt{n}_{\texttt{cp}}^{j}+1)(\texttt{m}_{\texttt{cp}}^{j}+1)$ control points, the undisturbed boundary is characterised by the coordinates of the control points: \text{\boldmathB\unboldmath}_{il}^{j}, for $i=0,\dotsc,\texttt{n}_{\texttt{cp}}^{j}$ and $l=0,\dotsc,\texttt{m}_{\texttt{cp}}^{j}$. The boundary of the distorted domain, \partial\Omega^{\text{\boldmath\mu\unboldmath}}, is defined by the position of the displaced control points, namely \text{\boldmathB\unboldmath}_{il}^{j}+\delta\!\text{\boldmathB\unboldmath}_{il}^{j}. The displacement range for each control point is characterised by

[TABLE]

In fact, each displacement of a control point on the $j$-th NURBS, \delta\!\text{\boldmathB\unboldmath}_{il}^{j}, might depend upon the parameters and can be written as

[TABLE]

where \text{\boldmathe\unboldmath}_{i}, for $i=1,2,3$, are the unit coordinate vectors. Then \text{\boldmath\mu\unboldmath}\in\text{\boldmath\mathcal{I}\unboldmath}:=\text{\boldmath\mathcal{I}\unboldmath}^{1}\times\text{\boldmath\mathcal{I}\unboldmath}^{2}\times\dotsb\times\text{\boldmath\mathcal{I}\unboldmath}^{M}, where \text{\boldmath\mathcal{I}\unboldmath}^{j}=\text{\boldmath\mathcal{I}\unboldmath}^{j}_{1}\times\text{\boldmath\mathcal{I}\unboldmath}^{j}_{2}\times\cdots\text{\boldmath\mathcal{I}\unboldmath}^{j}_{(\texttt{n}_{\texttt{cp}}^{j}+1)(\texttt{m}_{\texttt{cp}}^{j}+1)} is the range of variation of the coordinates of the control points of the curve \text{\boldmathS\unboldmath}_{\!j}.

A.3 Separated representation of the boundary displacement

The variation of a control point \text{\boldmathB\unboldmath}_{il}^{j} of a NURBS surface \text{\boldmathS\unboldmath}_{\!j}, namely \delta\!\text{\boldmathB\unboldmath}_{il}^{j}, changes the definition of the original curve only in the support of the basis function $R_{il}^{j}$, given by the subspace of the parametric space $[\lambda_{i},\lambda_{i+q^{j}+1}]\times[\kappa_{l},\kappa_{l+r^{j}+1}]$. The modified NURBS surface is parametrised by

[TABLE]

As in the two dimensional case, the displacement of the boundary mesh node \text{\boldmathX\unboldmath}_{k}=\text{\boldmathS\unboldmath}_{\!j}(\lambda_{k},\kappa_{k}) that belongs to the NURBS curve \text{\boldmathS\unboldmath}_{\!j} can be written in separated form as

[TABLE]

where the dependence of the displacements of the control points in terms of the parameters described in Equation (26) has been used. Moreover, since the NURBS parameterd $(\lambda_{k},\kappa_{k})$ is only dependent on the spatial coordinates \text{\boldmathX\unboldmath}_{k}, and not on the geometric parameters $\mu$, the previous equation can be written as,

[TABLE]

which characterises the displacement of the boundary nodes and has the desired separated form, as in the two dimensional case given by Equation (19).

A.4 Separated representation of the geometric mapping

The strategy to obtain a separated representation of the geometric mapping is not dependent on the dimensionality of the problem. Therefore, the strategy described in Section 4.4 is also valid in three dimensions. The implementation details are given in B.

Appendix B Implementation details of the separated representation of the geometric mapping

To obtain a separated representation of the displacement function in the whole domain, the solid mechanics problem (21) is considered. Its discretisation leads to the system of equations (22).

Assuming that the mesh nodes are ordered so that the boundary nodes are first, the vector $\mathbf{D}$ is given by

[TABLE]

where $D_{k}^{X_{l}}$ is the imposed displacement of node \text{\boldmathX\unboldmath}_{k} in the $X_{l}$ direction.

As usual in a FE context, the modified system of linear equations to be solved, after accounting for the Dirichlet boundary conditions, is

[TABLE]

The solution of this system of linear equations can be written as

[TABLE]

where $\mathbf{D}_{k}=\left(D_{k}^{X_{1}},\dotsc,D_{k}^{X_{\texttt{n}_{\texttt{sd}}}}\right)^{T}$ and $\mathbf{B}^{k}$ denotes the block of the matrix $\mathbf{B}:=-\mathbf{A}_{11}^{-1}\mathbf{A}_{12}$ containing the columns from $\texttt{n}_{\texttt{sd}}(k-1)+1$ to $\texttt{n}_{\texttt{sd}}k$, with dimension $\texttt{n}_{\texttt{sd}}|\mathcal{S}|\times\texttt{n}_{\texttt{sd}}$.

Using the separated representation of the imposed displacement (19) for each boundary node, the following separated representation of the nodal values of the displacement is obtained

[TABLE]

leading to the separated representation of the approximation of the displacement function given by Equation (23).

It is worth noting that the matrix $\mathbf{B}$ in Equation (27) only depends on the spatial discretisation of the original configuration, $\Omega$, and the selected material parameters ($E$ and $\nu$) and it is independent on the geometric parameters $\mu$. Therefore, it is possible to pre-compute and store the matrix $\mathbf{B}$ so that the separated representation of the displacement of Equation (27) can be computed with a negligible cost for different imposed boundary displacements (i.e. for different configurations \Omega^{\text{\boldmath\mu\unboldmath}} induced by different variations of the geometric parameters).

It is important to recall that the dimension of the matrix $\mathbf{B}$ is $\texttt{n}_{\texttt{sd}}(\texttt{n}_{\texttt{mn}}-|\mathcal{S}|)\times\texttt{n}_{\texttt{sd}}|\mathcal{S}|$, which, in practical applications, is much lower than the size of a standard FE matrix, namely $\texttt{n}_{\texttt{sd}}\texttt{n}_{\texttt{mn}}\times\texttt{n}_{\texttt{sd}}\texttt{n}_{\texttt{mn}}$.

Acknowledgements

This work is partially supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska–Curie actions (Grant number: 675919) and the Spanish Ministry of Economy and Competitiveness (Grant number: DPI2017-85139-C2-2-R). The first author also gratefully acknowledges the financial support provided by EPSRC (Grant number: EP/P033997/1). The second and third authors are also grateful for the financial support provided by the Generalitat de Catalunya (Grant number: 2017-SGR-1471).