# Parametric solutions involving geometry integrated with computer-aided design

###### 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.

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 (with number of spatial dimensions), whose boundary is characterised by a set of geometric parameters (with number of parameters characterising the geometry) and is partitioned into Dirichlet, , and Neumann, , frontiers such that and . Note that the set , which characterises the admissible range for parameters , can be defined as the Cartesian combination of the range for each parameter, namely, with for .

For each set of parameter , the objective is to determine the parametric solution , with , of the boundary value problem

(1) |

where is the thermal conductivity (symmetric and positive definite) matrix, is a source term, is the imposed temperature, is the imposed heat normal flux and is the outward unit normal vector. The standard variational form of the previous problem reads: find for all such that

(2a) | |||

where the space of trial functions is and its corresponding test functions space is . The parametric bilinear and linear forms and are defined by | |||

(2b) |

where

denote, respectively, the product of scalar/vector functions in and its traces over .

### 2.2 The multi-dimensional parametric problem

There are different alternatives to obtain, for any given set of parameters , an approximation to the solution of problem (1). The obvious option of solving a new problem for every instance of 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 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 , but it is seen now as a function in a larger dimensional space and it is written as with . Consequently, formally lies in a tensor product space, namely, . A standard weighted residuals approach, with integrals in and the usual integration by parts only in produces a weak form in this multi-dimensional setup. Namely, find such that

(3a) | |||

with the following definitions of the bilinear and linear forms | |||

(3b) |

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 . Note that the tensor product space 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

Note that the forms 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

## 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 and , see (2b), are integrated on spacial parametrised domains (domains depending on parameters ). As suggested in [2, 34], a mapping is necessary (not sufficient) in order to have an affine parameter dependence. This mapping transforms a reference domain into the geometrically parametrised (“deformed”) domain , namely

(4) | ||||

As classically in computational mechanics, the reference configuration is associated to a reference coordinate system denoted by , whereas the distorted domain will be associated to the spatial description . Following this analogy, a displacement field is used to relate both coordinate systems.

The mapping 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 and its corresponding boundary , partitioned into the Dirichlet, , and Neumann, , boundaries, all independent of parameter . Namely,

(5) |

where is the Jacobian matrix of the mapping and, more importantly, it is the only element in the above equations that depends upon parameters . Moreover, in order to compact the notation a new matrix is introduced,

(6) |

where the definition of the adjoint operator, , 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

where , and is the kinematic viscosity. Note that denotes the component of a matrix .

### 3.1 Separated displacements

The methodology proposed here allows to obtain an affine parameter dependence (namely, a separable expression for ) 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

(7) |

which induces a separated Jacobian matrix

(8) |

where for .

### 3.2 Affine parameter dependence

The affine parameter dependence of as defined in (5), is in practice determined by obtaining a separated expression of . In order to obtain a separated expression for matrix , both and are analytically separated. Then, the Higher-Order PGD-Projection proposed in [22] is used to obtain a compact separation for .

The separated representation for the determinant, , can be obtained using Leibniz formula from equation (8), whereas for , 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 in terms of its trace and its powers, namely

(9) |

where and . In practice, the two cases of interest are 2D () and 3D () problems, which are detailed next.

### 3.3 Two-dimensional approach

In 2D the Jacobian is a matrix. From (8), the separated expression for its determinant is

(10) |

again denotes the component of the matrix .

### 3.4 Three-dimensional approach

Following the previous rationale, in 3D the Jacobian is a matrix. The separated expression for the determinant, , is obtained using Leibniz formula as

(12) |

where is the set of the six permutations of the integers , where the element in position after the reordering is denoted , and denotes the signature of (i.e. for even and for odd ). Note that such a separation will induce sums of order .

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

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

(13) |

As done in 2D, once the adjoint is separated, from (6), a separated expression for can be computed. In this case, it consists of terms. As previously noted for 2D, the separated expression of is obtained with the Higher-Order PGD-Projection described next.

### 3.5 The Higher-Order PGD-Projection

Finally, to obtain a separated approximation of , namely

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 projection, say

for all 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 is obtained by solving

(14) |

where now, as usual in PGD, 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 ) 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, , is assumed to be described by a set of NURBS curves , being the total number of curves, namely

Next, the necessary concepts about NURBS curves are briefly recalled and the proposed strategy to build a geometric mapping that can be written in the separated form (7) is presented in detail.

### 4.1 NURBS curves

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

(15) |

where are the coordinates of the control points (forming the control polygon) and are rational basis functions defined as

In the above expression are the control weights associated to the control points and are the normalized B-spline basis functions of degree , which are defined recursively by

for and where (for ) are the knots or breakpoints, which are assumed ordered . They form the so-called knot vector,

which uniquely describes the B-spline basis functions. The number of control points, , and knots, , are related to the degree of the parametrisation, , by the relation , see [25] for more details.

Figure 1 shows an example of a two dimensional domain where the boundary is described by five NURBS curves.

### 4.2 Geometric parameters

The geometric parameters 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 , with , having control points, the undisturbed boundary is characterised by the coordinates of the control points: , for . This configuration will be used as the reference one in and will be associated to a reference coordinate system denoted by . The distorted domain, , will be associated to the spatial description . The boundary in the spatial domain, , is defined by the position of the displaced control points, namely . The displacement range for each control point is characterised by

In fact, each displacement of a control point on the -th NURBS, , might depend upon the parameters and can be written as

(16) |

where , for , are the unit coordinate vectors. Then , where is the range of variation of the coordinates of the control points of the curve .

Consequently, the maximum number of geometric parameters is

(17) |

but in practical problems, the number of geometric parameters is drastically lower than 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 of a NURBS curve , namely , changes the definition of the original curve only in the support of the basis function , given by the subspace of the parametric space . The modified NURBS curve is parametrised by

(18) |

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 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 associated to the control point .

Given a computational mesh for the reference configuration, , 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 , the index and parametric coordinate such that are computed. Then, for a deformed configuration, induced by the variation of a control point of a NURBS curve , namely , the new position of each boundary node is computed as .

Therefore, the variation of the control points of a boundary curve induces a displacement of the boundary mesh nodes, namely

where is the index set of mesh nodes on the boundary of the computational domain and 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 that belongs to the NURBS curve can be written in separated form as

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 is only dependent on the spatial coordinates , and not on the geometric parameters , the previous equation can be written as,

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

(19) |

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

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 , 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 of a NURBS curve can be defined as

(20) |

### 4.4 Separated representation of the geometric mapping

The proposed strategy to build a mapping between the reference configuration, , and the current configuration, , 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 is considered

(21) |

where is an external force defined by the user. The stress tensor is given by

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

with 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, , 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

(22) |

where and are vectors containing the nodal values of the approximated displacement and the imposed displacement of boundary nodes respectively.

The matrix in (22) is symmetric and positive definite. Thus, solving for induces a linear application applied on , see B for more details. Since 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

(23) |

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 denotes the degree of approximation used for the velocity field and the parametric functions, whereas a degree of approximation 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 and respectively, with , as represented in Figure 4. The boundary conditions correspond to known angular velocities, and , at and , 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

(24) |

where . The magnitude of the velocity of the analytical solution is also depicted in Figure 4 for and .

The inner radius 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 . 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 representing the variation of the radius of the inner circle.

In this example the reference configuration corresponds to and and . 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 , namely , and . 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 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 , but in all cases the minimum quality is always higher than 0.6. Furthermore, it is worth noting that the influence of the parameter 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 defined in Equation (6). As discussed previously, 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 .

The results suggest that the first two modes capture the main global features of , whereas the rest of modes capture the local variations near the inner circle. The results for the second diagonal component, , 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 .

The results show, again how the first modes capture the global behaviour of the component , whereas the last modes show relevant spatial variations in the close vicinity of the inner circle.

The first eight normalised parametric modes of are represented in Figure 11.

It is worth recalling that the same parametric modes are associated to all the components of the matrix .

The amplitude, , corresponding to the mode of the separation of , is computed as the product of the Euclidean norms of the spatial and parametric functions. Figure 12 shows the amplitudes of 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 using the higher-order PGD-projection is completely independent on the spatial discretisation. In all cases, 12 modes provide a decrease in the amplitude 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 , 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 and a degree of approximation .

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 , 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 modes for different values of the parameter .