Isogeometric analysis-based reduced order modelling for incompressible linear viscous flows in parametrized shapes

In this work we provide a combination of isogeometric analysis with reduced order modelling techniques, based on proper orthogonal decomposition, to guarantee computational reduction for the numerical model, and with free-form deformation, for versatile geometrical parametrization. We apply it to computational fluid dynamics problems considering a Stokes flow model. The proposed reduced order model combines efficient shape deformation and accurate and stable velocity and pressure approximation for incompressible viscous flows, computed with a reduced order method. Efficient offline-online computational decomposition is guaranteed in view of repetitive calculations for parametric design and optimization problems. Numerical test cases show the efficiency and accuracy of the proposed reduced order model. 1 Focus and motivation The capability to perform fast simulations is becoming increasingly relevant for several applications in engineering sciences, related for instance to naval and aeronautical engineering, as well as biomedicine. To this end, reduced basis methods [1, 2], proper orthogonal decomposition [3, 4, 5], proper generalized decomposition [6, 7], hierarchical model reduction [8, 9, 10], or more in general reduced order modelling (ROM) techniques [11], have received considerable attention in the last decades. ROMs do not replace, but rather build upon as an add-on, high-fidelity methods such as finite element, finite volume or discontinuous Galerkin methods. Indeed, the choice of the high-fidelity solver can be made depending on the particular problem at hand and on pre-existing expertise and software availability. Current literature has explored a broad variety of options, including reduced order models based on a finite element high-fidelity discretization (e.g. [12, 13, 14, 2, 15]), finite volume (e.g. [16, 17, 18, 19]) and finite difference methods (e.g. [20, 21, 22]). More recently, investigations towards the coupling with discontinuous Galerkin methods for multiscale problems [23] or domain-decomposition approaches [24, 25, 26], spectral element methods [27, 28], and extended finite element methods [29, 30] have been carried out. The aim of this work is to embed isogeometric analysis (IGA) [31, 32] as a high-fidelity discretization option in a ROM setting, for the simulation of incompressible linear viscous flows [33, 34, 35, 36] and to propose a complete workflow (pipeline) integrated with Free From Deformation (FFD) as efficient geometrical parametrisation. The latter is enhanced into an IGA context ready to be used within reduced order method (POD). A considerable advantage of IGA with respect to classical finite element analysis is the possibility to avoid any


Focus and motivation
The capability to perform fast simulations is becoming increasingly relevant for several applications in engineering sciences, related for instance to naval and aeronautical engineering, as well as biomedicine. To this end, reduced basis methods [1,2], proper orthogonal decomposition [3,4,5], proper generalized decomposition [6,7], hierarchical model reduction [8,9,10], or more in general reduced order modelling (ROM) techniques [11], have received considerable attention in the last decades. ROMs do not replace, but rather build upon as an add-on, high-fidelity methods such as finite element, finite volume or discontinuous Galerkin methods. Indeed, the choice of the high-fidelity solver can be made depending on the particular problem at hand and on pre-existing expertise and software availability. Current literature has explored a broad variety of options, including reduced order models based on a finite element high-fidelity discretization (e.g. [12,13,14,2,15]), finite volume (e.g. [16,17,18,19]) and finite difference methods (e.g. [20,21,22]). More recently, investigations towards the coupling with discontinuous Galerkin methods for multiscale problems [23] or domain-decomposition approaches [24,25,26], spectral element methods [27,28], and extended finite element methods [29,30] have been carried out.
The aim of this work is to embed isogeometric analysis (IGA) [31,32] as a high-fidelity discretization option in a ROM setting, for the simulation of incompressible linear viscous flows [33,34,35,36] and to propose a complete workflow (pipeline) integrated with Free From Deformation (FFD) as efficient geometrical parametrisation. The latter is enhanced into an IGA context ready to be used within reduced order method (POD). A considerable advantage of IGA with respect to classical finite element analysis is the possibility to avoid any geometrical approximation error and to perform direct design-to-analysis simulations by replacing classical mesh generation, and employing the same class of functions used for geometry parameterization in CAD packages during the analysis process. Even though most modern CAD tools are based on boundary representation (B-Rep) objects, it is still possible to use them in three-dimensional isogeometric analysis, by extending the computational domain inside (or outside) the enclosing (or enclosed) CAD surface (see, for example, [37]). A robust and reliable solution for such passage is still lacking, making this step an open question. However, the superior approximation properties of IGA methods make their adoption appealing also in biomedical and bioengineering applications [38], notwithstanding the fact that in this case the geometry is normally obtained through an approximate NURBS reconstruction of medical images.
Once the three-dimensional tensor product representation of the geometry is available, there is no distinction in computational cost or implementation complexity, with respect to simulations done on elementary geometries.
Preliminary related IGA-ROMs have been applied to steady potential flows [39,40], parabolic problems [41] or shell structural models [42]. In this work offline-online IGA-ROM is applied for the development of stable computational reduction strategies for viscous flows problems in parametrized shapes by FFD means. We investigate IGA-ROMs in a different context with respect to earlier works [39,40]. In [40] the authors neglect viscous terms and formulate the highfidelity discretization in terms of boundary integral equations and boundary element methods (BEM) to study external flows. The main novelty of the present work, besides the investigation of the other side of the spectrum of incompressible regimes (that is, when the Reynolds number tends to zero), is the coupling of FFD techniques applied to IGA geometries, for internal flows, and using finite element based IGA, in view of studies dealing with nonlinear viscous flows, for which BEM is not suited.
We would like to remark here that, although the background idea is the same as the one presented in [40], several technical issues are fundamentally different. One of the most obvious one is that the discrete systems obtained through boundary integral formulations are in general full, which implies that higher order and higher continuity finite element spaces do not influence the bandwidth of the resulting matrix. In finite element formulations of IGA methods, however, this is an important issue, and it may result in reduced performances also of the final reduced order model. In this work we show how the increased bandwidth of the high fidelity solver does not influence negatively on the combination IGA-ROM, provided that stable approximations are used for the high fidelity solver.
The proposed integrated approach is composed of the following numerical techniques: (i) isogeometric analysis, that integrates the geometrical representation of the domain and the finite dimensional approximation of the fluid dynamics problem [32], (ii) free-form deformation to efficiently deform the computational domain by means of few geometrical parameters [43], and (iii) proper orthogonal decomposition-based reduced order modelling to generate a stable reduced basis to be queried to cut down the computational cost of numerical simulations [44]. This integration has been introduced in a preliminary version in [45].
The approach we present is completely integrated and automatic from CAD to simulation, taking advantage of IGA and FFD perspectives for the accurate and efficient management of parametrized domains and shapes. The split between offline and online computational steps is crucial and it allows the versatility of bringing this proposed computational approach on very different devices, scenarios and situations in design and optimization, for instance.
The structure of the work is as follows. The parametrized formulation and the IGA method are introduced in Section 2; necessary assumptions related to the offline-online decomposition are also summarized. Section 3 summarizes the free-form deformation map which is employed to prescribe geometrical variations. The proposed stable POD-Galerkin ROM is introduced in Section 4, and 2D and 3D numerical tests are performed in Section 5 into an optimisation framework. Finally, conclusions and perspectives follow in Section 6.
2 Problem formulation and isogeometric analysisbased high-fidelity approximation

Parametrized formulation
The problem of interest throughout this work is a parametrized incompressible steady Stokes problem, obtained as a simplification of Navier-Stokes equations when inertial forces can be neglected, compared to viscous forces. Parameters of interest, denoted by µ ∈ D ⊂ R G , are related to the geometrical representation of the domain Ω = Ω(µ) ⊂ R d . The parametrized Stokes problem reads: find (u(µ), p(µ)) such that with boundary conditions representing essential and natural boundary conditions for Stokes equations, respectively. Here ν is a constant kinematic viscosity, while f (µ), g and h are prescribed forcing terms, boundary velocity profiles, and boundary tractions, respectively. For simplicity we assume that the sections Γ D and Γ N do not depend on the geometrical parameters, while the remaining part of the boundary Γ W (µ) = ∂Ω(µ) \ (Γ D ∪ Γ N ) may depend on µ.
Isogeometric formulations of Stokes flows have been extensively studied in the literature. We refer to [46] for a comprehensive analysis of stable choices of isogeometric finite element spaces, and to [47] for an alternative formulation based on boundary integral equations.

Isogeometric description of the parametrized domain
A CAD representation of the domain is usually obtained through a set of control points {P i } Ng i=1 , where in general P i ∈ R d is a d-dimensional IGA control point a , whose position depends on the geometrical parameters µ.
A d-dimensional geometrical representation is obtained by tensor product of d one-dimensional B-spline basis functions, denoted by ξ d i (s) , ξ d i : [0, 1] → R, and defined recursively as where is the d-th knot vector, a non-decreasing set of coordinates in the s parameter space, whereas p d is the polynomial order of the basis functions along the direction d.
Multivariate B-spline basis functions in R d (see for example Figure 1) can then be defined by tensor product as For simplicity of exposition we work on single patch geometries, where the reference domain is [0, 1] d , and we refer to [32,48] and the references therein for possible generalizations to multipatch geometries. The reference domain Ω = [0, 1] d can be deformed into the computational (parametrized) domain Ω(µ), by introducing a parametrized map c(·; µ) : Ω → R d : that depends on the parameter vector µ through the set of N g IGA control , where the subscript g indicates geometry. Different parameter values will produce different IGA control points and, thus, different computational domains. We will characterize in Section 3 how to efficiently prescribe the dependence of P i on µ to obtain a broad range of admissible shapes. This parametrization is crucial to embed IGA in a ROM setting dealing with parametric shapes.

Weak formulation on the reference domain and discrete problem
In order to derive a discrete approximation of the parametrized Stokes problem (1)-(2), we introduce a weak formulation on the reference domain Ω. For simplicity of exposition, we will use the same notation we used in equation (1) for the velocity and pressure fields, even though here the domain is different. Denote by V = [H 1 (Ω)] d and Q = L 2 (Ω) the velocity and pressure spaces. Multiplying (1) by test functions v • c and q • c (for the velocity and pressure field, respectively), integrating by parts and pulling back to the reference domain, we obtain the following problem: given µ ∈ D, find u ∈ V and p ∈ Q such that where the bilinear forms appearing in (5) are: Here, J (µ) is the Jacobian of the mapping c(s; µ). The linear form F (v; µ) encodes forcing terms, essential boundary conditions (by divergence-free lifting) and natural boundary conditions. Isogeometric approximations of Stokes flows (see, for example [46]) violate somehow the isogeometric paradigm, in the sense that we require two different B-spline spaces for the velocity and pressure fields in order to satisfy the infsup condition, and only one of the two is usually taken to be the same as the geometrical B-spline space. We introduce V N ⊂ V and Q N ⊂ Q, of dimensions N u and N p respectively. To differentiate w.r.t. to the geometric basis functions (which are always taken to be scalar, since we encode the dimensional information in the control points), we use the following, more general, notation respectively. An alternative notation, that allows one to distinguish between the properties of the different isogeometric spaces (see, for example, [36,46]) is given by the following: (8) where p i and α i represent respectively the degree and the maximal regularity in the i-th direction.
If one chooses to use the same basis functions for the geometry and the velocity (for example), then φ i are vector versions of B i , and N u = dN g , where N g is the number of the geometry basis functions. For an extensive discussion on the choices of stable pairs of isogeometric finite element approximations of Stokes flows, we refer the reader to [46] and the references therein. In this work we used a Taylor-Hood approximation (as presented, for example, in [36]), in which the pressure space is taken to be one degree less of the velocity space, maintaining the same knot vectors of the geometry and velocity spaces, i.e., we consider pairs of spaces given by (S p ,...,p p−1,...,p−1 − S p−1,...,p−1 p−2,...,p−2 ) which satisfy the inf-sup condition and represent a good balance between attainable accuracy and computational efficiency.
The isogeometric Galerkin formulation of the problem becomes: given µ ∈ D, find u N ∈ V N and p N ∈ Q N such that where u N = u N (µ) ∈ V N and p N = p N (µ) ∈ Q N denote the high-fidelity velocity and pressure solutions, respectively. Equation (9) can be written in matrix form as where and we indicate with u(µ) and p(µ) the R Nu and R Np vector of coefficients of the discrete, high-fidelity, velocity and pressure fields respectively.

Affine parametric dependence assumption
In this work we seek for an offline-online decomposition of the computational stages, as required in the reduced order modelling context for an efficient evaluation of the ROM [1]. During the offline stage, which we will summarize in Section 4.1, we carry out all expensive computations (related to the IGA highfidelity model); in contrast, we look for an online phase (related to the ROM) which is extremely fast (see Section 4.2). In order to achieve this, we require that matrices and vectors in (11) fullfil the following affine parametric dependence assumption: We employ the empirical interpolation method (EIM) [49] to approximate this assumption up to a desired tolerance. See also [50,51,52,53] for the application of EIM to viscous flows in parametrized domains.

Shape parametrization by free-form deformation
In this section we show how to relate geometrical parameters µ to the IGA control points position P i (µ). Unfortunately, choosing the IGA control points position as geometrical parameters (i.e. G = dN g and [P i (µ)] j = µ (i−1)d+j , i = 1, . . . , N g , j = 1, . . . , d) results in an extremely high parameter space dimension G 1 which, in turn, may lead to poor performance of the reduced order model (e.g. due to an intractable number of terms in the affine expansions (12)). The aim of this section is to introduce an efficient representation of the deformation of parametrized domains described by the IGA transformation (4).

Free-form deformation map
Free-form deformation (FFD) techniques, introduced in [43] in the late 80s, are a powerful tool for the deformation of a computational domain by means of a small number of displacements. FFD maps have been employed in the reduced order modelling framework for the first time in [54], as well as applied to shape optimization problems in [55], in both cases considering an underlying finite element high-fidelity discretization. FFD has been exploited in [54,55] to handle the deformation of Ω into Ω(µ) as the result of the application of the FFD map to each node of the finite element mesh. In contrast, in this work, we apply FFD to IGA control points to obtain their deformed position , and then rely on the map c(s; ·) in (4) to describe the deformed domain Ω(µ). To further highlight the sequential nature between the highfidelity IGA spatial description and the application of FFD map to its control points we will follow the original derivation in [43], that uses a different set of basis functions (Bernstein polynomials) than the more general ones employed in Section 2. In any case, further extensions to B-splines or NURBS can also be pursued [56].
Denote by D ⊂ R d a box that contains all IGA control points obtained (e.g.) for µ = 0. Moreover, in order to apply Bernstein polynomials defined on the reference hypercube b D = [0, 1] d , let ψ(p) be the affine function that maps D to D. A (second) set of equispaced control points {Q j } Ng j=1 , namely the FFD control points is introduced, where N g := d k=1 N g,k being N g,k the number of FFD control points in the coordinate direction k. The deformed position of the j-th control point is then obtained as Q j + µ j . Since it is possible for some FFD control points to be fixed or to be allowed to move only in some prescribed coordinate direction, the parameter vector µ ∈ R G will contain only the non-zero displacement components, so that G ≤ dN g . Effective computational reduction is obtained if N g N g ; numerical tests will show that only a small number of FFD control points will be necessary to obtain a large range of admissible shapes.
The Free-Form Deformation map T (·; µ) : D → R d is defined as the composition and b j (p) is the tensor product of one-dimensional Bernstein polynomials Finally, the parametrized position of each IGA control point is obtained applying the FFD as follows:

More practical geometrical parameters in channel configurations
One of the drawbacks of FFD from practical point of view is the lack of immediate interpretation of its parameters. Indeed, FFD is not interpolatory, so the magnitude of the displacement of a control point is not exactly equal to the actual deformation obtained at that spatial location. Recent works have improved the versatility (from the user point of view) of complex shape parametrization techniques thanks to the automatic prescription of control points position based on more intuitive geometrical parameters [57,58]. In particular, for the test cases of Sections 5.2 and 5.3, we take advantage of similiar ideas to propose a meanline FFD based on two (and four) intuitive geometrical parameters related to two (and four) admissible rotations of the meanline of a channel configuration. For the case of two rotations, a summary of the meanline FFD is shown in Figure 2. For the four rotation case, the extension is straightforward. Starting from a reference mesh Ω and associated IGA control points figure  (a)), a bounding box D and a lattice of FFD control points {Q j } Ng j=1 are introduced (figure (b)). The reference meanline of D is divided in four intervals (figure (c)). In particular, in view of obtaining different channel configurations, we employ two intuitive geometrical parameters θ 1 and θ 2 related to the rotation of the second and third interval ( figure (d)). FFD geometrical parameters {µ j } Ng j=1 are then automatically updated, being zero for all FFD control points in the first and last section, and rotated by θ 1 (θ 2 , respectively) in the second (third, respectively) section, as shown in figure (e). Finally, the position of IGA control points {P i (µ)} Ng i=1 is updated and (4) is applied to get the deformed domain Ω(µ) (see figure (f)). Since the relation between (θ 1 , θ 2 ) and the FFD parameters {µ j } Ng j=1 can be automatically obtained in Sections 5.2 and 5.3 we will refer to the former as geometrical parameters. Nevertheless, for the sake of exposition in the next section we still mantain the more general notation µ to denote them.
In a similar way, FFD can also be employed to perform local variations to the section area. In particular, FFD is applied in Figure 3 to enlarge the outlet section of 3D channel with rectangular section. This requires one geometrical parameter in 2D (related to the height of the outlet section) and two geometrical parameters in 3D (related to the width and height of the outlet section).

A POD-Galerkin ROM for parametrized Stokes equations
In this section we summarize a reduced order model (ROM) for parametrized Stokes equations based on a Proper Orthogonal Decomposition method and a Galerkin projection (see [59] for a deeper insight in the subject).

Reduced basis construction through Proper Orthogonal Decomposition
In the offline stage, denote by Ξ train = {µ 1 , . . . , µ Ntrain } ⊂ D a (usually large) training set of N train points. For each sample point µ i the high-fidelity IGA solver is queried to obtain truth velocity and pressure solution. The following snapshot matrices are then considered   A POD basis for the velocity and pressure reduced spaces are then obtained by a thin singular value decomposition (SVD) of the snapshot matrices, i.e. respectively) is the matrix representing the velocity (pressure, respectively) inner product; • U u ∈ R Nu×Ntrain (U p ∈ R Np×Ntrain , respectively) contains the velocity (pressure, respectively) left singular vectors of S u (S p , respectively); • W u ∈ R Ntrain×Ntrain (W p ∈ R Ntrain×Ntrain , respectively) is an orthogonal matrices of the velocity (pressure, respectively) right singular vectors of S u (S p , respectively); • Σ u ∈ R Ntrain×Ntrain (Σ p ∈ R Ntrain×Ntrain , respectively) is a diagonal matrix, containing the singular values of S u (S p , respectively) sorted in descending order.
Moreover, the so-called supremizer enrichment is employed in this work in order to satisfy the inf-sup stability also at the reduced order level [59,60,61]. Thus, for each training sample the following elliptic problem is solved The resulting supremizer snapshots s(µ i ), i = 1, . . . , N train are then stored in a snapshot matrix S s , on which a thin SVD is performed as described previously.
Finally, the reduced spaces dimensions N u are chosen such that the retained energy I u , given by the sum of the squares of the singular values up to N u normalized by the sum up to N train , is larger than a prescribed treshold. A similar procedure is applied to choose N s and N p . The basis functions of the reduced velocity space V N are then obtained as the union of the first N u left singular vectors of X 1/2 u S u to the first N s left singular vectors of X 1/2 u S s . Similarly, the basis functions of the reduced pressure space Q N are given by the first N p left singular vectors of X 1/2 p S p . The corresponding basis function matrices, that hold the basis functions as column vectors, are denoted by Z u,s and Z p , respectively.

Reduced order approximation through Galerkin projection on the reduced spaces
In the online stage, we let µ ∈ D be a new value and we seek an approximation of the form through a Galerkin projection over the reduced spaces V N and Q N . Therefore, the following problem has to be solved: where (see e.g. [1] for a detailed description) and we indicate with u N (µ) and p N (µ) the R Nu+Ns and R Np vectors of coefficients of the reduced order approximation of velocity and pressure fields. Moreover, thanks to the affine dependence assumption (12), during the online stage each block the ROM linear system (14) can be assembled as where the following matrices have been built at the end of the offline stage and stored in memory: u,s f q , resulting in very efficient (N independent) online queries. We refer to Figure 4 for a summary of the proposed reduced order model, where the offline stage is shown in red, while the online phase is displayed in green.
The exact solution is imposed as a Dirichlet boundary condition c at ∂Ω. In Figure 5 we plot the convergence test for the solution over several refinement cycles on a uniform grid. The rate of convergence is the one predicted by an a priori analysis, as shown in [46]. In Figure 6 the numerical solution for the last iteration is shown.
As a second test, we consider the two-dimensional Poiseuille flow in a rectangular channel Ω = [0, L] × [−l, l] (see Figure 7), whose solution is analytical:  In this case, the error on the numerical solution reaches the machine epsilon already for a single IGA element, that is, for 18 DoFs for u and 4 for p. This behaviour is related to the fact that the solution is quadratic in the velocity and linear in the pressure and the fact that we are using (S 2,2 1,1 − S 1,1 0,0 ) for the solution of the problem for both the preliminary tests. In Figure 8 we plot the solution for Poiseuille flow on a rectangular domain with L = 10 and l = 0.5.

Reduced order approximation of Poiseuille-like flows with meanline FFD
Once the code for the Poiseuille flow has been validated, we keep the same model and boundary conditions and deform the original rectangle (for the two dimensional problem) or parallelepiped (for the three dimensional case) domain through FFD, obtaining a family of possible different configuration of Poiseuillelike flows, such as the one depicted in Figures 10 and 16, and provide main results regarding the ROM framework explained in Section 4. A summary of the computational details is given in Table 1. In Figure 9 we provide the geometry of the four problems we treat during the model order reduction: problem 1 characterized by two rotations and 2D (S 2,2 1,1 −S 1,1 0,0 ) elements, problem 2 featuring two rotation problem and an approximation by high-order     In Figure 13 we perform an error analysis on the solution of the reduced model compared to the high-fidelity one for the geometry configuration of problem 1. In particular, they show that 10 basis functions are enough to have an error lower than 10 −3 for both pressure and velocity. For the sake of visualization, we also report the reconstructed velocity and pressure fields in Figure 11 obtained for 10 basis functions. We can compare it with the visualization of the high-fidelity solution of Figure 10. Similar considerations apply for the other problems: see Figure 15 for problem 2, Figures 19, 16 and 17 for problem 3, and Figure 23 for problem 4. Table 1 also highlights several factors that slightly affect the online performance in terms of CPU time. A first point to take into account is related to the number of terms resulting from the EIM approximation of parametrized tensors: comparing problems 1 to 2 and 3 we can see that both an increased high-fidelity discretization order and an higher number of parameters result in a larger number of EIM terms. A second factor to take into account is related to the reduced space dimension. This can be observed comparing problems 1 and 4, where the latter requires a larger reduced space due to a slower decay of POD singular values. In any case, computational speedups are of at least an order of magnitude. Moreover, problem 4 is characterized by a speedup of order 10 2 .

Shape optimization of Poiseuille-like flows with ROM and meanline FFD
We now present the results of the shape optimization routine for the deformable pipe. Motivated by the error analysis of the previous section, we choose N = 10.
The aim is to find the parameter values that minimize the pressure drop in the pipe, for prescribed inflow section and parametrized meanline variation (and outlet section, in case of problem 4). For prescribed outlet section, the exact result of the optimization procedure is the straight pipe, obtained for null value of the angles; for parametrized outlet section, the exact solution is characterized by null angles and maximum outlet area. The optimal parameter is denoted by µ * . Details about the optimization algorithm are summarized in Table 2. In Table 3 we summarize the main results for the optimization process, both for the high fidelity solver and for the reduced order model. The error on the angles and on the pressure drop is negligible in the case of the high fidelity solver. The error for the ROM is of the order of 10 −4 (10 −4 , respectively), and we obtain a computational speedup of about 36, for the two rotation case. Interestingly, such speedup is considerably higher than the speedup for a single simulation (which is around 20), most likely because it is generally easier for optimization software to explore a smaller state space, and some smarter procedure may be used internally to save computational effort. This behaviour is less evident for the four rotation case (problem 3). We expect that also in the nonlinear case the computational speedup would increase more considerably.
This simple shape optimization test case highlights the capability of the   In future more complex applications will deal with the optimal design process of aero-hydrodynamic components.

Conclusions and future work
We have presented a complete parametric design pipeline from CAD to accurate and efficient numerical simulation, by introducing geometrical parametrization based on FFD, high order simulations based on IGA and efficient and stable computational reduction strategies based on proper orthogonal decomposition, after the enrichment of the velocity space with suited supremizers. This setting is motivated and developed by industrial applications in mechanical, nautical and naval engineering at low Reynolds number (e.g. microfluidics devices characterized by low velocity flows and in small geometrical configurations). Results look promising to continue with the implementation of a viscous nonlinear model and more complex physical and geometrical problems in order to deal with more advanced fluid mechanics indexes (vorticity, viscous stresses, viscous energy dissipation), derived from the state equations. For example, we mention the project UBE (Underwater Blue Efficiency) whose goal is the shape optimization of immersed parts of motor yachts, including exhaust flow devices, for the reduction of emissions and vibrations , in order to increase on-board comfort. This parametric design automatic embedded pipeline is motivating also the investigation and improvement of some computational aspects related with FFD and the already mentioned EIM.