Finite transformation rigid motion mesh morpher; A grid deformation approach

In any shape optimization framework and specifically in the context of computational fluid dynamics, a robust and reliable grid deformation tool is necessary to undertake the adaptation of the computational mesh to the updated boundaries at each optimization cycle. Grid deformation has its share of challenges, namely, to maintain high mesh quality (avoid distorted elements and tangles) even when dealing with extreme deformations. In this work a novel grid deformation algorithm, the finite transformation rigid motion mesh morpher (FT‐R3M) is proposed. FT‐R3M is essentially a mesh‐free grid deformation approach, since it does not require any inertial quantities and it gracefully propagates the movement of the boundaries (surface mesh) to the internal nodes of the mesh (volume mesh), by keeping the motion of its parts (referred to as stencils) as‐rigid‐as‐possible. It is an optimization‐based method, which means that the interior nodes of the computational mesh are displaced to minimize a distortion metric related to the elastic deformation energy, by favoring rigidity in critical directions, thus being able to handle mesh anisotropies very efficiently. Results are presented for three test cases; a rotated airfoil with a mesh appropriate for viscous flow; a simulation of a low Reynolds duct case; a beam.


INTRODUCTION
A variety of stochastic and gradient-based methods have been devised to solve aerodynamic shape optimization problems such as the design of ducts or manifolds for minimum power losses, cars with optimal combination of drag and lift, 1 aircraft/wings for maximum lift 2,3 or minimum drag, 4 and so forth.During the optimization, a technique to deal with the necessary changes of the boundaries (surface mesh), namely, to adapt the computational mesh (volume mesh) to the updated geometry, is required.
A well-known method to handle the boundary changes during optimization, is remeshing, in which after every optimization cycle a computational mesh is generated anew.This is a time-consuming process and especially in complex geometries manual intervention might be needed, something which is not desired in a fully automated optimization process.In addition, in the case of a gradient-based optimization method, inconsistencies are produced in the gradient since the latter is computed at isoconnectivity, which makes this technique neither robust nor reliable.
On the other hand, a rival and a promising way to deal with the necessary changes is to use a grid deformation tool (a.k.a.mesh morpher) which propagates the movement of the boundaries to the interior computational mesh.The prerequisite is to generate a mesh only once, at the beginning of the optimization process and then, the mesh morpher will undertake the adaptation of the computational mesh to the updated boundary shapes at each optimization step without changing the mesh topology.
The challenge when a grid deformation tool is used, is to maintain a good mesh quality 5,6 of the deformed mesh, namely, to avoid tangled or distorted elements, highly skewed cells, changes in the aspect ratio of the elements, and so forth, which may cause nonconvergence or even divergence of the computational fluid dynamics (CFD) software, for example, in CFD-based problems.][9][10][11] Most of them require a trade-off between the attained mesh quality and CPU cost.A popular grid deformation method is the spring analogy method. 12,13In this approach, each edge of the mesh is "replaced" by a tension spring with the spring stiffness as the inverse of the edge length.Unfortunately, this method does not prevent inverted elements and it fails when the local mesh motion exceeds significantly the local mesh size.Many improvements have been achieved to prevent element inversion by introducing torsional springs 14 or semitorsional springs, 15 but handling mesh anisotropies is still a challenge.Another simple and easy way to deform a mesh, in terms of implementation complexity, is by using the Laplacian smoothing, [16][17][18] or the so-called Laplacian coordinates.The computation of the displacement of the internal nodes of the computational mesh requires the solution of a linear system the dimension of which scales with the number of nodes of the mesh.This technique proves to be efficient but not very robust in case of significant rotations and complex transformations.
A competitive grid deformation tool is the linear elasticity morpher, in which the computational mesh is handled as an elastic solid body. 19,20In this method, mesh deformation is accomplished by solving the linear elasticity equations for each mesh point inside the geometric domain.Since the elasticity equations contain material properties (Young's modulus and Poisson's ratio), these should be related to the mesh characteristics.In the context of a finite volume CFD solver, implementing this method is not trivial, since it is usually solved with the finite element method.A very promising and well-established grid deformation technique is this one based on the fundamental work of de Boer et al., 21 Rendall and Allen, 22 and Jakobson and Amoignon 23 on mesh morphing based on radial basis functions (RBF) interpolation.Mesh morphing based on RBF is a robust and reliable grid deformation technique but computationally very heavy, as the matrices involved in the computations are very dense.The cost of solving such a system with an iterative method scales with the number of surface points (N s ) that cause the deformation as O(N 2 s ).Improvements have been achieved by using appropriate preconditioners 24 and the fast multipole method 25,26 for the acceleration of the matrix-vector multiplications that can reduce the cost to O(N s logN s ), or even less, but at the cost of increasing the complexity.The work of Rendall and Allen 27,28 is based on data reduction algorithms for selecting the surface points in such an efficient way so that the cost of training the RBF network is reduced up to eight times in comparison to previous works and the cost of the mesh update up to 55 times.Their work has been followed by Michler,29 where he is focused on the deflection of control surfaces of aircraft configurations where the meshes are very dense and big.In his work, locally confined RBF-based mesh deformation has been applied, the mesh deformation is restricted in the vicinity of the deflecting control surfaces, leading thus to a reduced computational cost of the method.A big challenge in this technique is the selection of nonuniform spacing of points and handling of mesh anisotropies.
The rigid motion mesh morpher has already been introduced. 30It has been shown that it can handle mesh rotation and mesh anisotropy very efficiently.An improvement upon this concept, is the development of the finite transformation rigid motion mesh morpher (FT-R3M), which is presented in this article that eliminates the need for subcycling and keeping track of the "rigid-motion" history of the stencils for all subcycles.FT-R3M is a nonlinear variant of the rigid motion mesh morpher, since, as it is shown in this article, a nonlinear constraint is imposed.It employs, as it will be explained below, the Gauss-Newton algorithm to solve the nonlinear least squares problem that occur and the Polar decomposition method in order to project the estimation of the rotation matrix that is computed, to the special orthogonal space SO(n) (rotation group), where n is the spatial dimensions of the problem.
The basic idea of FT-R3M is to adapt/deform the internal nodes of the computational mesh, the displacement of which is not prescribed, to a given displacement field of the prescribed nodes (in most cases the prescribed nodes correspond to the boundary nodes of the shape), by keeping parts/elements (a.k.a.stencils) of the mesh as-rigid-as-possible.To do so, the energy functional of the stencils between the initial and the final state of the shape, should be minimal.RIGID MOTION AS A BUILDING BLOCK TOWARDS MESH MORPHING FT-R3M is a mesh-free minimization-based grid deformation technique, in which a target functional is defined, related to the elastic deformation energy, which needs to be minimized in order to propagate the movement of the boundaries to the internal nodes of the CFD mesh.In the subsections that follow, the rigid motion and the total energy functional to be minimized will be introduced.

Rigid motion
A rigid motion (or an isometry) consists of rotations, translations, reflections, and glide reflections (or a combination of them) such that the distance between every pair of points/vectors is preserved.More precisely, a rigid motion is defined as a map  ∶ R n → R n that conserves the following inner product where  ∶ R n → R n is an affine transformation with a rotation matrix.
The idea of the morpher, which is proposed in this article, is to try to mimic such kind of motion.This means that the internal nodes of the computational mesh, which are free to be displaced, are grouped into groups and are forced to follow a motion which is as-close-as-possible to its rigid motion.In the context of the FT-R3M, only in rotations and translations (and the combination of them) are allowed, since reflections and glide reflections produce inverted elements which is undesirable in most software dealing with the solution of partial differential equations.Thus, assuming a rotation matrix R ∈ SO(n) (R ⋅ R T = I n , where I n is the identity matrix) and a translation vector ⃗ t ∈ R n , a rigid motion when acting on any vector ⃗ x, produces the transformed vector (⃗ x) of the form

The energy functional
Let us denote by N the set nodes, N the boundary nodes, Ñ the internal nodes and  the mesh.From now on, let ⃗ X i be the initial position vector of node with index i and ⃗ x i the final one.If the boundary nodes follow a rigid body motion (prescribed motion), we also want to a priori impose to the interior nodes of the mesh the same rigid body motion.The aforementioned can be expressed as To find the solution in such a problem, the pair (R, ⃗ t ) ∈ SO(n) × R n and ⃗ x i ∈ Ñ is needed.This leads to an optimization problem in which the following "energy functional" should be minimized Solving an optimization problem like this in Equation ( 4) is "translated" as follows: Find the closest rigid motion between the initial and the final state.Without loss of generality, Equation (4) can be written on a pair of nodes that form an edge, not necessarily a geometrical edge, instead of the nodes themselves, as follows F I G U R E 1 An example of two neighboring stencils.A simple way to form stencils is by creating a stencil on each node.Then, with the use of a user-defined radius and distance-based criteria, the stencil will consist of the nodes that belong inside a circle on 2D or sphere on three-dimensional meshes (dashed circles, one per stencil).The bigger the radius it is, the bigger the overlapping between the stencils it is, leading to smoother deformations, but increasing the cost of solving the system which is created (to be shown in following chapters) [Colour figure can be viewed at wileyonlinelibrary.com]By minimizing Equation ( 5), the whole mesh is handled as one body, namely, there is no definition of locality.Therefore, in order to make the method more flexible and robust, nodes of the mesh are grouped into cloud of points, which in our context are called "stencils."The role of the stencil is to introduce the locality on the proposed method and capture with more flexibility the deformation of different parts/areas of the computational mesh.Local transformation problem from the initial to the final state are defined and then, all those local problems are glued together by summing them up into a global problem.A stencil is practically a collection of nodes or in other words, cloud of points (graph).There is a freedom in the selection/creation of the stencils.The criteria of forming the stencils are distance-based, meaning that the topology of the mesh plays no role of their formation, thus the proposed method belongs to the family of the mesh-free (mesh-less) methods.An example of two stencils is shown in Figure 1.The overlapping of the stencils is necessary, which means that a node should belong to more than one stencil in order to achieve smooth deformations.The reason is that if the stencils are not overlapping, then the movement of the nodes of two neighboring stencils might be inconsistent, since the displacement of each node is computed by minimizing a global metric consisting from the summation of the local problems.The updated position of an internal node of the computational mesh is a compromise between the rigid motion of each stencil where this node belongs to.By defining the global energy functional as a weighted summation of this of each stencil belonging to the mesh, the final expression of the energy functional to be minimized can be written as In Equation ( 6), S is the stencil set, s is a stencil belonging to this, (i, j) ∈ s denotes a a pair of nodes belonging to stencil s and R s is the rotation matrix of stencil s; w s is a scalar weight per stencil that stresses the importance of some stencils as higher than that of some others' (e.g., in case of a boundary layer), or reflects the volume size of the stencil to mimic a domain integral.Regarding  s,ij , this is a scalar weight per stencil-edge that accounts for mesh anisotropy.It prevents the distortion of a stencil by favoring/penalizing most the edges that are compressed, handling thus mesh anisotropies.Furthermore, by denoting as the total energy functional in Equation ( 6) to be minimized can be rewritten as In the final expression of the energy functional stated in Equation ( 6) and, consequently also in Equation ( 8), there is a "hidden" nonlinear constraint; namely, every matrix R s (one for each stencil of the mesh) should be a "real" rotation matrix that belongs in the special orthogonal SO(n) group.This constraint can be expressed as In addition, if the constraint equation is denoted by then the problem to be minimized (including the constraint) so as to propagate the movement of the boundaries to the interior computational mesh, can be expressed as Equation ( 11) express the energy functional which is required to be minimized in order to to compute the final position of the internal nodes of the computational mesh, when a prescribed displacement field of the boundaries is given.The unknowns in this problem are the coefficients of each matrix R s , ∀s ∈ S and each vector ⃗ x i , ∀i ∈ Ñ.In the chapter that follows, it is shown how the minimization of the nonlinear problem, stated in Equation (11), is treated and solved.

TREATMENT OF THE NONLINEAR PROBLEM
Since the problem to be solved is nonlinear, due to the quadratic constraint of Equation (11), it has to be treated appropriately.In this case, the Gauss-Newton algorithm, which is appropriate for solving nonlinear least squares problems, is employed.According to the Gauss-Newton algorithm, the latter can be derived by linearly approximating the vector ⃗ e s,ij .
Assuming that a current guess (⃗ ) is available and by introducing the increments (⃗ x i , R s ), the energy expressed in Equation (11), using also Equation (7), can be written using the Taylor's theorem as In our case, since ⃗ e s,ij is linear with respect to (w.r.t) its free variables, the linear approximation using Taylor's theorem is exact.The term (⃗ e s,ij ∕R s ∶ R s ) in Equation ( 12) denotes the directional derivative of the vector valued function ⃗ e s,ij w.r.t. the second-order tensor R s in the direction of R s .Considering that at the current step, the matrix R  s is a rotation matrix (to be explained later on how this constraint is ensured) it yields that Taking into account Equation ( 13) and the fact that and some mathematical manipulations, Equation (12) becomes From the linearized set of equations stated in Equation ( 14) it is obvious that the matrix Finally, since , where ⃗ b s is a vector.By taking into consideration the latter and by rearranging some terms in Equation ( 15) we end up with or by denoting as Equation ( 16) is written as Once the displacement field of the prescribed nodes (namely, the displacement of the nodes in the surface mesh) is known, the energy functional in Equation ( 18) is minimized in a least squares sense, by finding the stationary points w.r.t. the corresponding unknowns by satisfying The derivation of the linearized total energy functional w.r.t. the degrees of freedom, namely, ⃗ x i and ⃗ b s , as stated at Equation ( 19), yields a symmetric positive definite system to be solved of the form where N and H are the multipliers of ⃗ x i and ⃗ b s , respectively, which derive from the first term of Equation (19), whereas H T and M are the multipliers of ⃗ x i and ⃗ b s , respectively, deriving from the second term of Equation (19).Thus, if the left-hand matrix is denoted as A, we have where ⃗ u is the vector that consists of the free variables.Using the condensation 31 technique, also known as static condensation, we can eliminate ⃗ b s in terms of ⃗ x i so as to decrease the size of the system that is finally solved.To do so, the second equation of the system (20) needs to be written in terms of ⃗ b s , namely, then, the first equation of the expression (20) becomes After having solved the linear system and (⃗ x i , R s ) are available, (⃗ The solution of the linear system (21) provides the displacement of every node of the mesh (⃗ x i ) and the increment of the rotation matrix (R s ).In this case, the quadratic constraint in Equation ( 9) is not satisfied, which means that R +1 ⋆ s it is not a rotation matrix that belongs to the special orthogonal group SO(n) (more details in chapter 4; that is the reason for denoting the update as R +1 ⋆ s .

PROJECTION ON THE ORTHOGONAL GROUP
Linearizing the quadratic constraint and then solving the linear system in Equation ( 21) does not satisfy the initial constraint, namely that R +1 s should belong to the special orthogonal space SO(n).The linearized constraint, meaning that the matrix Ω +1 s is antisymmetric, ensures that the calculated matrix R +1 ⋆ s lies on the space which is tangent at R  s .Thus, to satisfy the initial constraint in Equation ( 9), the R +1 ⋆ s needs to be projected, for every stencil, on the special orthogonal group.
A way to deal with such a problem, is to solve the so-called Wahba's problem. 32The latter seeks for a rotation matrix between two coordinates systems from a set of weighted vector observations.Assuming that all ⃗ x +1 i are known at the current step, which in our case are the nodal position vectors computed after the minimization of Equation ( 16), the rotation matrix R s of any stencil s ∈ S is computed by minimizing the stencil energy According to Wahba's problem, the minimization of Equation ( 26) is related to the orthogonal Procrustes problem, which is to find the orthogonal matrix R s that is closest to B, where Since R +1 s is a priori imposed to be an orthogonal matrix, with the proviso that the determinant is +1 (otherwise if the determinant is −1, it will still be orthogonal, but a reflection matrix), tr[R +1 s (R  s ) T ] = n, where n is the dimension of the problem and matrix B is known.Hence, the term −2tr(R +1 s B) in Equation ( 27) should be minimized.In conclusion, R +1 s in nothing more than the polar factor of the matrix denoted as B.
In the literature there are many ways to compute the polar factor of a matrix. 33,34A well-known method is the singular value decomposition (SVD), but is not recommended for performance reasons.After performing a wide set of simulations, it has been noticed that 5%-6% of the total runtime is consumed in computing the polar factor of the matrix B. In our context, an iterative method is used based on Heron's method 35 for the square root of 1 (a.k.a.Newton method), which has been proved to be 2%-4% marginally faster CPU-wise than the SVD.These thoughts are also confirmed by Chao et al. 36 The Heron's method does not ensure that the resulting orthogonal matrix will be of positive determinant.This means that, if the determinant is negative, the orthogonal matrix is a reflection matrix which is not desired in our case, since this leads to element inversion.Thus, SVD is used only when the determinant of the orthogonal matrix is negative so as to avoid such a case.

ALGORITHM OF THE MORPHING FRAMEWORK
At this point, it is worth describing the algorithm to be followed in order to obtain the new/updated position of every internal node constituting the computational mesh.Since the problem expressed in Equation ( 11) is nonlinear, an iterative process needs to be followed to achieve convergence.Briefly, the necessary steps to be done in order to compute the final position of the internal nodes of the computational nodes are: 1. Create a graph which contains the stencils and the corresponding nodes belonging to them.A simple way to create the stencils is to define as many number of stencils as the number of the nodes and use distance-based criteria to find the neighbors.The number of neighbors that belong to the stencil is defined by a user defined radius.The bigger the radius it is, the more neighbors in the stencil are included.This step needs to be done only once, before starting the iterative process of computing the displacement of the nodes of the computational mesh.

Initialization of ⃗
x , where I n is the identity matrix.3. Computation of the derivatives (Equation 19) required to build the system described in Equation ( 20).Solving the above-mentioned linear system (or the one that has been condensed, Equation ( 23)) with the use of least-squares method and computation of (⃗ , namely, the nodal positions and an estimation of the rotation matrix, respectively.This step is the so-called prediction step.4. Since this solution is obtained by linearizing the constraint in Equation ( 9), R +1 ⋆ s is not a rotation matrix that belongs to the orthogonal group SO(n), but it belongs in the space which is tangent to this.Thus, it has to be projected back to SO(n).In this step, the position vectors ⃗ x +1 i are kept fixed and, the R * s for every stencil are projected to the orthogonal group SO(n) to obtain the closest rotation matrix R +1 s .This step is called correction step.

Check if steady state has been achieved, namely, if ⃗
x i , R s have converged, namely, if the mesh is not changing anymore or in other words when ||⃗ x n − ⃗ x n−1 || 2 < , where  is a very small threshold.Since the rotation matrix is analytically computed, it will also be converged if the vector ⃗ x i is also converged.If not, the iterative process continues returning to step 2 until convergence is achieved.

APPLICATIONS
In this section, the FT-R3M proposed in this article, is being demonstrated and tested as a standalone tool to show its efficiency in some test cases.Metrics are used to quantify the grid quality of the initial and the final computational mesh.More specific, the value of the skewness, the nonorthogonality 37 (those used by OpenFOAM ® ) and the aspect ratio are used.
For the sake of completeness, it is worth mentioning that the nonorthogonality metric (to be denoted as N-O) measures the angle between the line connecting two cell centers and the normal of their common face (the lower the better).The skewness (to be denoted as SK) metric measures the distance between the intersection of the line connecting two cell centers with their common face and the center of that face (smaller is better).The nonorthogonality and skewness metrics are analyzed more detailed in the thesis of Jasak 38 in pages 83 and 124, respectively.Last, for the aspect ratio (to be denoted as AR), the longest edge of an element is divided by the shortest one.
OpenFOAM ® is able to provide reliable simulations on meshes that the nonorthogonality of the cells of the computational mesh is smaller than 70.If the nonorthogonality is between 70 and 90, then special treatment is required, for example, nonorthogonal correctors.Regarding the skewness metric of the cells of the computational mesh, OpenFOAM ® provides solution with good/acceptable accuracy if the maximum value of the skewness metric in the mesh is smaller than 0.95.If the maximum value of the skewness in the computational mesh is bigger than the latter value, then difficulties in convergence may occur.

Rotation of an isolated airfoil
The first case deals with the rotation of a two-dimensional (2D) isolated airfoil anticlockwise around its center.In this case, the far-field boundary of the airfoil is kept fixed, namely, the displacement field on that is 0. The computational mesh is appropriate for viscous flow and consists approximately of 32K nodes and 60K elements, with a boundary layer near the walls (quadrilateral elements) and triangular elements everywhere else.
The initial mesh and shape of the airfoil is demonstrated in Figure 2. The rotation of the airfoil takes place every 5 • so as to keep track of the quality metrics for each position of the airfoil.However, the morphing procedure is performed between the initial shape and the desired position of the airfoil, since no subdivision of the displacement field is needed when the FT-R3M is used.
The airfoil is demonstrated in Figure 5 in various positions.A focus has been made in the leading and trailing edge on the resulting mesh in Figure 6, after 90 • of rotation, for demonstration purposes and to show that the orthogonality of the boundary layer is preserved.
The history of the aforementioned grid quality metrics is shown in Figure 3.It is shown that FT-R3M preserves the grid quality after the deformation, even in extreme cases like that.In the latter, each one from the aforementioned metrics is normalized with corresponding one in the initial state of the airfoil.From Figure 3, we can observe that, the aspect ratio (mean value, maximum value, and SD) and the skewness are not changing significantly during the deformation.On the other side, the maximum and the mean value of the nonorthogonality of the elements is increasing (becoming worse) as the degrees of rotation are growing.This is reasonable since the deformation of 90 • is extreme.However, even in such deformation there is no occurrence of any inverted element.Indicatively, the cost of the morpher on this specific case it is compared with the cost that is needed to run the equivalent flow equations (EFS).If M is the cost of solving the EFS, then the cost of the morpher (e.g., for 10 • of rotation) is 0.05M.At this point, it worth mentioning that the cost of the EFS w.r.t. the cost of morpher is big since the flow is turbulent.Moreover, the cost of the morpher also depends on the nature of the displacement field.In other words, the bigger the rotation of the elements caused by the movement of the boundaries it is, the heavier (computationally) the system to be solved becomes.
Regarding the quality metrics and more specific the aspect ratio (normalized value) of the cells belonging to the boundary layer, these are tabulated on Figure 4.In this figure it is shown that the aspect ratio of the cells of the boundary layer is preserved during the deformation.The mean value, the SD ant the maximum value of the aspect ratio of the cells of the boundary layer is changed approximately 0.0012% during the deformation of the airfoil in comparison to the initial one, whereas the minimum value of the aspect ratio has changed approximately 0.010% in comparison to the initial one (Figures 5 and 6).

Deformation of a duct
The second case is dealing with the deformation of an S-Bend duct (Figure 7) which is an air duct test case provided by VolksWagen AG.The displacement field has been provided by the NTUA group after performing an adjoint-based shape optimization for minimum total pressure losses across the duct. 39The flow is laminar, with inlet velocity U in = 0.1m/s, normal to the boundary, leading to a Reynolds number of approximately Re = 400, based on the inlet hydraulic diameter.The mesh is comprised of approximately 470K hexahedrals and 480K nodes.The final shape of the duct is the optimized one with approximately 60% less total pressure losses than the initial one.The challenge on this case is to deform the mesh from the initial to the optimized state and maintain a good quality.The final shape of the duct is shown in Figure 8 whereas in Figure 9, cut of the final shape is presented in order to demonstrate the displacement of a few of the internal nodes of the computational mesh.
In addition, quality metrics for the initial and the final mesh are tabulated in Table 1 to demonstrate the efficiency of FT-R3M.The CPU cost of the deformation of the duct is approximately 0.13M.Furthermore, in Figure 10, the streamlines of the initial and the deformed/optimized duct are presented.In the initial geometry, there is a high recirculation area that takes place in the part of the duct which is near to the outlet.This has the biggest contribution to the total pressure losses, which is not the case in the optimized geometry, where the streamlines are smoother, hence the reduced total pressure losses.

Deformation of a beam
The next case, in which FT-R3M is used as the mesh morphing tool, deals with the deformation of a beam with a rectangular cross-section (Figure 11).The beam fixed at a wall from the one side and the purpose of this application is to apply extreme deformations and adapt the mesh to the updated boundaries (new shape of the beam) and check its grid quality.The unstructured mesh has been generated using CFD-GEOM 40 (ESI's grid generator) and consists, approximately, of 195K elements (tetrahedrals) and 79K nodes.The deformations applied in this test case are extreme and the reason is to show the efficacy of FT-R3M.Thus, a twisting, a bending deformation and a combination of them is applied to the right surface of the beam.Quality metrics that used in previous chapters are also used here and are quantifying the quality of the resulting grid.The aforementioned metrics are normalized with the corresponding one of the initial mesh and tabulated in Table 2.In Figures 12,13

CONCLUSIONS AND DISCUSSION
The mesh morphing/movement tool for adapting the computational mesh to the updated boundaries was developed in this article and was shown to be capable of handling big deformation fields very efficiently.It was shown that, the computation of the final position of each node of the computational mesh requires the solution of a linear system, Equation ( 21)  and, a projection of the estimation of the rotation matrix of each stencil into the special orthogonal group SO(n) till convergence is achieved.The energy functional, related to the deformation energy, between the initial state and the final one is minimized so that, at the same time a good quality of the final mesh is ensured, according to several used quality metrics.Moreover, since FT-R3M requires the formation of the so-called stencils, which are practically cloud of points and do not depended on the grid connectivity, it makes the grid displacement tool a mesh-free/mesh-less method and can naturally handle any kind of mesh, consisting even of nonconvex cells and nonplanar faces.
The cases shown in this article demonstrate the efficiency of the FT-R3M in big/extreme deformation fields.Mesh anisotropies are intrinsically handled, by favoring rigidity in critical directions and, the orthogonality of the boundary layers are very well preserved.Such a tool is essential, for example, in aerodynamic shape optimization since during the optimization, the computational mesh needs to be adapted to the updated boundaries and most important to be of good quality.For instance, maintaining the mesh quality near the body, for example, the boundary layers, is very important in high-Reynolds number or in heat transfer problems in order to have accurate results of the analysis tool.
Coupling FT-R3M with a gradient-based optimization algorithm, for example, with the adjoint method, [41][42][43] might show the efficacy of the proposed mesh movement tool inside the optimization loop.For instance, since it preserves the quality of the mesh, the optimization can go further without the need of remeshing, since it is very common that the optimization stops because of low quality of the resulting mesh.

F I G U R E 2 3
Airfoil case: Initial mesh and shape of the isolated airfoil.A focus has been made on the right figure near the boundaries of the airfoil [Colour figure can be viewed at wileyonlinelibrary.com]Airfoil case: History of the quality metrics during the deformation of the airfoil.Each of the metrics are normalized with the corresponding metric of the initial mesh.AR, aspect ratio; N-O, nonorthogonality metric; SK, skewness [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 4
Airfoil case: History of the aspect ratio of the cells belonging to the boundary layer during the deformation of the airfoil.Each of the quantities are normalized with the corresponding one of the initial mesh [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 5 F I G U R E 6
Airfoil case: The resulting mesh in 25 • , 50 • , 75 • , and 90 • anticlockwise [Colour figure can be viewed at wileyonlinelibrary.com]Airfoil case: Focus on the leading (left) and trailing (right) edge of the rotated airfoil (90 • ).The orthogonality of the boundary layer in the rotated airfoil has been preserved which is very crucial for viscous flows [Colour figure can be viewed at wileyonlinelibrary.com] , and F I G U R E 7 Sbend case: Initial mesh and shape of the duct [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 8 Sbend case: Final/Deformed shape and mesh of the duct.Big displacements are presented in the central curved part of the duct [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 9 Sbend case: Cut of the final shape of the duct [Colour figure can be viewed at wileyonlinelibrary.com] 14, the beam is presented in various positions.The presence of the wall is just for demonstration reasons and it does not affect the morphing process.It indicates only which side of the beam is fixed at the wall.

F I G U R E 10
Sbend case: Velocity streamlines in the initial (up) and the deformed (bottom) shape of the duct [Colour figure can be viewed at wileyonlinelibrary.com]TA B L E 1 Sbend case: Quality metrics before and after the deformation of the duct Normalized quality metric State avg AR std AR max AR min AR avg N-O std N-O max N-O avg SK std SK max SK Before 1.000

F I G U R E 11
Beam case: Initial shape and mesh of the beam.The wall on the left of the figure appears only for demonstration purposes and has no effect on the morphing process [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 12 Beam case: The beam has been twisted 350 • over the axis which is parallel to its length [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 13 Beam case: Bending deformation of the beam [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 14
Beam case: Combination of twisting and bending deformation of the beam [Colour figure can be viewed at wileyonlinelibrary.com]

AR max AR min AR avg N-O std N-O max N-O avg SK std SK max SK
Each quality metric is normalized with the corresponding one of the initial mesh.Abbreviations: AR, aspect ratio; N-O, nonorthogonality metric; SK, skewness.Quality metrics for each case of the beam Each quality metric is normalized with the corresponding one of the initial mesh.Abbreviations: AR, aspect ratio; N-O, nonorthogonality metric; SK, skewness.
Note:TA B L E 2Note: