Finite Transformation Rigid Motion Mesh Morpher

In any optimization framework, a robust and reliable mesh morpher is necessary to undertake the adaptation of the CFD mesh to the updated boundaries at each optimization cycle. Morphing has its share of challenges, namely to maintain high mesh quality (avoid distorted elements and tangles) even during extreme deformations. In this work, the Finite Transformation Rigid Motion Mesh Morpher (FT–R3M) is presented, an improved version of the Rigid Motion Mesh Morpher (Eleftheriou and Pierrot in Rigid motion mesh morpher: a novel approach for mesh deformation, 2016), that eliminates the need for sub-cycling, making it more efficient in terms of CPU time. FT–R3M, which bears some similarities to Chal et al. (ACM Trans Graph 29(4):38, 2010), is a mesh–less mesh morphing tool, since it does not require any inertial quantities, that 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, namely the deformation energy. Since FT–R3M is minimizing the deformation energy between the initial and the final configuration, as opposed to R3M, in which the deformation energy is minimized from each sub-cycle to another, there is a significant gain in terms of the quality of the resulting mesh. The efficiency of the morpher proposed in this article will be demonstrated in small and medium–size cases.


Introduction
A variety of stochastic and gradient-based methods have been devised to solve aerodynamic shape optimization problems such as the design of airfoils or wings for optimal drag and/or lift, ducts for minimum power losses, cars with optimal combination of drag and lift etc. In these kind of optimization problems, it is required to use a technique so as to deal with the necessary changes of the boundaries (surface mesh), namely to adapt the computational mesh to the updated geometry in order to proceed with the optimization process.
One well-known method to handle the changes in the boundaries during optimization, is remeshing, in which after every optimization cycle a new grid is generated. This is a time-consuming process and especially for gradient-based methods, in which the information needed for the shape changes is retrieved from the gradient with respect to free design variables, the gradient consistency is lost. Moreover, for complex geometries, manual intervention in the mesh generation will be needed, which makes this method neither robust nor reliable. On the other hand, a promising way to propagate the movement of a shape to the interior mesh is to make use of a mesh deformation tool (mesh morpher). This requires to generate a mesh only once, at the beginning of the optimization process and then, the mesh morpher will undertake to deform this mesh at each optimization cycle. In this article, the Finite Transformation Rigid Motion Mesh Morpher will be presented and demonstrated.
The challenge when a mesh morpher is used, is to maintain a good mesh quality of the deformed mesh, namely to avoid negative or distorted elements, highly skewed cells etc, which may cause divergence of the CFD software, in CFD-based problems. In literature, there is a variety of mesh deformation techniques. Most of them require a trade-of between attained mesh quality and CPU cost. A simple and easy way, in terms of implementation complexity, to deform a mesh is by using the Laplacian smoothing [7,15], or the so-called Laplacian coordinates. The computation of the displacement related to the nodes of the CFD mesh requires the solution of a linear system whose dimension scales with the number of nodes of the mesh. This technique proves to be efficient but not very robust, since it cannot handle mesh rotations and generally complex transformations. A rival mesh deformation tool is the Linear Elasticity morpher, in which the computational mesh is handled as an elastic solid body [10,14]. In this method, mesh deformation is accomplished by solving the linear elasticity equations for the mesh point inside the geometric domain. Since the elasticity equations contain material properties (Young's modulus and Poisson's ratio), these are related to the mesh characteristics. Another popular mesh deformation method is the Spring Analogy Method [1]. In 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 [6] or semi-torsional springs [2], but still mesh anisotropies cannot be handled. A very promising mesh deformation technique is the mesh morpher based on Radial Basis Functions (RBF) [9]. Mesh morphing based on RBF is a robust and reliable mesh deformation technique but computationally very "heavy", as the matrices involved in the computations are very dense.
The Rigid Motion Mesh Morpher has already been introduced [5]. It 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, that eliminates the need for sub-cycling and keeping track of the 'rigid-motion' history of the stencils for all sub-cycles. It employs, as it will be explained below, the Polar decomposition method in order to project an estimation of the rotation matrix 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 nodes, whose displacement 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 of the mesh asrigid-as-possible, hence keeping the deformation energy of these parts/elements, between the initial and the final state of the shape, minimal. There is a significant gain in terms of CPU and morphing efficiency.

Rigid Motion as a Building Block Towards Mesh Morphing
The Finite Transformation Rigid Motion Mesh Morpher is a minimization-based approach. There is a target functional, namely the deformation energy in our context, which has 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 deformation energy to be minimized will be introduced.

Rigid Motion
A rigid motion (or an isometry) consists of rotations, translations (or a combination of them) such that the distances between every pair of points/vectors is preserved. In particular, a rigid motion is defined as a map φ : R n → R n that conserves the inner product ∀ x, y ∈ R n , φ (x), φ (y) = x, y 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 Deformation Energy
Let us denote by N the set nodes, ∂ N the boundary nodes and N the internal nodes. 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, we also want to apriori impose to the interior nodes of the mesh the same rigid body motion To this end, the pair (R, t)∈SO(n)×R n and x i ∈ N is needed, that minimizes the 'energy functional' (deformation energy) By minimizing expression 4, the whole mesh is handled as one body. In order to make the method more flexible and robust, we group some parts/elements of the mesh into 'stencils'. A stencil can be a collection of edges, not necessarily geometrical edges. In our context, an edge is a pair of nodes sharing a common cell. There is a freedom in the selection of the stencils; in the simplest case, a 1:1 ratio between the number of stencils and the number of elements is considered. Decomposing the final deformation energy as a weighted summation of the deformation energy of each stencil belonging to the mesh, the final expression of the deformation energy to be minimized can be written as In Eq. 5, S is the stencil set, s is a stencil belonging to this and (i, j) ∈ s denotes an edge belonging to stencil s; w s is a scalar weight per stencil that stresses the importance of some stencils as higher that some others' (e.g. in case of a boundary layer) and, µ s,i j a scalar weight per stencil-edge that accounts for mesh anisotropy, by preventing the distortion of a stencil by favouring rigidity in directions in which distortion is imminent. Over and above, by denoting as the total deformation energy in eq. 5 to be minimized can be rewritten as In the final expression of the deformation energy stated in eq. 5, there is a hidden non-linear constraint; namely, every matrix R s (one for each stencil of the mesh) should be a 'real' rotation matrix inside the orthogonal SO(n) group, which can be expressed as ∀s ∈ S, R T s R s = I n Finally, the Energy to be minimized can be expressed as

Algorithm of the Morphing Framework
After having introduced the total deformation energy in eq. 9, it is worth describing the algorithm to be followed in order to minimize it and, obtain the new/updated position of every internal node constituting the CFD mesh. Since the problem expressed in eq. 9 is non-linear, it has to be linearized and follow an iterative process to achieve convergence. Initially, the quadratic constraint stated in eq. 9 should be linearized. After that follows the derivation of the linearized deformation energy w.r.t. the unknowns, namely the nodal positions and a skew-symmetric matrix for each stencil (to be explained later on) and the assembly of the linear system to be solved. The necessary steps to be undertaken by the morpher are: 1. Solution of the linear system and computation of (x κ+1 i , R * s ) from (x κ i , R κ s ), the nodal positions and an estimation of the rotation matrix respectively. Since this solution was obtained by linearizing the constraint in eq. 9, R * s will not be a rotation matrix that belongs to the orthogonal group SO(n), thus it has to be projected back to SO(n) in order to find the closest rotation matrix to R * s . This step is called prediction step. 2. Computation of the closest rotation matrix to R * s . 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. 3. Check if steady state has been achieved, namely if x i , R s have converged. If not, the iterative process continues returning to step 1.

Building the System of Equations
Before going to the linearization of the total deformation energy, it is worth mentioning a few things about the structure of the orthogonal group and its tangent space, that will also be used during the linearization of the energy. Furthermore, in this section, the linearized energy will be constructed, as well as the necessary steps to build the final system to be solved whose solution is the position of each node of the mesh.

Structure of the Orthogonal Group and its Tangent Space
To better understand the implicit quadratic constraint in eq. 9, the structure of the orthogonal group and its tangent space will be briefly explained. Let us denote where GL(n, R) is a general linear group of degree n (the multiplication of two invertible matrices is also an invertible matrix). In addition, let us denote as Since ϕ(R) = 0 n is singleton (unit set; contains exactly one element), it is also a closed set. Moreover, ϕ(R) is continuous, hence SO(n) is a closed set. Let R ∈ SO(n) and B ∈ GL(n, R) be given and ε ∈ R, then we have Hence, the Gâteaux derivative of ϕ is Since dϕ(R, B) is linear and continuous, it follows that SO(n) is a differentiable manifold and that its tangent space at R is the null space of dϕ(R, Ω ): which means that the space tangent to the orthogonal group SO(n) at R is the space of linear transformations B such that R T B is skew-symmetric.

Linearization of the Problem
Assuming that a current guess (x κ i , R κ s ) is available and by introducing the increments (δ x i , δ R s ), the energy expressed in eq. 9 (using also eq. 6) can be rephrased as Now, the linearization of the constraints takes place at this point, by neglecting the second order terms (assuming that δ R s δ R T s = 0). Then, introducing Ω s = δ R s R T s , a skew-symmetric matrix (see also section 4.1), eq. 15 becomes Finally, since Ω s is antisymmetric, it implies that the product Ω s (X j − X i ) can be replaced by b s × (X j − X i ) with b s a vector. By taking into consideration the latter and by rearranging some terms in eq. 16 we end up with Once the displacement of the prescribed nodes (namely the displacement of the nodes in the surface mesh) is known, the deformation energy in eq. 17 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 deformation energy w.r.t. the degrees of freedom, namely δ x and b s , as stated at eq. 18, yields a symmetric positive definite system to be solved of the form where u is the vector that consists of the free variables. Using the condensation [13] technique, also known as static condensation, we can eliminate one degree of freedom in terms of the other so as to decrease the size of the system.

Projection on the Orthogonal Group
The linearization of the constraint in eq. 15 leads to obtain a matrix R * s which lies on the tangent space of the orthogonal group SO(n) after minimizing the eq. 17, which practically means that it is not a 'real' rotation matrix, thus R * s R * T s = I n . The latter makes the projection of the estimated R * s in the orthogonal group, essential. A way to deal with such a problem, is to solve the Wahba's problem [16]. The latter seeks for a rotation matrix between two coordinates systems from a set of weighted vector observations. Assuming that all x i are known, which in our case are the nodal position vectors computed after the minimization of Eq. 17, the rotation matrix R s of stencil s ∈ S is computed by minimizing the stencil energy According to Wahba's problem, the minimization of eq. 20 is related to the orthogonal Procrustes problem, which is to find the orthogonal matrix R s that is closest to B, where B = ∑ k µ s,i j (x j − x i )(X j − X i ) T , in the sense of Frobenius norm Since R 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) = n, where n is the dimension of the problem and matrix B is known, the term −2tr(R s B) in eq. 21 should be minimized. Hence R 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 [11,8]. Herein, an iterative method is used based on Heron's method for the square root of 1, which has been proved to be very efficient.

Applications
In this section, the Finite Transformation Rigid Motion Mesh Morpher, proposed in this work, is being demonstrated and tested as a stand-alone tool to show its efficiency in three cases, in particular in the moving boxes, in a rotating airfoil and in a beam fixed at the wall.

Moving Boxes Case
This case deals with the movement of a box, which is placed inside a larger box. It is a simple case, at first sight, though it can provide useful information for the efficiency and the behaviour of the morpher proposed in this work.
The mesh in between the two boxes consists of 2400 hexahedral elements and approximately 5000 nodes. In the first case, the inner box is rotated 80 degrees anticlockwise around its center and in the second one, it is translated in x-direction (by one time its box width). In both cases, the outer box is kept fixed (undeformable) and the nodes in between the two boxes are adapted to each movement. The initial mesh is presented in Fig. 1 whereas the deformed meshes for both cases, are presented in Fig. 2. The purpose in the first case is to test the FT-R3M in mesh rotation whereas in the second one, is to show how the proposed mesh morpher treats the squeezed elements. In order to quantify the mesh quality, two quality metrics are used. Table 1 tabulates these quality metrics (maximum value of skewness and non-orthogonality metric) for each mesh (initial and deformed). For the sake of completeness, it is worth mentioning that the non-orthogonality metric measures the angle between the line connecting two cell centers and the normal of their common face (the lower the better) whereas the skewness 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). In this simple case, it is illustrated that the Finite Transformation Rigid Motion Mesh Morpher is able to handle mesh rotations and, at the same time, maintain a good mesh quality in such extreme deformations. FT-R3M can also prevent the distortion of an element (squeezed elements) by keeping these stencils more rigid than others for which the distortion is not imminent.

Deformation of a 2D airfoil
The second problem is dealing with the rotation of a 2D airfoil. The computational mesh is appropriate for viscous flows and consists of approximately 32K nodes and 60K elements, with a boundary layer at walls (quadrilateral elements) and triangular elements everywhere else. The airfoil is rotated anti-clockwise and FT-R3M undertakes the CFD mesh adaptation to the rotated boundaries. Moreover, quality metrics (in particular skewness and non-orthogonality) are used to quantify the quality of the resulting mesh. The purpose is to demonstrate the utility of the proposed mesh morpher in extreme deformations in an anisotropic mesh. In Fig. 3, the deformation energy on the initial and the deformed mesh of the airfoil are presented. In Fig. 4 a focus has been made on the leading (left) and the trailing (right) edge of the airfoil on the resulting mesh, after 90 o degrees of rotation for demonstration purposes. The skewness and the non-orthogonality metrics for the initial and the resulting mesh are tabulated in table 2. The rotation of 90 o degrees of the airfoil is just before an inverted element appears. In Fig. 4 it can be observed that FT-R3M preserves the orthogonality of the boundary layer on the deformed airfoil.  Table 2 Mesh quality metrics before and after the deformation of the airfoil.

Deformation of a Beam
The last case, in which FT-R3M is tested, deals with the deformation of a beam with a rectangular cross section which is fixed at a wall. The aim is to apply extreme deformations in the free cross-section and adapt the mesh in the new shape of the beam. The initial shape and mesh of the beam is demonstrated in Fig. 5; the mesh has been generated using CFD-GEOM and consists of approximately 195K elements (tetrahedra) and 79K nodes. The deformed beam is demonstrated in Fig. 6 and 7 for different displacement fields (bending and twisting deformation respectively). In addition, the same quality metrics that have been used in previous chapters are also used here to monitor the quality of the deformed mesh and are tabulated in Table 3. Because of the extreme deformation field that is applied, the non-linear iterations needed in this case are more in number than these on the airfoil and the moving boxes cases.    Table 3 Quality metrics for the bending deformation (left) and the twisting deformation (right).
In this particular case, the quality of the mesh is kept almost the same during the deformation (Table. 3). FT-R3M is capable of handling extreme deformations and at the same time ensuring a good mesh quality.

Conclusion
In this work, a mesh deformation tool, the Finite Transformation Rigid Motion Mesh Morpher, has been developed and introduced. FT-R3M belongs to the family of optimization-based methods, since the internal nodes of the mesh are displaced to minimize a distortion metric. It has been demonstrated that FT-R3M can handle intrinsically mesh anisotropies and mesh rotations even in extreme deformations and maintains a good quality of the resulting mesh.
In this article, FT-R3M has been demonstrated as a stand-alone tool, thus as ongoing work, it will be implemented as a part of a fully automated discrete adjoint solver loop [12]. Future work may include smoothing of the surfaces, which is necessary in an adjoint solver loop, since the adjoint sensitivity vector usually contains numerical noise. Moreover, the functionality of the FT-R3M may be extended in order to be used as static deformation tool, namely as a mesh optimization tool.