Rigid Motion Mesh Morpher: A Novel Approach For Mesh Deformation

OPT-i An International Conference on Engineering and Applied Sciences Optimization M. Papadrakakis, M.G. Karlaftis, N.D. Lagaros (eds.) Kos Island, Greece, 4-6, June 2014 RIGID MOTION MESH MORPHER: A NOVEL APPROACH FOR MESH DEFORMATION George S. Eleftheriou1,2 , Guillaume Pierrot1 1 ESI-Group 99 rue des Solets, Parc d’Affaires SILIC, 94513 Rungis CEDEX, France 2 National Technical University of Athens P.O. Box 64069, 15710 Athens, Greece e-mail: {george.eleftheriou,guillaume.pierrot}@esi-group.com Keywords: rigid motion mesh morpher, as-rigid-as-possible, mesh morphing, mesh deformation, mesh anisotropy, mesh rotation, distortion, mesh-less, optimization-based, cfd, computational fluid dynamics, adjoint based optimization Abstract. In adjoint based shape optimization problems, after the sensitivities have been computed, there are two ways of dealing with the necessary changes to the mesh. The first one is re-meshing based on the new shape, while the second one is adapting the existing mesh (by moving the nodes) to fit the new shape. Re-meshing may be time consuming, as it is introduced as a separate step inside the optimization loop, and tedious as it may require by-hand manipulation. It also introduces inconsistencies in the process as the sensitivities have been computed at isoconnectivity. Morphing also has its share of challenges, namely maintaining the mesh quality (avoiding distorted and negative cells) while deforming it. Within this context, various mesh morphing techniques have been developed. The Spring Analogy [1] is simple but may suffer robustness issues. Laplacian smoothing [11] is suitable for translation but does not account for rotation. The linear elasticity approach [3] does not account for mesh anisotropy and is difficult to implement for general meshes because finite elements are used to solve the equations. Finally, the Radial Basis Functions [7] are promising but computationally heavy as the matrices involved are full, restricting the mesh size and complicating the implementation. The Rigid Motion Mesh Morpher approach, proposed in the present study, aims at overcoming the above-mentioned limitations, being more flexible and essentially mesh-less, since it does not require any inertial quantities or cell connectivities related to the mesh. Firstly, the set of boundary nodes (nodes defining the shape) of the mesh is identified. The prescribed motion of these nodes (their velocities) is known. Then, all nodes are grouped into “stencils” which are required to deform in an as-rigid-as-possible way. Hence, every stencil has a rotation velocity and a translation velocity. Those, along with the velocities of all internal nodes, form the set of unknowns. The as-close-to-rigid-as-possible motion is ensured by attempting to minimize a metric representing the difference of the actual deformation from a perfectly rigid motion (a translation plus a rotation). It will be shown that this quantity is related to the anisotropic deformation energy. Next, the resulting system of equations is solved, having the prescribed motion of the boundary nodes as boundary conditions. George S. Eleftheriou, Guillaume Pierrot 1 INTRODUCTION Let us consider the general shape optimization problem which, in gradient-based optimization problems, can be described as ~ min J(U, d) ~ d∈D ~ =0 R(U, d) (1) (2) With • J: the objective function • U : the system state given by the natural problem R (Navier-Stokes equations) defined on the domain Ω (if steady) or Ω × (0, T ] (if unsteady) • d~ ∈ D ⊂ RD , D ≥ 1 the design variables of the model. In our study this vector is the nodes coordinates of the mesh Mh related to the physical domain Ω. Let us denote by {Ni }N i=1 the set of these nodes and by ΩB and ΩI respectively the boundary and the interior of the domain Ω. II = {i/ Ni ∈ ΩI } and IB = {i/ Ni ∈ ΩB } are the sets of node indices belonging respectively to ΩI and ΩB . In the case of a gradient-based method, the −∇d~J at each optimization iteration gives the best descent direction. The computed gradient concerns the boundary nodes and the interior nodes must somehow follow their movement. This can be achieved by trashing the existing mesh and creating a new one (re-meshing) based on the new shape. However it is timeconsuming and introduces inconsistencies, as the connectivities of the sequence of meshes produced differ from iteration to iteration. It also requires by-hand manipulation. Another approach is to just make the internal nodes follow the movement of the boundary nodes defining the shape, thereby deforming the existing mesh. The cell connectivities are kept intact during the process and there is no new mesh to be built. This is called mesh morphing or mesh smoothing. The problem of this approach is that the movement does not necessarily end up well for the resulting mesh quality. It might be severly distorted, having negative cells. Within this context, several mesh morphing methods have been developed. They are classified in four principal categories described in the subsections 1.1, 1.2, 1.3 and 1.4. We attempt to quickly cite the principles of the techniques along with their weaknesses and merits. To facilitate notation we can globally define as ~g = (~ui )i∈IB the displacement prescribed over ΩB . 1.1 Spring analogy The spring analogy method [4, 2, 12, 1] is the most simple and classical approach, performing quite well in a number of cases. Unfortunately, it doesn’t prevent cell collapse, and fails when the local mesh motion exceeds significantly the local mesh size. Things become even worse, when it has to deal with mesh anisotropy and mesh rotation. Its principal idea is the introduction of a spring energy Es and the requirement for it to be minimum Es = where Lij is the length of the spring. X 1 (uj − ui )2 L ij ij (3) George S. Eleftheriou, Guillaume Pierrot 1.2 Laplacian smoothing Also known as “Laplacian coordinates” or “Laplacian based”, this approach is based on the solution of a simple Laplace equation [10, 13, 11, 6] such that the displacement ~u of the nodes’ coordinates is the solution of the following problem:  ∇ · (µ(~x)∇~u) = 0 in Ω (4) ~u = ~g on ΩB with µ(~x) as the diffusion coefficient selected over Ω preserving the mesh anisotropy. The computation of the displacement ~u related to the internal nodes (4) comes down to solving a linear system whose dimension is equal to the number of nodes. This method seems efficient when the mesh is subject to translation but suffers robustness issues in cases of more complex transformations such as rotation. 1.3 Linear elasticity In this method [3, 10, 5], in order to handle rotation, a linear elasticity model is introduced to compute the displacements ~u such that:  ∇ · σ(~u) = 0 in Ω (5) ~u = ~g on ΩB σ being the Cauchy stress tensor  σ = λ tr(ε(~u)I + 2µε(~u) ε(~u) = 12 (~uT + ~u) (6) with λ and µ being the Lamé constants. This method reaches its limits with large rotations or mesh anisotropy as it doesn’t handle them intrinsically. Furthermore it is difficult to implement for general polyhedral meshes because finite elements are used to solve the equations. 1.4 Radial Basis Function interpolation This is an interpolation method [8, 7] searching for the displacement ~u(~x) X ~u(~x) = γl φ(k~x − ~xl k) + h(~x) (7) l ~u is the smooth interpolant which minimizes a certain norm. φ is the radial basis function [7]. The constraint is to be exact over the boundary nodes: ~u(~ xi ) = ~v(~ xi ) i ∈ IB . (8) B The search for {γl }N l=1 satisfying (8), requires the solution of a linear system, the matrix of which is unfortunately dense. This causes difficulties in cases of large NB . (NB +1)×(NB +1) B In order to determine the coefficients {γl }#I l=1 we have to solve a linear system in R where NB = #IB . Compared to the previously cited techniques, RBF morphing proves quite robust with respect to both mesh rotation and mesh anisotropy. The problem is that the matrices involved in the computations are always dense by design. Therefore, it’s a matter of choice between a simple implementation at high computational cost, or a complex implementation at the cost of solving the system almost as if it was sparse. George S. Eleftheriou, Guillaume Pierrot 2 THE RIGID MOTION MESH MORPHER CONTRIBUTION The method belongs to the family of optimization-based methods according to the classification cited in [9]. This means that the nodes are not moved based on a heuristic algorithm such as Laplacian smoothing, they are moved to minimize a given distortion metric. More specifically, it tries to favour rigidity in the critical directions of imminent distortion. There is an inherent contradiction starting from the very name of this method. Rigid Motion Mesh Morpher. Why Rigid Motion? Afterall, how could something (a mesh in particular) be kept rigid while being deformed (morphed)? The stated objective of this method and its main philosophy is to deform a mesh, meanwhile keeping elements/parts of it as-rigid-as-possible. Obviously they cannnot be kept entirely rigid, as this would defeat the purpose, however there are ways to favour rigidity in certain critical directions when distortion of the element/part becomes imminent. These elements/parts are essentially groups of nodes called stencils. During the mesh deformation, there is a kind of prevention mechanism that kicks-in whenever a particular stencil is about to distort or degenerate to a line segment. The stencils can be defined with freedom as to the choice of nodes belonging to them. In the simple case, a 1:1 ratio between the number of stencils and nodes could be considered as follows; a stencil may consist of one node, plus all of its adjacent/neighbouring nodes sharing a cell. The main qualities of the method is that it is mesh-less and handles intrinsically mesh anisotropy and mesh rotation. Let us assume that an arbitrary rigid stencil s is in motion. Its motion comes as the combination of a translation, under velocity a~s and a rotation under velocity b~s . In this case, the velocity of a node j belonging to it, denoted as υ~sj would be υ~sj = a~s + b~s ∧ (x~j − c~s ) (9) where x~j is the coordinate vector of the node and c~s is the coordinate vector of the stencil’s center of gravity. In the simple case, its coordinates can be computed as the mean value of the coordinates of all nodes constituting the stencil (center of gravity of a point cloud). Moving further, we no longer require that the stencil s be rigid; it is now allowed to deform. Therefore the velocity of the node, υ~j would no longer be equal to υ~sj υ~j 6= a~s + b~s ∧ (x~j − c~s ) = υ~sj (10) Let e~sj be a vector representing this difference, between the actual deformation of the stencil and its perfectly rigid motion. This vector will be a function of the stencil s and the nodes j belonging to it: e~sj = υ~j − υ~sj = υ~j − a~s − b~s ∧ (x~j − c~s ) (11) Taking control over how this vector evolves/develops for every single stencil of the mesh and trying to bound/minimize it, to the extent that this is possible, gives the first idea of how the optimization problem will be posed. Therefore we enhance/enrich it with some weighting coefficients e~sj = √ ws µsj (υ~j − a~s − b~s ∧ (x~j − c~s )) (12) The square root is there just for convenience and the need for it will become clearer in section 5.3; ws is a scalar weight per stencil that may arbitrarily stress the importance of some stencils as higher than others’. It can be objective when set to 1.0. To better explain the need for a per stencil-node scalar weight µsj though, an example will be provided. Consider, for instance, George S. Eleftheriou, Guillaume Pierrot an anisotropic stencil, depicted in figure 1b, that is in a critical condition, meaning that if it is allowed to further deform and squeeze, it risks degenerating to a line segment or distort. This sounds like exactly the case where it would be useful to favour rigidity in the direction of squeeze, freezing any further deformation. This is where µsj comes in handy serving as a limiter to further distortion. For instance, it could be something like   ||x~j − c~s ||2 exp − h2 (13) µsj = ||x~j − c~s ||2 + ε h being the average size of the mesh and ε an arbitrary positive constant just to avoid division by zero when the node under consideration coincides with the stencil’s center of gravity. The coefficient and its underlying logic will be detailed in section 4. The quantity e~sj T · e~sj (scalar), henceforth denoted as ke~sj k2 , will be shown (later) that relates to the anisotropic deformation energy of the stencil. 3 6 9 2 5 8 1 4 7 (a) Isotropic stencil 3 6 9 2 5 8 1 4 7 (b) Squeezed/anisotropic stencil (we favour rigidity in the direction of squeeze) Figure 1: Isotropic and anisotropic stencils George S. Eleftheriou, Guillaume Pierrot 3 TOTAL DISTORTION ENERGY At this point, a quantity called “total distortion energy” of the mesh will be introduced, to serve as the objective function and at the same time the distortion metric, bound to be minimized by the morphing algorithm. It is XX E(u~j , a~s , b~s ) = ke~sj k2 (14) s j∈s which is quadratic. Hence, to find the stationary points of the total distortion energy and minimize it we require for all nodes j and stencils s ∂E ∂E ∂E = = =0 ∂ υ~j ∂ a~s ∂ b~s 4 (15) ASYMPTOTIC DEVELOPMENT We will see in this section how to relate the distortion metric of a stencil S: X Es = µsj n o2 ~vj − α ~ s − β~s ∧ (~xj − ~cs )   j∈S (16) to the strain rate tensor: ε̇ = The first remark is that if (15) then we surely have: n 1 2 ∂v α ∂v β + ∂xβ ∂xα   (~vi )i∈II , α ~ s , β~s 1≤β≤d s∈S o is solution of the optimisation problem n o2   X ~ ~ µsj ~vj − ~as − bs ∧ (~xj − ~cs ) α ~ s , βs = argmin ~as ,~bs (17) 1≤α≤d (18) j∈S And subsequently: X j∈S µsj ! α ~s = X j∈S µsj n o ~vj − β~s ∧ (~xj − ~cs ) (19) and: Ys β~s = X j∈S with: µsj (~xj − ~cs ) ∧ (~vj − α ~ s) (20) George S. Eleftheriou, Guillaume Pierrot Ys = X j∈S µsj {~eδ ∧ (~xj − ~cs )} · {~eγ ∧ (~xj − ~cs )} !1≤γ≤d (21) 1≤δ≤d where (~eδ )1≤δ≤d denotes the canonical basis of Rd . Up to the change of variable: α ~s ← α ~ s + β~s ∧ (~xs − ~cs ) (22) where: equation (19) and (20) become: P xj j∈S µsj ~ xs = P j∈S µsj and: P vj j∈S µsj ~ α ~s = P j∈S µsj Ỹs β~s = X j∈S (23) (24) µsj (~xj − ~xs ) ∧ (~vj − α ~ s) (25) with: Ỹs = X j∈S µsj {~eδ ∧ (~xj − ~xs )} · {~eγ ∧ (~xj − ~xs )} !1≤γ≤d (26) 1≤δ≤d Let’s now assume that (~vj )j∈S is the interpolated of a smooth vector field ~v at nodes location: (~vj )j∈S = (~v (~xj ))j∈S We have the following Taylor expansion: ~v (~xj ) = ~v (~xs ) + D~v (~xj − ~xs ) + O(h2 ) ∀j ∈ S (27) where h is a characteristic size of the stencil and: D~v = is the vector field differential at ~xs .  ∂v α ∂xβ 1≤β≤d 1≤α≤d (28) George S. Eleftheriou, Guillaume Pierrot By the definition of ~xs (23) and (24) it comes then that: α ~ s = ~v (~xs ) + O(h2 ) (29) And (25) becomes: Ỹs β~s = X j∈S µsj (~xj − ~xs ) ∧ (D~v (~xj − ~xs )) + O(h3 ) (30) Now, by introducing the classical splitting of D~v between its symmetric and anti-symmetric part: D~v = ε̇ + Ω̇ (31) where: ε̇ =  1 2  ∂v α ∂v β − ∂xβ ∂xα 1≤β≤d (32) 1≤α≤d is the so-called rotation rate tensor, it comes: Ỹs (β~s − ω ~ s ) = O(h3 ) where (33) 1~ ω ~s = ∇ ∧ ~v 2 is the so-called rotation rate vector. Hence we have: β~s = ω ~ s + O(h) (34) Using now (27),(29) and (34), and substituting back in equation (16) it comes: Es = = X µsj j∈S X µsj j∈S  ~vj − ~v (~xs ) − ω ~ s ∧ (~xj − ~xs ) + O(h2 )  ε̇s (~xj − ~xs ) + O(h2 ) 2 2 = (Xs ε̇s ) : ε̇s + O(h3 ) (35) where : denotes the Frobenius product of 2 matrices and: Xs = X j∈S µsj (~xj − ~xs )α (~xj − ~xs )β !1≤β≤d 1≤α≤d (36) George S. Eleftheriou, Guillaume Pierrot is the so-called stencil anisotropy tensor. We will now consider 2 practical cases. The first one corresponds to a pure isotropic stencil (Fig. 1a). With the weights µsj being chosen as indicated in (13), it is easy to see that: ~cs = ~xs = ~x5 where ~x5 denotes the central node of the stencil and: Xs = µs Id with: µs )  (x1j − x15 )2 kx~j − x~5 k2 = exp − h2 ||x~j − x~5 ||2 + ε j∈S ) (   X (x2j − x25 )2 kx~j − x~5 k2 = exp − h2 ||x~j − x~5 ||2 + ε j∈S ( X  Consequently, there is no favoured direction on the stencil distortion metric. In case of a squeezed mesh, on the other hand (Fig. 1b), as it happens in boundary layers for instance, things are quite different. Indeed, we have again that ~cs = ~xs = ~x5 but the stencil anisotropy tensor is no longer scalar: Xs = DIAG (µαs )1≤α≤d with: µαs = ( X j∈S kx~j − x~5 k2 exp − h2  
2 ) xαj − xα5 ||x~j − x~5 ||2 + ε Because this time the stencil dimension in the y direction hy is way smaller than in the x direction hx and we have that: hx c(Xs ) = hy where c(Xs ) denotes the condition number of the matrix. Consequently it comes:  Es = hx c(Xs ) ε̇2yy + ε̇2xx + O(h3 ) and we see that we put more rigidity in the direction of squeezeness, which is is to be preferred as it is the direction of most imminent breakdown of the mesh. In this sense we can say that the Rigid Motion Mesh Morpher algorithm handles naturally mesh anistropy. George S. Eleftheriou, Guillaume Pierrot 5 BUILDING THE SYSTEM OF EQUATIONS The quadratic minimization problem of the total distortion energy of the mesh, as stated in (15), brings us to the following symmetric positive definite system       Pu Auu Au(a|b) u = · (37) A(a|b)u A(a|b)(a|b) P(a|b) (a|b) where the RHS consists of the boundary conditions, namely the prescribed nodes’ velocities. Attempting to solve it using the Schur complement leads to two different cases of elimination: Either u = f (a|b) (a|b) = g(u) or (38) So for instance −1 u = −A−1 uu · Au(a|b) · (a|b) + Auu · Pu → −1 (−A(a|b)u · A−1 (39) uu · Au(a|b) + A(a|b)(a|b) ) · (a|b) = −A(a|b) · Auu · Pu + P(a|b) Specializing the above analysis, for ∂ ∂E = ∂ υ~j ∂ υ~j XX s j∈s ke~sj k2 ! ∂E = 0 we have: ∂ υ~j ∂ = ∂ υ~j X j∈s1 ke~s1 j k2 + X j∈s2 ke~s2 j k2 + ... + X j∈sN kes~N j k2 ! (40) So for a stencil:  X ∂ke~s j k2  k  ∂ X ∂ υ ~ ke~sk j k2 = j j∈sk  ∂ υ~j j∈s  k 0 Considering the above, equation (40) becomes if j ∈ sk if j ∈ / sk X ∂ke~sj k2 X ∂E ∂ e~sj = = =0 2e~sj T ∂ υ~j ∂ υ~j ∂ υ~j s∋j s∋j But X s∋j 2e~sj T (41) X X  ∂ e~sj T ∂ e~sj T ∂ e~sj e~sj = 0 =0⇔ e~sj =0⇔ ∂ υ~j ∂ υ~j ∂ υ~j s∋j s∋j (42) (43) Since e~sj is linear with respect to υ~j , a~s , b~s , it could also be expressed as follows: e~sj = ∂ e~sj ∂ e~sj ~ ∂ e~sj υ~j + a~s + bs ∂ υ~j ∂ a~s ∂ b~s (44) Hence  X  ∂ e~sj T  ∂ e~sj ∂ e~sj ∂ e~sj ~ ∂E = υ~j + a~s + bs = 0 ∂ υ~j ∂ υ~j ∂ υ~j ∂ a~s ∂ b~s s∋j (45) George S. Eleftheriou, Guillaume Pierrot We denote by X  ∂ e~sj T ∂ e~sj Nj = ∂ υ~j ∂ υ~j s∋j  T   ∂ e~sj ∂ e~sj ∂ e~sj Hjs = ∂ a~s ∂ b~s ∂ υ~j the equation becomes Nj υ~j = − For ∂E = 0 we have ∂ a~s ∂ ∂E = ∂ a~s ∂ a~s XX s j∈s ke~sj k2 ! ∂ = ∂ a~s X s∋j X j∈s1 (46) (47)   a~s Hjs ~ bs ke~s1 j k2 + X j∈s2 (48) ke~s2 j k2 + ... + X j∈sN kes~N j k2 ! (49) depending on the stencil  X ∂ke~s j k2  k  ∂ X 2 ∂ a~s ke~sk j k = j∈sk  ∂ a~s j∈s  k 0 equation (49) becomes if s = sk (50) if s 6= sk X ∂ke~sj k2 ∂E = ∂ a~s ∂ a~s j∈s  X  ∂ e~sj T  ∂ e~sj ∂ e~sj ∂ e~sj ~ = υ~j + a~s + bs = 0 ⇒ ∂ a~s ∂ a~s ∂ υ~j ∂ a~s ∂ b~s j∈s j∈s ! ! X  ∂ e~sj T ∂ e~sj X  ∂ e~sj T ∂ e~sj X  ∂ e~sj T ∂ e~sj ~ bs = − a~s + υ~j ~s ∂ a ~ ∂ a ~ ∂ a ~ ∂ a ~ ∂ υ ~ s s s s j ∂ b j∈s j∈s j∈s (51) X ∂ke~sj k2 ∂E = 0 we end up with: ∂ b~s ! ! X  ∂ e~sj T ∂ e~sj X  ∂ e~sj T ∂ e~sj X  ∂ e~sj T ∂ e~sj a~s + υ~j b~s = − ~s ~s ~s ~s ∂ a ~ ∂ υ ~ s j ∂ b ∂ b ∂ b ∂ b j∈s j∈s j∈s (52) Similarly for Combining equations (52) (53) into a more compact form:    T ∂ e~sj ∂ e~sj  ∂ a~ ∂ υ~j    s   X X  a~s   υ~j = − Ms ~ = − HTjs υ~j   bs   j∈s   j∈s  ∂ e~sj T ∂ e~sj  ∂ υ~j ∂ b~s (53) (54) George S. Eleftheriou, Guillaume Pierrot Where  X  ∂ e~sj T ∂ e~sj  ∂ a~s  j∈s ∂ a~s   Ms =    X  ∂ e~sj T ∂ e~sj  ∂ a~s ∂ b~s j∈s is a symmetric 6x6 matrix. 5.1  X  ∂ e~sj T ∂ e~sj  ∂ a~s ∂ b~s  j∈s       T X ∂ e~sj ∂ e~sj   ~s ~s ∂ b ∂ b j∈s (55) eliminating the node velocities υ~j From (48) we have: υ~j = − X N−1 j Hjs s∋j   a~s b~s (56) and replacing υ~j in equation (54) (the two relations do not necessarily concern the same stencil, that’s why an s′ is introduced),  !   X X a~s′ a~s −1 T Nj Hjs′ ~ ⇒ Hjs − Ms ~ = − b s′ bs j∈s s′ ∋j  !   X X a~s′ a~s − HTjs N−1 ⇒ Ms ~ = − j Hjs′ ~s′ b bs ′ j∈s s ∋j  !   X X a~s′ a~s HTjs N−1 ⇒ Ms ~ = j Hjs′ b~s′ bs s′ j∈s∩s′       X X  a~s′  a~s −1 T  ′ Ms ~ − N H H js js j  b~s′ = 0  bs   ′ ′ s j∈s∩s (57) ass′ s6=s′ 5.2 eliminating stencil velocities a~s , b~s Starting from equation (54)   X a~s T = − M−1 ~j s Hjs υ ~ bs j∈s (58) Changing the node index to i (instead of j so far), since it’s not the same nodes we are referring to, equation (48) becomes   X a~s (59) Ni υ~i = − His ~ bs s∋i George S. Eleftheriou, Guillaume Pierrot   a~s Replacing ~ as it is expressed in (58) in the above equation yields bs     X X  −1 T   H M H Ni υ~i − is s js   υ~j = 0  j s∋{i,j} (60) aij i6=j 5.3 Elaborating on some expressions So far we referred generically to some expressions without getting into the details of computing them or making them more explicit. Finally, this subsection is detailing this aspect. Let us first consider the partial derivatives of e~sj : ∂ e~sj √ = ws µsj I(3×3) ∂ υ~j ∂ e~sj √ = − ws µsj I(3×3) ∂ a~s   0 cs3 − xj3 xj2 − cs2 ∂ e~sj √ 0 cs1 − xj1  = ws µsj xj3 − cs3 ~ ∂ bs cs2 − xj2 xj1 − cs1 0 It is worth noting that compact notation (61) (62) (63) ∂ e~sj ∂ e~sj ∂ e~sj , are diagonal and is antisymmetric. Using a more ∂ υ~j ∂ a~s ∂ b~s  ∂ e~sj ∂ b~   s −1 0 0 0 cs3 − xj3 xj2 − cs2 √ 0 cs1 − xj1  = ws µsj  0 −1 0 xj3 − cs3 0 0 −1 cs2 − xj2 xj1 − cs1 0  ∂ e~sj ∂ e~sj = ∂ a~s ∂(a~s , b~s ) (64) In light of the differentiation of e~sj the expressions of Nj and Hjs now become clearer  ∂ e~sj Hjs = ∂ υ~j T  ∂ e~sj ∂ a~s  ∂ e~sj ⇒ ∂ b~s  −1 0 0 0 cs3 − xj3 xj2 − cs2 0 cs1 − xj1  Hjs = ws µsj  0 −1 0 xj3 − cs3 0 0 −1 cs2 − xj2 xj1 − cs1 0  X  ∂ e~sj T ∂ e~sj Nj = ⇒ ∂ υ ~ ∂ υ ~ j j s∋j X Nj = ws µsj I(3×3) s∋j (65) (66) George S. Eleftheriou, Guillaume Pierrot We denote M = diag(Ms ) N = diag(Ni ) (67) (68) ≤ j ≤ nnodes H = (Hjs )11 ≤ s ≤ nstencil (69) As an example, eliminating u as shown in (39), would mean M = A(a|b)(a|b) H = Au(a|b) N = Auu (70) (71) (72) Thus, the coefficient matrices indexed according to the two different cases of elimination (38) are given: Au = N − HM−1 HT T −1 A(a|b) = M − H N H 6 (73) (74) CONCLUSIONS • In the present article, a novel method, Rigid Motion Mesh Morpher has been introduced. We provided indications, according to which, it handles intrinsically, at the same time, the rotations and the mesh anisotropy of any mesh, since its design is mesh-less. • It necessitates the solution of a sparse, symmetric positive definite (SPD) problem, proportional to the size of the mesh. • Ongoing work: we need to validate the approach on a variety of test cases and compare it to existing approaches with the hope to demonstrate its superiority. • It will be implemented as a part of an adjoint solver loop. • Future work may include smoothing of the surface node velocity field as provided by the adjoint solver or equivalently a CAD-free parametrization of the surface. 7 ACKNOWLEDGEMENTS This work was funded by the EU through the FP7-PEOPLE-2012-ITN ”AboutFlow” Grant agreement number:317006. REFERENCES [1] Carlo L. Botasso, David Detomi, and Roberto Serra. The Ball-Vertex Method: a New Simple Spring Analogy Method for Unstructured Dynamic Meshes. Computer Methods in Applied Mechanics and Engineering, 194:4244–4264, 2005. [2] Christoph Degand and Charbel Farhat. A three-dimensional torsional spring analogy method for unstructured dynamic meshes. Computers & structures, 80(3):305–316, 2002. George S. Eleftheriou, Guillaume Pierrot [3] Richard P. Dwight. Robust Mesh Deformation using the Linear Elasticity Equations. In Computational Fluid Dynamics 2006, pages 401–406. Springer, 2009. [4] Ch Farhat, C Degand, B Koobus, and M Lesoinne. Torsional springs for two-dimensional dynamic unstructured fluid meshes. Computer methods in applied mechanics and engineering, 163(1):231–245, 1998. [5] David A Field. Laplacian smoothing and Delaunay triangulations. Communications in applied numerical methods, 4(6):709–712, 1988. [6] Peter Hansbo. Generalized laplacian smoothing of unstructured grids. Communications in numerical methods in engineering, 11(5):455–464, 1995. [7] Stefan Jakobsson and Olivier Amoignon. Mesh deformation using radial basis functions for gradient-based aerodynamic shape optimization. Computers & Fluids, 36(6):1119– 1136, 2007. [8] AM Morris, CB Allen, and TCS Rendall. Cfd-based optimization of aerofoils using radial basis functions for domain element parameterization and mesh deformation. International journal for numerical methods in fluids, 58(8):827–860, 2008. [9] Nilanjan Mukherjee. A Hybrid, Variational 3D Smoother for Orphaned Shell Meshes. In IMR, pages 379–390. Citeseer, 2002. [10] Marco Andrea Pischedda. Mesh Morphing techniques based on continuum deformation models. Master’s thesis, Politecnico di Milano, 2008. [11] Z. Su, S. Wang, C. Yu, F. Liu, and X. Shi. A Novel Laplacian Based Method for Mesh Deformation. Journal of Information and Computational Science, 7(4):877–883, 2010. [12] Anh H Truong, Chad A Oldfield, and David W Zingg. Mesh movement for a discrete-adjoint Newton-Krylov algorithm for aerodynamic optimization. AIAA journal, 46(7):1695–1704, 2008. [13] Shaoting Zhang, Junzhou Huang, and Dimitris N Metaxas. Robust mesh editing using Laplacian coordinates. Graphical Models, 73(1):10–19, 2011.