* Zhuo
email: zhuoli@ee.tamu.edu

Weiping Shi
email: wshi@ee.tamu.edu

An Time Algorithm for Optimal Buffer Insertion with b Buffer Types

Buffer insertion is a popular technique to reduce the interconnect delay. The classic buffer insertion algorithm of van has complexity

where n is the number of buffer positions.

Cheng and Lin extended van Ginneken's algorithm to allow b buffer types in time For modern design libraries that contain hundreds of buffers, it is a serious challenge to balance the speed and performance of the buffer insertion algorithm. In this paper; we present a new algorithm that computes the optimal buffer insertion in time. The reduction is achieved by the observation that the C )pairs of the candidates that generate the new candidates must form a hull. On industrial test cases, the new algorithm is faster than the previous best buffer insertion algorithms by orders of magnitude.

Introduction

Delay optimization techniques for interconnect are increasingly important for achieving timing closure of high performance designs. One popular technique for reducing interconnect delay is buffer insertion. A recent study by Saxena et al.

projects that 35% of all cells will be intra-block repeaters for the 45 node. Consequently, algorithms that can efficiently insert buffers are essential for the design automation tools.

In 1990, van proposed an optimal buffer insertion algorithm for one buffer type. His algorithm has time complexity

where n is the number of candidate buffer positions. Lillis, Cheng and Lin [START_REF] Lillis | Optimal wire sizing and buffer insertion for low power and a generalized delay model[END_REF] extended van Ginneken's algorithm to allow buffer types in time O(n n)for multi-pin nets, for one buffer type. Several works have built upon van Ginneken's algorithmand its extension for multiple buffer types to include wire sizing simultaneous tree construction [8, 6, 5, 9, noise constraints [2] and resource minimization [7, Modern design libraries may contain hundreds of different buffers with different input capacitances, driving resistance, intrinsic delay, power level, etc. If every buffer available for the given technology is supplied, it is stated in [START_REF] Alpert | Buffer library selection[END_REF] that the current algorithms could possibly take days or even weeks for large designs since all these algorithms are quadratic in terms of b. Alpert et. [START_REF] Alpert | Buffer library selection[END_REF] studied how to reduce the size of the buffer library with a clustering algorithm. Though the buffer library size is reduced, the solution quality is often degraded accordingly.

In this paper, we propose a new algorithm that performs optimal buffer insertion with buffer types in time. Our speedup is achieved by the observation that the candidates that generate new buffered candidates must lie on the convex hull of (Q,C). Experimental results show that our algorithm is significantly faster than previous best algorithms.

Section 2 formulates the problem. Section 3 describes the new algorithm. Simulation results are given in Section4 and conclusions are given in Section 5.

Preliminary

A net is given as a routing tree T = (V,E ) , where V = and E V x V. Vertex is the source vertex and also the root of T , is the set vertices, and is the set of internal vertices.

New Algorithm

The previous best algorithm for multiple buffer types by Lillis, Cheng and Lin consists of three major operations: 1) adding buffers at a buffer position in time, 2) adding a wire in time, and 3) merging two branches in + time, where and are the numbers of buffer positions in the two branches. As a result, their algorithm has time complexity

In this section, we show that the time complexity of the first operation, addingbuffers at a buffer position, can be reducedto and thus our algorithm can achieve total time complexity Assume we have computed the set of nonredundant candidates for and now reach a buffer position see Fig. 1. Wire (v, has 0 resistance and capacitance. Define as the slack if we add a buffer type at v for any candidate a in If we do not insert any buffer at then every candidate for is a candidate for If we insert a buffer at then for every buffer type Proof: This procedure is known as Graham's scan in computational geometry It finds the convex hull of a set of points in sorted order in linear time.

It is well known that a set of points form a convex hull if and only if there are no consecutive a and that satisfy Eq. (2). Therefore, Convexpruning is correct since it checks all consecutive candidates.

To analyze the time complexity, consider the number of forward and backward moves. Each time moves backward, it deletes a candidate. Therefore, there can be at most k backward moves.

The number of forwardmoves is the size of the list plus the number of backward moves. Therefore the number of forward moves is at most 2k. Hence the time complexity is Proof: Since are in the nondecreasing order of capacitance and the given set of nonredundant candidates are in nondecreasing order of it takes

+

Since the operation of adding a buffer is reduced to time from Theorem 1 and 2, it is easy to see that buffer insertion with b buffer types can be done in worst case time with our new algorithm.

the index in time.

= time to merge the two sorted lists.

Simulation

Both the algorithm of Lillis et al. [START_REF] Lillis | Optimal wire sizing and buffer insertion for low power and a generalized delay model[END_REF] and the new algorithm are implemented in C and runon a Sun SPARC workstations with 400 and 2 GB memory. The device and interconnect parameters are based on TSMC 180 nm technology. We have 4 different buffer libraries, with the size 8, 16, 32 and 64 respectively. is chosen from 180 to 7000 Q, is chosen from 0.7 to 23 and is chosen from 29 ps to 36.4 ps. The sink capacitances range from 2 to 41 The wire resistance is 0.076 and the wire capacitance is 0.118 Table shows for large industrial circuits, the new algorithm is up to 11times faster than Lillis'

The memory usage is not shown in the table, but there is only almost 2% memory overhead due to the double linked list used by the new algorithm. When b is small, algorithm has a little time overhead compared to Lillis' algorithm. due to function Fig. 3 shows the time complexity curve of two algorithms for the net with 1944 sinks and 33133 buffer positions with respect to the size of buffer library b. In the fig-ure, the y axis is normalized to the running time of the case when the buffer library size is 8. Though the worst case time complexity of Lillis' algorithm is quadratic in terms of b, it behaves more like a linear function of as observed in

The time complexity curve of our algorithm is also linear, but has a much smaller slope. of the two algorithms for the net with 1944 sinks, with respect to the number of buffer positions n. The buffer library size is 32. In the figure, the y axis is normalized to the running time of the case with 1943 buffer positions. We can see that while Lillis' and our algorithms both behave quadratically, our algorithm shows much slower growing trend since the operation of adding a buffer becomes more dominant among three major operations when n increases.

Conclusion

We presented a new algorithm for optimal buffer insertion with buffer types of worst case time This is an improvement of the previous best algorithm Simulation results show our new algorithm is significantly faster than algorithms for large industrial circuits with large buffer libraries. Our algorithm can also be applied to reduce buffer cost. We leave the details to the journal version.

=

  Figure 1. consists of buffer position and

Figure 2 .

 2 Figure 2. (a) Nonredundant candidates on (Q,C) plane. (b) Nonredundantcandidates after convex pruning.

>

  Accordingto the definition of convex pruning, is not pruned.Lemma4Let the set of nonredundant candidates ter Convexpruning be and assume are sorted in increasingQ and order. Considerany three candidates a, = 1,.. . ,b. Nonredundant candidates in are stored in increasing order using a double link list pointed by header.Buffer types are sorted in non-increasing driver resistance order and stored in array B. Function P a) computes as defined in Eq. (1).

  with zero resistance and capacitance,nonredundantcandidates of N are stored in increasingQ and order, then function AddBuff er generates all new candidates in time. Proof: Let the set of nonredundant candidates after Convexpruning be From Lemma 3, we know that all best candidates are in From Lemma and Lemma 4, starting from the first candidates in function AddBuf fer can find all in the increasing order of i. Now consider the time complexity. Function Convexpruning takes time according to Lemma 2. The for loop takes + b) = time. It takes only time to sort the entire buffer library in terms of the input capacitance and establish the order from buffer index i to the order in Each time function AddBuf f er is called, the new candidates can be sorted in nondecreasing order by using Theorem 2 Given a set nonredundantcandidates sorted in increasing Q and all b new candidates can be inserted in time.

Figure 3 .

 3 Figure 3. Comparison of normalized running time with respect to buffer library size among two algorithms. Number of sink is and number of buffer positions is 33133.

Fig. 4

 4 Fig.4shows the time complexity of the two algorithms for the net with 1944 sinks, with respect to the number of buffer positions n. The buffer library size is 32. In the figure, the y axis is normalized to the running time of the case with 1943 buffer positions. We can see that while Lillis' and our algorithms both behave quadratically, our algorithm shows much slower growing trend since the operation of adding a buffer becomes more dominant among three major operations when n increases.

ReferencesC.

  Alpert and A. Devgan. Wire segmenting for improved buffer insertion. In DAC, pages 588-593, 1997. [2] C. J. Alpert, A. Devgan, and S. T. Quay. Buffer insertion for noise and delay optimization. In DAC, pages 362-367, 1998.

Figure 4 .

 4 Figure 4. Comparison of normalized running time with respect to buffer positions among two algorithms. Number of sink is and number of buffer types is 32.

  tion. IEEE Trans. CAD, to appear. W. Shi and Z. Li. An time algorithm for optimal buffer insertion. In pages 580-585,2003. W. Shi, Li, and C. J. Alpert. Complexity analysis and speedup techniques for optimal buffer insertion with minimum cost. In ASPDAC, pages 609414,2004. [ L. P. P. van Buffer placement in distributed tree network for minimal delay. In ISCAS, pages 868, 1990. H. Zhou, D. F. Wong, I. M. Liu, and A. Aziz. Simultaneous routing and buffer insertion with restrictions on buffer locations. IEEE Trans. CAD,

slack of v under a is a ) = min - a ) } . Buffer Insertion Problem: Given routing tree T = (V,E ) ,sink capacitance and for each sink capacitance and resistance for each edge e , possible buffer position and buffer library find a can- didate a for T that maximizes

	each buffer type	the intrinsic delay is	driv-
	ing resistance is	and input capacitance is	A
	function f :		specifies the types of buffers al-
	lowed at each internal vertex. Each edge e E is associ-
	ated with lumped resistance	and capacitance
	Following previous researchers [14, 7, 9, 15,	we use
	the Elmore delay for the interconnect and the linear delay
	for buffers. For each edge e =	signals travel from
	to	The	delay of is
	where		is the downstream capacitance at	For any
	buffer type	at vertex	the buffer delay is
				.	+
	where		is the downstream capacitance at	When a
	buffer is inserted, the capacitance viewed from the upper
	stream is
	For any vertex	V , let	be the	down-
	stream from and with being the root. Once we decide
	where to insert buffers in	we have a candidate a for
	The delay from to sink	under a is
					a ) .
	The effect of a candidate to the upstream is described
	by slack Q and downstream capacitance	Define
		as the downstream capacitance at node under
	candidate a. For any two candidates	and	of
	we say	dominates	if	and
			The set of nonredundant candidates
	of	which we denote as	is the set of candidates
	such that no candidate in	dominates any other candi-
	date in		and every candidate of	is dominated by
	some candidates in	Once we have	the can-
	didate that gives the maximum

Each sink vertex is associated with sink capacitance and required arrival time A buffer library contains differenttypes of buffers and its size is representedby b. For = + where the sum is over all edges e in the path from to s. If is a buffer in then is the buffer delay. If is not a buffer in a,then = The a ) can be found easily. The number of total nonredundant candidates is at most 1for one buffer type and 1for b buffer types [ where n is the number of candidatebuffer positions.

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

This research was supported by the NSF grants CCR-0098329, 01 13668, 512-0266-2001.


Shuvam Sen

Guillaume De Nayer

Michael Breuer
email: breuer@hsu-hh.de

A fast and robust hybrid method for block-structured mesh deformation with emphasis on FSI-LES applications

Keywords: Mesh deformation, Radial basis function, Inverse distance weighting function, Greedy algorithm, Transfinite interpolation, Fluid-structure interaction, Large-eddy simulation

The present work introduces an efficient technique for the deformation of block-structured grids occurring in simulations of fluid-structure interaction (FSI) problems relying on large-eddy simulations (LES). The proposed hybrid approach combines the advantages of the inverse distance weighting (IDW) interpolation with the simplicity and low computational effort of transfinite interpolations (TFI), while preserving the mesh quality in boundary layers. It is an improvement over the state-of-the-art currently in use. To reach this objective, in a first step three elementary mesh deformation methods (TFI, IDW and radial basis functions) are investigated based on several test cases of different complexities analyzing their capabilities but also their computational costs. That allows to point out the advantages of each method but also demonstrates their drawbacks. Based on these specific properties of the different methods, a hybrid methodology is suggested which splits the entire grid deformation into two steps: First the movement of the block-boundaries of the block-structured grid and second the deformation of each block of the grid. Both steps rely on different methodologies, which allows to work out the most appropriate method for each step leading to a reasonable compromise between the grid quality achieved and the computational effort required. Finally, a hybrid IDW-TFI methodology is suggested which best fits to the specific requirements of coupled FSI-LES applications. This hybrid procedure is then applied to a real-life FSI-LES case.

INTRODUCTION

To compile an optimal methodology for the simulation of dynamically coupled fluid-structure interactions (FSI) of lightweight flexible structures in turbulent flow, a robust and efficient method combining eddy-resolving methods such as large-eddy simulation (LES) with flexible structures was developed and validated [START_REF] Breuer | Fluid-structure interaction using a partitioned semiimplicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF]. It relies on a partitioned Arbitrary Lagrangian-Eulerian (ALE) method applying block-structured grids. The FSI-LES applications tested were three-dimensional, however, with nominally two-dimensional geometries [START_REF] De Nayer | Flow past a cylinder with a flexible splitter plate: a complementary experimental-numerical investigation and a new FSI test case (FSI-PfS-1a)[END_REF][START_REF] De Nayer | Numerical FSI investigation based on LES: Flow past a cylinder with a flexible splitter plate involving large deformations (FSI-PfS-2a)[END_REF]. In these FSI-LES cases a simple mesh deformation method based on transfinite interpolation (TFI) was sufficient to obtain a reasonable quality of the deformed mesh. However, for practically relevant FSI applications with fully three-dimensional geometries and large deformations, this mesh deformation method is not usable anymore. Moreover, improvements of the FSI-LES method have been achieved concerning S. SEN ET AL. updated by the disk relaxation algorithm using repulsive and attractive interactions to give rise to a high-quality mesh.

Apart from the above-mentioned methods, hybrid procedures combining different techniques are also available in the literature. The so-called moving sub-mesh approach (MSA) is a mixture between a physical analogy-based pseudo-structure method carried out on a coarser grid and an interpolation from this coarse grid onto the original grid [START_REF] Lefranc | A simple mesh deformation technique for fluid-structure interaction based on a submesh approach[END_REF]. A similar idea based on a spring method was successfully applied by Xuan et al. [START_REF] Xuan | Spring-interpolation approach for generating unstructured dynamic meshes[END_REF]. Another hybrid algorithm proposed by Landry et al. [START_REF] Landry | Robust mesh-mover algorithms for hybrid meshes[END_REF] moves the coarse mesh by the IDW method and then interpolates it to the fine mesh using a linear interpolation. In general, hybrid methods have the advantage that an optimal compromise between computational speed and quality of the grid can be found.

In the present work a new hybrid method is introduced with special emphasis on FSI-LES applications. The aim is to arrive at an appropriate grid deformation algorithm, which is the best compromise between high mesh quality and low CPU-time for block-structured grids with orthogonal cells in the vicinity of walls. To reach this goal, three elementary mesh deformation methods are selected, evaluated and then combined. TFI is chosen for its fast speed and its high quality of the deformed mesh obtained in a single block. RBF is taken into account, because of its reputation of robustness and the very high grid quality, which can be achieved. Finally, IDW is also selected for its simplicity, speed and good control of the mesh deformation in the near-wall region preserving the grid orthogonality.

The paper is organized as follows: Section 2 mathematically discusses the TFI, RBF and IDW methods. Section 3 deals with the numerical evaluation and comparison of these elementary methods. In Section 4, the hybrid algorithm is proposed. Section 5 is concerned with the numerical validation of the proposed hybrid algorithm based on a real-life FSI-LES application.

MATHEMATICAL DESCRIPTION AND IMPLEMENTATION

Transfinite interpolation (TFI)

The transfinite interpolation was introduced by Gordon and Hall [START_REF] Gordon | Construction of curvilinear coordinate systems and applications to mesh generation[END_REF]. They used the term "transfinite" to describe the general class of interpolation methods, which match the original primitive function to be interpolated at a non-denumerable, i.e., transfinite number of points. Later Eriksson [START_REF] Eriksson | Generation of boundary-conforming grids around wing-body configurations using transfinite interpolation[END_REF] used TFI for grid generation in CFD. When used for this purpose, TFI has the advantage to provide complete conformity to surface boundaries of the domain to be meshed. Mathematically TFI is the Boolean sum of the tensor products formed using univariate interpolations. In general, these univariate interpolation functions are restricted to the linear Lagrangian combination of surface displacements based on the arc-length from the boundary points to the node.

For a known surface displacement, the mesh deformation d i,j,k = (d x i,j,k , d y i,j,k , d z i,j,k ) at an arbitrary volume point (i, j, k), 1 ≤ i ≤ I, 1 ≤ j ≤ J, 1 ≤ k ≤ K in three dimensions is given by

d i,j,k = F 1 ⊕ F 2 ⊕ F 3 = F 1 + F 2 + F 3 -F 1 F 2 -F 2 F 3 -F 3 F 1 + F 1 F 2 F 3 , (1)
where F i are univariate projectors. For two-dimensional cases the interpolation is only done along two directions and thus d i,j = (d x i,j , d y i,j ) reads

d i,j = F 1 ⊕ F 2 = F 1 + F 2 -F 1 F 2 . (2)
The linear interpolation and the tensor products of univariates are given by [START_REF] Thompson | Numerical Grid Generation: Foundations and Applications[END_REF][START_REF] Smith | Transfinite interpolation (TFI) generation systems. Handbook of Grid Generation[END_REF]:

F 1 = (1 -ξ i,j,k ) d 1,j,k + ξ i,j,k d I,j,k , (3)
F 2 = (1 -η i,j,k ) d i,1,k + η i,j,k d i,J,k , (4)
F 3 = (1 -ζ i,j,k ) d i,j,1 + ζ i,j,k d i,j,K , (5)
F 1 F 2 = (1 -ξ i,j,k )(1 -η i,j,k ) d 1,1,k + ξ i,j,k (1 -η i,j,k ) d I,1,k +(1 -ξ i,j,k )η i,j,k d 1,J,k + ξ i,j,k η i,j,k d I,J,k , (6)
F 2 F 3 = (1 -η i,j,k )(1 -ζ i,j,k ) d i,1,1 + η i,j,k (1 -ζ i,j,k ) d i,J,1 +(1 -η i,j,k )ζ i,j,k d i,1,K + η i,j,k ζ i,j,k d i,J,K , (7)
F 3 F 1 = (1 -ζ i,j,k )(1 -ξ i,j,k ) d 1,j,1 + ζ i,j,k (1 -ξ i,j,k ) d 1,j,K +(1 -ζ i,j,k )ξ i,j,k d I,j,1 + ζ i,j,k ξ i,j,k d I,j,K , (8)
F 1 F 2 F 3 = (1 -ξ i,j,k )(1 -η i,j,k )(1 -ζ i,j,k ) d 1,1,1 + ξ i,j,k (1 -η i,j,k )(1 -ζ i,j,k ) d I,1,1 +(1 -ξ i,j,k )η i,j,k (1 -ζ i,j,k ) d 1,J,1 + ξ i,j,k η i,j,k (1 -ζ i,j,k ) d I,J,1 +(1 -ξ i,j,k )(1 -η i,j,k )ζ i,j,k d 1,1,K + ξ i,j,k (1 -η i,j,k )ζ i,j,k d I,1,K
+(1ξ i,j,k )η i,j,k ζ i,j,k d 1,J,K + ξ i,j,k η i,j,k ζ i,j,k d I,J,K .

The arc-length based computational coordinates (ξ, η, ζ) used to maintain the grid distribution laws of the initial grid are

ξ 1,j,k =0, η i,1,k =0, ζ i,j,1 =0, ξ i,j,k = i n=2 r n,j,k -r n-1,j,k I n=2 r n,j,k -r n-1,j,k ,η i,j,k = j n=2 r i,n,k -r i,n-1,k J n=2 r i,n,k -r i,n-1,k ,ζ i,j,k = k n=2 r i,j,n -r i,j,n-1 K n=2 r i,j,n -r i,j,n-1 , (10)
where r i,j,k denotes the coordinate at (i, j, k) and . describes the Euclidian norm.

In order to reduce the CPU-time consumption of a CFD simulation, the coefficients (ξ i,j,k , η i,j,k , ζ i,j,k ) can be computed and stored at the beginning. Indeed, they are only based on the initially undeformed grid and do not vary during the time marching. TFI is a very fast method since its computational complexity linearly scales with the number of volume grid points N v . However, it requires a certain amount of memory for storing all coefficients.

Radial basis function (RBF) interpolation

Applying RBFs displacements at the coordinate points denoted by d(r) are interpolated using the following function

d(r) = N s n=1 α n φ( r -r s n ) + P (r). (11)
Here, r s n are the surface nodes at which the displacements are known, φ is the form of the RBF being used and P (r) is a set of three polynomials. Following the work of De Boer et al. [START_REF] Boer | Mesh deformation based on radial basis function interpolation[END_REF] P (r) is restricted to linear polynomials so that rigid-body displacements can be recovered. Given a surface displacement vector d s n , n ∈ {1, 2, ..., N s } for the three-dimensional case, three different sets of coefficients α n appearing in Eq. ( 11) and three linear polynomials of the form P = β 0 + β 1 x + β 2 y + β 3 z can be obtained by solving the linear symmetric indefinite system

            d s 1 d s 2 . . . d s N s 0 0 0 0             =             φ s 1 s 1 φ s 1 s 2 . . . φ s 1 s N s 1 x s 1 y s 1 z s 1 φ s 2 s 1 φ s 2 s 2 . . . φ s 2 s N s 1 x s 2 y s 2 z s 2 . . . . . . . . . . . . . . . . . . . . . . . . φ s N s s 1 φ s N s s 2 . . . φ s N s s N s 1 x s N s y s N s z s N s 1 1 . . . 1 0 0 0 0 x s 1 x s 2 . . . x s N s 0 0 0 0 y s 1 y s 2 . . . y s N s 0 0 0 0 z s 1 z s 2 . . . z s N s 0 0 0 0                         α 1 α 2 . . . α N s β 0 β 1 β 2 β 3             . ( 12
)
The first N s vector equations of the above system correspond to N s surface nodes whose displacements are known. The additional four vector equations originate from the compatibility S. SEN ET AL. requirement of a zero displacement vector at the far field given by

N s n=1 α n f (r s n ) = 0 (13)
with f (r) = {1, x, y, z}. Hence, the system finally is of the order (N s + 4) × (N s + 4) with φ s i s j = φ( r s i -r s j ). Once the values of α n and P are obtained, Eq. ( 11) can be used to estimate displacements at the interior nodal points.

In the literature various RBFs have been used to deform a mesh [START_REF] Boer | Mesh deformation based on radial basis function interpolation[END_REF][START_REF] Rendall | Unified fluid-structure interpolation and mesh motion using radial basis functions[END_REF][START_REF] Wang | Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions[END_REF][START_REF] Rendall | Evaluation of RBFs for volume data interpolation in CFD[END_REF][START_REF] Bos | Radial basis function based mesh deformation applied to simulation of flow around flapping wings[END_REF]. These basis functions can be divided into two categories: (a) RBFs with compact support and (b) RBFs with global support. It has been noticed that the approximation properties of RBFs with global support are generally better than those with compact support [START_REF] Boer | Mesh deformation based on radial basis function interpolation[END_REF][START_REF] Wang | Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions[END_REF]. On the other hand, RBFs with compact support lead to a well-conditioned sparse system and are computationally more efficient [START_REF] Rendall | Unified fluid-structure interpolation and mesh motion using radial basis functions[END_REF][START_REF] Wang | Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions[END_REF][START_REF] Rendall | Evaluation of RBFs for volume data interpolation in CFD[END_REF]. Some popular RBFs are listed in Table I for reference. Note that RBFs with compact support are scaled with a support radius R and are set to zero outside R.

Name

Abbreviation φ(r) Support Thin plate spline TPS

r 2 ln( r ) Global Wendland's C 0 W C 0 (1-r ) 2 Compact Wendland's C 2 W C 2 (1-r ) 4 (4 r +1) Compact Continuous TPS C 1 CTPS C 1 1 + 80 3 r 2 -40 r 3 +15 r 4 Compact -8 3 r 5 +20 r 2 ln( r ) Continuous TPS C 2 a CTPS C 2 a 1 -30 r 2 -10 r 3 +45 r 4 Compact -6 r 5 -60 r 3 ln( r ) Continuous TPS C 2 b CTPS C 2 b 1 -20 r 2 +80 r 3 -45 r 4 Compact -16 r 5 +60 r 4 ln( r )
For solving the resulting system of equations an iterative Krylov subspace based solver called minimum residual (MINRES) [START_REF] Paige | Solution of sparse indefinite systems of linear equations[END_REF] is preferred. MINRES is particularly suitable since the associated system of equations is symmetric. The system is assumed to be converged if the residual falls below a predefined tolerance ε c . Note that for small problems (see Section 3.1) a direct solver such as LU decomposition is advantageous in terms of required CPU-time. Furthermore, no dependency on the matrix characteristics of different RBFs is observed. However, for practical applications involving large systems of equations, direct solvers become unfeasible.

The solution of the system of equations ( 12) is computationally intensive. To reduce the computational effort Rendall and Allen [START_REF] Rendall | Efficient mesh motion using radial basis functions with data reduction algorithms[END_REF] introduced a data reduction strategy called the greedy algorithm. They discussed three different implementation perspectives for this algorithm and were able to significantly reduce the computational cost. The single point greedy algorithm [START_REF] Rendall | Efficient mesh motion using radial basis functions with data reduction algorithms[END_REF] starts with a single boundary point randomly chosen and calculating the solution of this system based on prescribed displacements. It then calculates the largest error between the exact and the approximate system. Based on this largest error and a predefined tolerance ε g the algorithm updates the coefficients of the reduced system, which may lead to the addition of further points in the already chosen list of so-called greedy points. Note that a greedy point is added when the corresponding coefficient is updated for the first time. However, the main deficiency of the greedy algorithm is that it fails to take the consequence of correcting a single equation of the system on the remaining equations of this system into account.

The greedy algorithm was speed up by Wang et al. [START_REF] Wang | Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions[END_REF] using a double-edged strategy, where a second search for the largest error is required in an approximately conjugate direction. In the present work a hybrid greedy algorithm [START_REF] Rendall | Efficient mesh motion using radial basis functions with data reduction algorithms[END_REF] is implemented. This algorithm starts by carrying out a greedy point selection procedure for m times, where m is some predefined number. As noted by Rendall and Allen [START_REF] Rendall | Efficient mesh motion using radial basis functions with data reduction algorithms[END_REF] the greedy algorithm visits a point many times and hence the number of executions is not equal to the total number of points chosen. Next the system consisting of so-called active greedy points is inverted. The process continues iteratively until the maximum absolute error reduces to the predefined tolerance ε g .

Here it should be noted that the greedy algorithm does not work with some global RBFs, in particular the thin plate spline (TPS). A point reduction algorithm for TPS, if it exists, is not known.

Inverse distance weighting (IDW) interpolation

IDW interpolation, first introduced by Shepard [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF], was adopted for mesh deformation by Witteveen [START_REF] Witteveen | Explicit and robust inverse distance weighting mesh deformation for CFD[END_REF]. In its generality this method interpolates a function at a given point r based on the weighted average of known values of scattered data points. Witteveen proposed that for a given boundary displacement vector d s n , n ∈ {1, 2, ..., N s }. The displacements at the interior nodes can then be calculated as

d(r) = N s n=1 w n (r) d s n N s n=1 w n (r) , (14)
where w n (r) is the weighting function often expressed as

w n (r) = r -r s n -p . ( 15
)
Witteveen used different values of p ∈ {1, 2, 3, 4} depending on the nature of the boundary node, i.e., whether the boundary node is dynamic or static as well as translating or rotating.

Luke et al. [START_REF] Luke | A fast mesh deformation method using explicit interpolation[END_REF] modified the IDW method by applying the improved weighting function

w n (r) = A n L def r -r s n a + αL def r -r s n b , (16)
where A n is the area weight of the boundary node n. They also suggested the best values for the parameters a and b as 3 and 5, respectively. L def and α are defined as

L def = N s max n=1 r s n -r centroid where r centroid = 1 N s + N v N s +N v m=1 r m , (17)
α = η L def N s max n=1 d s n (r) -d mean with d mean = N s n=1 a n d s n (r) and a n = A n N s m=1 A m . (18)
Following Luke et al. [START_REF] Luke | A fast mesh deformation method using explicit interpolation[END_REF] the value of the parameter is set to η = 5. L def denotes the maximum distance of any mesh node from the mesh centroid and α is a parameter, which provides the relative weight of the nearby nodes in comparison to the distant ones. Equation ( 18) is only one option and the user is free to set a fixed value of α based on his need.

Furthermore, Luke et al. [START_REF] Luke | A fast mesh deformation method using explicit interpolation[END_REF] improves the quality of their deformed mesh by interpolating not the displacements, but the rotations and the translations resulting from the splitting of the displacements. However, it is worthwhile to point out that the IDW method loses its extreme conserving property once rotations are included. Denoting the introduced rotation matrix M s and the associated translation vector as t s for a surface boundary point having the coordinate r in the original mesh, the displacement field d s (r) is defined as

r + d s (r) = M s r + t s . ( 19
)
Following the work of Maruyama et al. [START_REF] Maruyama | High-quality mesh deformation using quaternions for orthogonality preservation[END_REF] an ideal choice to estimate the rotation matrix is to define quaternions at each boundary mesh point. In terms of a quaternion Q defined at a boundary node, Eq. ( 19) can be rewritten as

R d = Q R u Q * + T , (20)
S. SEN ET AL.

where R u = [0, r], R d = [0, r + d] and T = [0, t] are the quaternions associated with the undeformed position vector, the deformed position vector and the translation vector, respectively. The subscripts u and d are used to indicate the undeformed and deformed states, respectively. The subscript s used to denote surface boundary points is dropped for the sake of brevity.

Q and T at each boundary node can be determined by using the steps outlined by Maruyama et al. [START_REF] Maruyama | High-quality mesh deformation using quaternions for orthogonality preservation[END_REF] and briefly reproduced in Appendix A. Note that T is coordinate dependent since the rotation is calculated with respect to a local coordinate system (see example in Appendix B). To eliminate this coordinate dependence, a strategy called transformation division algorithm was proposed by Maruyama et al. [START_REF] Maruyama | High-quality mesh deformation using quaternions for orthogonality preservation[END_REF] (see Appendix A). For large deformations this transformation division process is an absolute necessity as seen in Appendix B. However, if the IDW procedure relies on the last deformed mesh and the new displacements at the boundary nodes are small, this procedure can be omitted (see Appendix B). For example, in case of unsteady FSI-LES computations, very small time steps are typically used and thus marginal deformations appear at the boundaries of the structure between two time steps, so that the transformation division steps can be skipped.

Once the rotation quaternions and the translation vectors are known at the surface nodes, they can be interpolated to the volume nodes by using the IDW interpolation methodology given by Eq. ( 14). Note that the components of the rotation quaternions at the volume grid points have to be normalized before finally estimating the displacement vectors.

In the context of block-structured grids and a parallelization strategy based on domain decomposition, the application of IDW is done in the following steps: At the beginning of the simulation all fixed and moving boundary points are marked and saved in two different arrays for optimization purpose. These arrays have to be known by all processors. Then, at each time step the quaternions and translations of all boundary points N s are evaluated as explained above and shared by all processors. Finally, each processor executes a loop over all volume grid points of its own block evaluating the displacements based on Eq. ( 14). This explains the fact that the computational effort of IDW scales according to N s × N v , since the sum of all volume grid points contained in all blocks is N v .

EVALUATION OF THE ELEMENTARY METHODS BASED ON TEST CASES

In this section the elementary mesh deformation methods described above are evaluated on a few test cases of rising complexity. To provide a global measure of the mesh quality, the harmonic mean value

f skew =      1 N cv N cv -1 k=0 1 f skew k -1 if 0 < f skew k ≤ 1 for all k = 0, 1, ..., N cv 0 if f skew k = 0 for any k = 0, 1, ..., N cv
is defined following the work of Knupp [START_REF] Knupp | Algebraic mesh quality metrics for unstructured initial meshes[END_REF] where f skew k is the skew quality metric of the k-th control volume of the mesh having N cv control volumes. f skew = 1 will denote perfect orthogonality of the mesh.

For comparison of the computational effort of the different methods the floating-point operations and CPU-times are measured on the same empty machine with identical compiler options. The floating-point operations are evaluated by the open-source tool likwid † [START_REF] Treibig | Likwid: A lightweight performance-oriented tool suite for x86 multicore environments[END_REF] and given in Giga-FLOP (GFLOP). The CPU-times reported below are scaled by the CPU-time of the fastest successful method leading to relative computational costs nearly independent of the machine. 1 a square with the center at (0.5, 0.5) and a prescribed deformation on all sides of the domain is considered. Obviously, a sharp angle is generated at the bottom. On the other three sides the deformations are according to trigonometric functions in order to vary the angles and the overall shape. On the left and right boundaries identical displacements given by 0.1 sin (2πy) are applied, whereas on the top a deformation according to 0.2 sin (πx) is used. Moreover, on all boundary points a unit translation along the x-direction is prescribed. The original grid consisting of 21 × 21 grid points is clustered at the bottom and left wall. The undeformed grid is shown in Fig. 1(a). Computations are carried out for this coarse grid and a finer one with 301 × 301 grid points. The main objective is to evaluate the performance of TFI, RBF and IDW in a situation, which involves both contraction and relaxation. Also of interest is to see the effect on grid clustering. Furthermore, the orthogonality preserving nature of the schemes and its effect on the interior nodes is analyzed.

3.1.2. TFI mesh deformation For this type of deformation involving a convex undeformed domain, TFI is expected to perform well. The deformed grid obtained by TFI is depicted in Fig. 1(b). The color contours show the skew quality metric of the grid f skew . Values of unity define the optimum and decreasing values represent a deterioration of the orthogonality. Apparently, the translation is recovered with ease. However, the orthogonality is deteriorated, which is particularly evident in the lower left portion. Although high skew metric values are preserved at a large portion of the interior nodes, overall it is clear that TFI does not maintain the high level of orthogonality of the original grid. Nevertheless, the clustering of the undeformed mesh is preserved. A benefit of the extreme preserving property of TFI is that all grid points are retained inside the boundary without any cross-over. The main advantage of TFI is its very low computational effort (see Table II). Table II. Problem 1: Comparison of the floating-point operations (GFLOP), the relative computational costs and global measure of the mesh quality f skew . All computations done using the fine grid with 301 × 301 grid points.

Method

RBF mesh deformation

Since the main interest is to compare the number of floating-point operations and the computational costs of various RBFs having both global and local support, a uniform radius of influence of R = 1.5 is chosen for all RBFs with local support. This choice covers the entire domain. For all cases considered, convergence is assumed when the residual falls below ε c = 10 -12 within the MINRES solver. The grid deformed by Wendland's C 0 (W C 0 ) function is depicted in Fig. 2(a), while Fig. 2(d) uses the same basis function but with the hybrid greedy algorithm. Similarly, Figs. 2(b) and 2(e) show the results of Wendland's C 2 (W C 2 ) function without and with the greedy algorithm, respectively. The grids deformed by the continuous thin plate spline (CTPS) C 1 function are displayed in Figs. 2(c) and 2(f), where the greedy algorithm is applied in the second case. All RBFs used here fully recover the translation. Concerning the resulting skewness f skew , Fig. 2 shows that the W C 2 function leads to slightly worse values than for the two other RBFs. That is also visible in Table II based on the global measure of the mesh quality f skew . This observation is in contradiction to De Boer et al. [START_REF] Boer | Mesh deformation based on radial basis function interpolation[END_REF] who identified W C 2 as the best choice for RBFs with local support. For the hybrid greedy algorithm the value of m is set to the number of boundary points and the tolerance limit ε g defining the maximum absolute error to 10 -2 . The hybrid greedy method works well for all three RBFs with local support. A strong reduction of the floating-point operations and of the CPU-time is observed in Table II. Furthermore, based on Fig. 2 and the global quality measure f skew listed in Table II it is also evident that the quality of the mesh deformation only marginally differs between the cases with or without greedy selection. To give an idea about the efficiency of this algorithm, the number of points selected by the greedy procedure is listed in parentheses in Table II. This has to be compared to the entire number of boundary points N s = 1200 for the fine grid. Among the three local RBFs the hybrid node selection algorithm works best with W C 2 .

Figure 3 represents the grid deformed using the global TPS basis function providing a comparison of the coarse (21 × 21) and the fine grid (301 × 301). Note that in Fig. 3(b) only every 15 th grid point in both directions is shown. The quality of mesh deformations based on RBFs is theoretically independent of the number of nodes and their clustering. Thus, the deformed meshes are nearly identical. The same is true for the other RBFs not shown here. Figs. 23 Another interesting issue is the computational effort required by different RBFs. W C 2 takes much more time, whereas W C 0 is computationally the fastest. It is also seen that the effort required by CTPS C 1 and TPS is comparable. The matrices associated with different RBFs have a higher condition number along with more clustering of the eigenvalues near zero for the matrix corresponding to W C 2 . This explains its higher computational effort using the MINRES solver. Replacing the iterative solver by a direct LU decomposition the CPU-time required for the solution of the system does not depend on the chosen RBF anymore. The relative computational costs (compare values in Table II) are 12.96, 13.08, 21.07 and 16.98, for W C 0 , W C 2 , CTPS C 1 and TPS, respectively. Thus, for the present two-dimensional problem the direct solver leads to a drastic gain in CPU-time for the same grid quality.

3.1.4. IDW mesh deformation Finally, the IDW method is applied. Three different cases are considered to understand various effects of IDW. The first case depicted in Fig. 4(a) is obtained by directly interpolating the displacements with the weighting function of Luke et al. [START_REF] Luke | A fast mesh deformation method using explicit interpolation[END_REF]. The two other cases (Figs. 4(b) and 4(c)) are based on the interpolation of the translations and rotations (represented by quaternions) with the same weighting function. Here the parameter α is kept fixed and set to 0.1 which is especially suited for producing grid deformations of high quality according to Luke et al. [START_REF] Landry | Robust mesh-mover algorithms for hybrid meshes[END_REF]. Similar results are obtained with α dynamically predicted based on Eq. ( 18). Obviously, IDW is able to recover the translation without any difficulty. Although the clustered nature of the original grid is apparently preserved in Fig. 4(a), it has to be noticed that the orthogonality near the boundaries is not conserved. Computation leading to Figs. 4(b) uses 25 transformation division steps. Figure 4(c) is obtained using even 500 transformation division steps in order to check the stability of the transformation division algorithm. It is found that the algorithm performs in a sound manner with minor improvements of the mesh rendering the algorithm reliable. Both grids make an effort to preserve the orthogonality near all boundaries. That is an advantage S. SEN ET AL.

in comparison with the case without rotation. Although at some interior points, for example at the lower left side, the orthogonality is not preserved, such a situation is inevitable due to the singular character of the deformation. Thus, the fundamental property of the IDW procedure with rotation is the preservation of the orthogonality and the clustered nature of the grid near boundaries even under strained conditions.

As expected, a substantial increase of the computational effort is noted for the application of a large number of transformation division steps in Table II. It is found that the number of floatingpoint operations for IDW is directly proportional to the number of transformation division steps employed.

3.1.5. Summary All methods applied recover the translation perfectly. A comparison of the global skewness f skew obtained by TFI, RBF and IDW reveals that all four RBFs deform the mesh slightly more appropriate than TFI. IDW possesses a good orthogonality preserving property near boundaries, better than TFI and RBF. TFI is faster than IDW, which is faster than pure RBFs. However, it is worthwhile to point out that by usage of the greedy algorithm the interpolations using compact RBFs become computationally much more efficient and reach a similar level as IDW. For the present small problem size a direct solver speeds up the RBF method significantly, but practical applications necessitate an iterative solver. The IDW computation based on rotation requires many steps of the transformation division technique, which decelerates the mesh deformation significantly. However, as discussed in Section 2.3 and shown in Appendix B, this drawback is eliminated for time-resolved predictions relevant for the present area of interest, i.e., LES of turbulent flows.

A remark concerning the two measures in Table II to evaluate the computational effort has to be added. The number of floating-point operations has the advantage that it is independent of the machine architecture. A clear disadvantage is that it is fully unclear what the compiler makes out of the code and whether some operations may be vectorized or running in parallel. The relative CPU-time, however, comprises the floating-point operations but also other important issues such as the memory access leading to the overall performance of the algorithm. For example, comparing TFI and W C 0 yields a ratio of floating-point operations of about 720. The relative computational costs, however, increase only by a factor of 20, since the RBF algorithm runs more efficiently with a significantly increased number of floating-point operations per second. Therefore for practical purposes, it represents a more reliable measure of the computational effort and thus is solely used below.

Problem 2:

Rotation and deformation of an airfoil 3.2.1. Description of the problem To evaluate the three elementary mesh deformation strategies in situations involving large rotation of a body with a sharp angle, a 90 • rotation of an elliptic airfoil around the center of the profile is considered. The major and minor axes of the profile have a length of 1 and 0.1, respectively. When this airfoil rotates about the origin of the coordinate system in a unit step, it deforms into a NACA 0012 airfoil of the same chord length c. The initial undeformed O-grid and a reference final grid are both generated using conformal transformations [START_REF] Sen | A robust implicit compact scheme for two-dimensional unsteady flows with a biharmonic stream function formulation[END_REF]. The grid points at the outer radius (r = 12 c) remain fixed during the deformation process. The grid consists of 401 points equally distributed around the airfoil and 241 points in the wall-normal direction clustered towards the body leading to 802 boundary points.

TFI mesh deformation

From a theoretical point of view this type of rotational grid deformation is a major challenge for the TFI interpolation method. Although the grid depicted in Fig. 5(a) looks promising, a closer view in Fig. 5(b) reveals the failure of TFI regarding almost all aspects of this specific test case. Although TFI possesses the extreme preserving property, the grid deformation yields degenerated cells with cross-overs at the two extremities of the airfoil. The orthogonality is destroyed also at boundary points, where no degenerated cells appear. Furthermore, an artificial clustering of grid points on the surface of the airfoil appears, which is unacceptable for any numerical simulation. Hence as expected, TFI is not suitable in case of a large rotation of a structure with sharp edges.

RBF mesh deformation

The computations are carried out using six different RBFs. The five compact RBFs have a radius of support R = 15 in order to cover a major part of the domain under investigation. This value offers the best compromise between grid quality and computational effort. For four of these six RBFs further computations using the hybrid greedy algorithm are performed.

Table III provides an overview of all cases considered including the computational effort required and the global measure of the mesh quality f skew . The costs are provided as relative values to the corresponding time required for the mesh deformation using IDW ‡ . Obviously, Wendland's C 0 function fails to deform the mesh successfully leading to degenerated cells with cross-overs. A variation of the support radius R between 10 and 20 does not solve this problem. In a further test documented in Table III the convergence criterion for MINRES is decreased from ε c = 10 -5 to 10 -10 . Obviously, that has no effect and the deformation algorithm continues to produce degenerated cells. To the best of the authors' knowledge such a failure of W C 0 has not been reported in the literature. As visible in Fig. 6 the degenerated cells mainly appear in the vicinity of the trailing edge. This may be attributed to the combination of the rotation and the sudden change of the shape of the trailing edge. Table III. Problem 2: Comparison of relative computational costs (IDW as reference), global measure of the mesh quality f skew and remarks concerning the grid quality using RBFs (support radius set to 15) and IDW (Degenerated cells correspond to cross-overs of the cells and bad first cell to a highly and unacceptably deformed first cell on the airfoil). W C 2 is in general able to deform the mesh in a correct manner. The deformed mesh possesses the best quality in the vicinity of the airfoil in comparison with the other RBFs depicted in Fig. 7. It also closely resembles the grid around the NACA 0012 generated by the conformal transformation. However, in regions away from the airfoil the function W C 2 does not reach the same level of grid quality as the other RBFs. As clearly visible in Figs. 7, 8 and 9, the function W C 2 produces more skewed cells compared with the other basis functions investigated. That explains the lowest values of f skew found for W C 2 .

CTPS C 1 , CTPS C 2 a , CTPS C 2 b and TPS perform almost in a similar manner and the outcome perceptibly differs from the deformations prescribed by W C 2 . For the latter the effect of rotation on the grid has been successfully kept away from the vicinity of the airfoil, whereas in case of the various types of thin plate spline functions the effect of rotation is somehow distributed starting at the wall-nearest grid points and resulting in a skewed grid in the vicinity of the airfoil. It is noteworthy that all basis functions are able to cope with the near-singularity situation close to the trailing edge. In the regions away from the airfoil the CTPS C 1 and TPS basis functions are found to perform best as visible based on the skew quality metric contours in Figs. 7 to 9 and the values of the global skewness f skew in Table III. To summarize, W C 0 fails to appropriately deform the grid in case of a sharp-edged rotating structure. W C 2 leads to a quasi rigid-body movement of the grid points near the rotating airfoil, whereas the major grid deformation is shifted in the direction of the fixed boundaries. Thus, the quality of the deformed grid suffers in this region leading to more skewed cells. The four RBFs (CTPS C 1 , CTPS C 2 a , CTPS C 2 b and TPS) offer the best grid quality. The computational costs of CTPS C 1 and TPS are much lower than for CTPS C 2 a and CTPS C 2 b . W C 2 is the most expensive one. This trend closely resembles the relative costs reported for problem 1 in Table II and is explained by the ill-conditioned W C 2 matrix. Note that the computational effort with compact RBFs obviously depends on the support radius. Indeed, for RBFs without any data reduction technique the computational effort can be reduced by decreasing the support radius. Since this also deteriorates the grid quality, a compromise on R has to be found.

Next, it is of interest to investigate the performance of the hybrid greedy algorithm concerning data reduction and computational effort. If ε g is set to a high value (here ε g = 10 -2 ), the greedy algorithm generates an useless grid with degenerated cells for all RBFs tested. Applying a lower value of ε g = 10 -3 , CTPS C 1 and CTPS C 2 a yield an appropriate mesh. However, for W C 2 and CTPS C 2 b only a carefully chosen value of ε c leads to appropriate deformations without degenerated cells. In most cases of failure the deformed mesh is degenerated at the trailing edge. From a geometrical point of view it is mandatory for an efficient data reduction strategy to choose greedy boundary points in the vicinity of this trailing edge. In all computations the initial guess needed for the greedy algorithm is on purpose taken at a node of the leading edge. This particularly harsh choice is made in order to see the response of the greedy algorithm to situations involving FSI-LES predictions, where the appearance of a sudden singularity can not be ruled out. Nevertheless, if the mesh is not degenerated, it is almost indistinguishable from the deformed mesh obtained without the data reduction technique. Therefore, the corresponding counterparts to Figs. 8 and9 using the RBFs CTPS C 2 a and CTPS C 2 b are omitted here for the sake of brevity. As can be seen in Table III the hybrid greedy algorithm is very time consuming. The procedure iterates many times and is more expensive than the algorithm without any data reduction in most cases. Therefore, the use of a greedy procedure often makes sense only for problems with several successive mesh deformations of a nearly identical grid such as encountered in time-dependent FSI problems. In this case the greedy algorithm is applied for a certain time step. The chosen greedy boundary points are stored and used for many successive deformation cycles without renewal of the chosen boundary points, so that the CPU-time is globally reduced. However, if the grid is strongly deformed, it might be necessary to renew the choice of the boundary points after a while.

IDW mesh deformation

Finally, problem 2 is evaluated using the IDW mesh deformation technique with the value of α calculated by Eq. [START_REF] Boer | Mesh deformation based on radial basis function interpolation[END_REF]. The design of this interpolation strategy by Luke et al. [START_REF] Luke | A fast mesh deformation method using explicit interpolation[END_REF] is specifically targeted towards a method, which can interpolate the rotation with high accuracy. Therefore, IDW theoretically should be able to provide a high-quality grid for this problem. Since the origin is located at the center of the airfoil and the rotation is carried out with respect to this origin, no transformation division technique is employed.

The results for the same initial grid as used before are depicted in Fig. 10. Note that a nearly identical distribution of the skew quality metric is found on a much coarser grid evidencing the grid independence of the algorithm. Based on the work of Luke et al. [START_REF] Luke | A fast mesh deformation method using explicit interpolation[END_REF] the computations are repeated using a fixed value of α = 0.1. Again, no significant variation of the grid properties can be noticed. Furthermore, the computational effort for a fixed or a predicted value of α are almost identical (see Table III). In comparison to the different RBFs without greedy algorithm, IDW is at least a factor of 3 faster even in the case of the worst tolerance ε c = 10 -5 . For stricter tolerance values ε c the CPU-time ratio turns out to be even more clearly in favor of IDW. Here it is worthwhile to mention that the computational complexity of both RBF and IDW scales with N v × N s as noted earlier and hence for sufficiently large systems the CPU-times for both RBF and IDW may be isomorphic. The mesh deformation using IDW is nearly perfect, at least in the vicinity of the airfoil. A comparison with the grid generated by the conformal transformation reveals the effectiveness and significance of IDW for grid deformations involving significant rotation. The deformed grid is nearly orthogonal at all points close to the boundary of the airfoil, including the trailing edge. In regions away from the airfoil the skewness level increases. Thus, it is obvious that the IDW grid deformation is well-suited for pushing skewed grid cells away from the rotating body even with sharp edges. This is an important advantage for FSI-LES computations.

Summary

Based on the performance of the three methods (TFI, RBF and IDW) it can be concluded that RBF and IDW perform significantly better than TFI. IDW delivers a nearly perfect mesh in the vicinity of the airfoil, whereas all RBFs perform better at a certain distance from the body. Among all RBFs the performance of W C 2 is the best as far as the orthogonality preservation in the vicinity of the body is concerned. Thus, for this test case the outcome is in close agreement with the observations of De Boer et al. [START_REF] Boer | Mesh deformation based on radial basis function interpolation[END_REF]. It is worthwhile to mention that the IDW technique without the transformation division procedure takes much less time than any RBF. However, in case of very large systems the computational costs for both may be comparable. Globally, the greedy algorithm does not perform very well for all RBFs applied to problem 2 and shows its limits. This is due to the presence of a sharp edge on the moving body combined with the severe rotation. The use of a limited number of greedy boundary nodes implies that the imposed deformation is only correctly represented at these selected boundary points. In order to get an appropriate grid deformation, a sufficient number of boundary points describing the sharp geometry has to be selected. In the original greedy algorithm that is not assured. Hence, a possible improvement of the greedy algorithm might be to take the first cell height in the error analysis into account, so that more boundary points are taken in highly resolved near-wall regions. Moreover, the greedy algorithm is more expensive than the mesh deformation without any data reduction technique and thus only makes sense if the choice of the boundary points can be used several times.

HYBRID MESH DEFORMATION METHODS

In order to find a reasonable compromise between the quality of the deformed grid and the computational effort required for the deformation, block-structured grids offer the opportunity to S. SEN ET AL. combine the favorable properties of different methods while avoiding their drawbacks as far as possible. The hybrid mesh deformation methodology relies on this block structure and consists of a combination of two elementary mesh deformation algorithms: One for the deformation of the block interfaces and one for the inner points. In the following, the applied hybrid methods are denoted at first by the name of the technique applied for the block interface and in the second place by the name of the method used for the inner points. Since the number of grid points on the block interfaces only scales with N 2/3 v , major savings can be achieved if expensive techniques such as RBF or IDW are solely applied to the block interfaces and not to the entire grid. These hybrid grid deformation strategies are investigated based on the experience gained by the evaluation of the three elementary strategies. Furthermore, they are designed in such a manner that important advantages, i.e., low computational effort and high skew quality metric preservation are combined.

Since in FSI-LES predictions the CPU-time consumption is of paramount importance, it is anticipated that grid deformations based on RBFs applied to the entire domain are not feasible. However, taking the favorable mesh quality preservation property of RBF into account, it can be concluded that RBF may be of utmost importance when applied in a sub-domain. In such a situation a much smaller system of equations has to be inverted and therefore a huge amount of CPU-time can be saved. Furthermore, it has been noticed that for convex domains without substantial rotation TFI performs well and is very fast. The strength of IDW is its capacity to accurately capture the rotation and to preserve the orthogonality in the vicinity of the body, whereas it is found to be fall behind the others for a clustered grid. Based on these observations the following two hybrid methodologies are proposed:

• IDW-RBF: IDW is applied at the block interfaces to preserve more or less the orthogonality.

For the inner points of each block it is exploited that RBF achieves a high grid quality. Furthermore, since the sub-domains are processed in parallel by different processors, the system of equations related to the RBF implementation of the inner points may not turn out to be very costly. • IDW-TFI: Again IDW is applied at the block interfaces to retain the advantages of IDW, whereas TFI is used for the displacements of the inner points. Once the orthogonal character of the mesh is maintained at the interfaces, it is exploited that TFI produces high-quality grid deformations of the inner nodes very fast.

Both hybrid variants bypass the disadvantages of IDW for clustered meshes (see problem 1) to a large extent, since it is only used at the block interfaces and not for the nodes in the vicinity of the interface.

Additionally, a hybrid RBF-TFI methodology is used for comparison purposes. Since the solution of a huge system of equations involving all boundary nodes is required, this method is certainly not very interesting due to the required CPU-time.

For the sake of brevity and based on the experiences on the previous test cases considered, the application of the RBF mesh deformation technique is restricted to the TPS basis function with global support not taking the greedy algorithm into account. The fact that TPS leads to reasonable mesh deformations for problem 1 and 2 which are among the best motivates this choice. The predictions based on IDW are carried out using quaternions but without transformation division steps. These are not required here since the origin of the coordinate system lies at the center of rotation and all boundary points of the structure undergo the same rotation.

Description of the problem

The proposed hybrid algorithms and their elementary counterparts are tested in a situation with substantial rotation and deformation. Furthermore, the analysis is extended to the three-dimensional case and for deformations, which resemble the most important features of a standard test case for FSI-LES simulations [START_REF] Breuer | Fluid-structure interaction using a partitioned semiimplicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF][START_REF] De Nayer | Flow past a cylinder with a flexible splitter plate: a complementary experimental-numerical investigation and a new FSI test case (FSI-PfS-1a)[END_REF][START_REF] De Nayer | Numerical FSI investigation based on LES: Flow past a cylinder with a flexible splitter plate involving large deformations (FSI-PfS-2a)[END_REF]. In analogy to most FSI-LES cases the integration domain is divided into several sub-domains allowing to parallelize the computation based on domain decomposition with explicit message passing.

The undeformed block-structured grid along with the partitioning into ten blocks is depicted in Fig. 11. The geometrical setup is as follows: A square cylinder with a unit edge is placed in a plane rectangular channel with a height of 13 units, a length of 12 units and a width of 0.5 units. A splitter plate with a length of 5 units and a thickness of 0.25 units is attached to the square cylinder. The center of the cylinder coincides with the origin of the coordinate system. Note that the associated flow problem is a three-dimensional generalization of the nominally two-dimensional geometric configuration. In this study the cylinder is exposed to a rotation of 30 • around the center. Furthermore, the splitter plate is deformed by a trigonometric function according to 0.2 sin (2πx), where x is the distance from the point of attachment. Hence, the splitter plate is subjected to two different kinds of deformations. The grid possesses 8 cells in the spanwise direction. Each of the spanwise planes is further sub-divided into square cells of dimension 1/16 units leading to 363,033 nodes (= 193 × 209 × 9 without subtracting the nodes inside the structure) in the entire domain under consideration.

Results of the elementary and hybrid methods

The results obtained by the elementary and hybrid methods are depicted in Fig. 12. For the grid deformed by TFI-TFI a C 1 -discontinuity is observed at the interfaces between the blocks (see Fig. 12(a)). This grid discontinuity is not taken into account by the skew quality metric. A jump in the f skew level appears but its value remains high. Therefore, the harmonic average f skew reported in Table IV is found to be high. That shows that f skew is not sufficient for the evaluation of the overall grid quality. This sudden change at the block interface is attributed to the differences in the interpolating polynomials used in distinct sub-domains. Thus as noted for problem 1, the pure TFI algorithm is unable to maintain the orthogonality and is inappropriate for FSI-LES predictions. As visible in Fig. 12(b) the RBF-RBF deformation strategy works quite well. The orthogonality is preserved in the deformed grid. Although at the far field, especially in the downstream direction, some nodes are skewed, their skew quality metric is still high and a reasonable deformation of the grid leading to an overall high-quality grid can be observed. As usual IDW generates a mesh, on a real-life FSI problem (see Fig. 13): A flexible membranous hemisphere (diameter D = 1.5 × 10 -1 m) is attached on a flat plate and placed in an approaching thick turbulent boundary layer (Re = ρ air D U ∞ / µ air = 50,000) (see Fig. 14). This FSI-LES case is derived from the turbulent flow past a rigid hemisphere [START_REF] Wood | Experimental investigation and large-eddy simulation of the turbulent flow past a smooth and rigid hemisphere[END_REF] with two modifications: First, the grid is coarsened and consists of about 4.5 × 10 6 control volumes. Second, the deformation of the hemisphere is imposed based on the following movement of the surface in spherical coordinates (r FSI , θ, φ):

r FSI = A def z D/2
sin(ω t t) cos(4 φ) with A def = 1.1 × 10 -2 m and ω t = 2π 10 3 rad s -1 . Hence, it presently describes only a one-way coupled simulation which, however, is sufficient here for the evaluation of the hybrid technique. Note that this imposed movement is not physical. However, the chosen parameters leads to severe deformations in critical parts of the geometry: On the top of the hemisphere, on the sides and at the junction between the membrane and the bottom wall. The amplitude of the deformation A def is set to a value much larger than what is expected in a real physical case in order to stress the robustness of the proposed hybrid method. A snapshot of the geometry and the flow is shown in Fig. 15 at the time instant of the maximum deformation (t = 2.5 × 10 -4 s). The solver (FASTEST-3D) and the techniques used to compute this FSI-LES problem are described in Breuer et al. [START_REF] Breuer | Fluid-structure interaction using a partitioned semiimplicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF] and were validated by complementary experimentalnumerical FSI investigations [START_REF] De Nayer | Flow past a cylinder with a flexible splitter plate: a complementary experimental-numerical investigation and a new FSI test case (FSI-PfS-1a)[END_REF][START_REF] De Nayer | Numerical FSI investigation based on LES: Flow past a cylinder with a flexible splitter plate involving large deformations (FSI-PfS-2a)[END_REF].

Figures 15(a) and 15(b) show the vorticity magnitude of the turbulent flow around the rigid and the flexible hemisphere in two slices (y/D = 0, z/D = 0.014). Non-physical numerical oscillations are not detected, which is a first rough indicator for a grid of high quality. To evaluate the grid quality, the skew quality metric is plotted in the slice y/D = 0 in Figs. 15(c) and 15(d). The undeformed grid possesses a high quality. The first cells are all perfectly orthogonal to the wall and the skew quality metric is close to unity in the whole mesh. After the deformation the quality of the grid is slightly altered in several regions around the hemisphere (see Fig. 15(d)). However, the skew quality metric remains high in these areas (> 0.9). Furthermore, it is important to mention that the high grid quality is maintained near the moving membrane (see Fig. 15(e)) and near the bottom wall (see Fig. 15(f)). As explained in Section 3.2.4 this is one of the main advantages of a method based on IDW: The deformation is pushed away from the FSI interface and the node distribution remains unchanged in the wall vicinity. This particularity ensures that the boundary layer can be correctly resolved by the LES solver.

The junction between the bottom wall and the flexible hemisphere is a critical region during the mesh deformation based on IDW. In case of a wall-resolved LES the mesh is very fine here. When the structure inflates, the area is reduced and degenerated cells can appear. In case of the hybrid IDW-TFI method, IDW is only applied on the block boundaries, which are far away from the FSI interface. Therefore, this is not an issue for the present flexible hemisphere. The mesh does not contain any bad cells and the orthogonality is more or less conserved also in this critical region. Concerning the CPU-time consumption the hybrid IDW-TFI method is about 6.7 times faster than the pure IDW method for exactly the same grid quality (see Table V). Note that the corresponding ratio of floating-point operations is about 7.3. Thus the pure IDW methods needs more floatingpoint operations but runs slightly more efficiently than the hybrid IDW-TFI method. No comparison with the elementary TFI and RBF methods is possible. The former leads to cross-overs, whereas for the latter the system of equations can not be solved in a reasonable time due to the large number of boundary points.

CONCLUSIONS

Three elementary mesh deformation methods are thoroughly investigated regarding their responses to situations involving substantial deformations. As expected, TFI is able to deform the grid quite fast, but is incapable of dealing with substantial rotation. Moreover, in case of a block-structured grid the deformed mesh is not smooth at the block interfaces. Nevertheless, in convex domains with prescribed deformations on the boundaries, its performance is reliable, although some loss of orthogonality is observed. RBF on the other hand produces a mesh of high quality in almost all situations, but is computationally expensive. Direct solvers for RBF are faster than iterative ones but limited to small problems. Thus, the latter is preferred for general utilization. At a glance of the test cases studied, the TPS basis function with global support and Wendland's C 2 function with local support are the best choices. W C 2 is able to preserve the orthogonality, but is timeconsuming. TPS represents a good compromise between quality and speed. A well-thought-out use of IDW is also an attractive alternative guaranteeing a reasonable computational effort. An IDW method split into a rotational and translational deformation maintains a perfect orthogonality in the near-wall region, which is mandatory for LES. However, it requires to set many problem-dependent parameters, where unfavorable choices can lead to a degenerated mesh with cross-overs. The IDW grid deformation relying on a rotation and a translation step introduces a coordinate dependence. To solve this problem the time-consuming transformation division steps can be applied. However, for FSI-LES simulations involving small time steps and thus small relative displacements between subsequent time steps, IDW can renounce this additional procedure. For this purpose is has to be implemented in such a manner that it relies on the last deformed grid mimicking the transformation division steps and thus eliminating the coordinate dependence and the high computational costs.

Hybrid methods represent an appropriate compromise between quality requirements and computational effort. Since for block-structured grids the deformation consists of two steps (deformation of the block interfaces [∼ N 2/3 v ] and displacement of the inner points [∼ N v ]), in principle six combinations of TFI, RBF and IDW are possible. Since the application of TFI for the block interfaces in not appropriate and the hybrid method RBF-IDW is canceled owing to the expected high effort, three possible hybrid variants remain: IDW-TFI, IDW-RBF and RBF-TFI. These are analyzed based on a realistic FSI scenario.

It is found that the deformation of the block boundaries plays the major role for the quality of block-structured grids. The block boundaries of a typical FSI-LES prediction are best adapted by IDW since it correctly captures the rotation of a deforming body. Once these block boundaries are deformed in an efficient manner, the deformation of the inner mesh is less critical and can be performed either by TFI or RBF. In case of substantial deformations RBF is favored for the inner deformation, but even in conjunction with the greedy algorithm the computational effort is still extremely high, which renders it unfeasible for a regular usage within each time step or even each FSI sub-iteration step of a coupled FSI-LES prediction. In many cases the application of TFI within the second step of the hybrid algorithm leads to deformed grids of high quality, while the CPU-time requirements stay at a reasonable level. Based on these arguments the hybrid IDW-TFI methodology is clearly favored and additionally tested based on a real-life FSI problem of a hemispherical membranous dome with an imposed time-dependent movement. It turns out that the proposed hybrid IDW-TFI method appropriately deforms the grid maintaining the original distribution and the orthogonality of the cells in the near-wall region. Thus, it generates high-quality grids in a reasonable time as mandatory for FSI-LES.

S. SEN ET AL.

41 × 22 and a finer one of size 401 × 202. Since the beam is immaterial, its thickness has been set to one grid spacing in y-direction. The beam is deformed by prescribing a parabolic displacement y = x 2 , where x is the distance from the center of the beam.

In IDW all the displacements are determined in terms of a global coordinate system, but the quaternions are calculated locally [START_REF] Maruyama | High-quality mesh deformation using quaternions for orthogonality preservation[END_REF]. This leads to a situation where the translation vector, computed at a surface point, depends on its global position, whereas the rotation, expressed by quaternions, does not. Hence, the displacements of the interpolated volume grid points loses symmetry.

The so-called transformation division procedure introduced by Maruyama et al. [START_REF] Maruyama | High-quality mesh deformation using quaternions for orthogonality preservation[END_REF] intends to make this interpolation independent of the coordinate system by dividing the deformation into small steps. Thus, it is a finite implementation of an essentially infinite procedure. In case of an unsteady computation it is found that with a small time step the deformations are most often small in each step. Hence the transformation division technique may not be necessary. To validate this procedure, the problem is tackled as a steady and an unsteady process, where the latter is completed in unit non-dimensional time.

In Fig. 17 the deformed mesh and its skew quality metric predicted by various strategies are presented. Figures 17(a), 17(b) and 17(c) are obtained without using the transformation division technique. In Fig. 17(a) a symmetric deformation is visible because the origin of the coordinate system is located at the center of the geometry. However, a non-symmetric grid is predicted if this origin is shifted away from the center (see Fig. 17(b)). In Fig. 17(c) a time-marching strategy is used dividing the entire deformation into 25 time steps. Obviously, this leads to a symmetric mesh. The transformation division procedure is applied in the fourth case shown in Fig. 17 a). Both procedures relying on a splitting of the entire deformation into a finite number of sub-steps yield a higher grid quality than the original procedure. Thus, it can be concluded that making the node displacement procedure independent of the coordinate system is fundamental.

3. 1 .

 1 Problem 1: Deformation of a two-dimensional square region 3.1.1. Description of the problem As depicted in Fig.

  Deformed mesh and f skew contour.

Figure 1 .

 1 Figure 1. Problem 1: Undeformed mesh and TFI deformation.

Figure 2 .

 2 Figure 2. Problem 1: RBF mesh deformation with compact support (21 × 21 grid) and f skew contour.

Figure 3 .

 3 Figure 3. Problem 1: RBF mesh deformation with global support (TPS) and f skew contour. (For the fine grid only every 15 th grid point in both directions is shown).

Figure 4 .

 4 Figure 4. Problem 1: IDW mesh deformations (α = 0.1) and f skew contours.

  (a) Full view. (b) Close view.

Figure 5 .

 5 Figure 5. Problem 2: TFI grid deformation (every 5 th grid point in each direction is displayed).

  (a) View of trailing edge. (b) View of leading edge.

Figure 6 .

 6 Figure 6. Problem 2: Mesh deformation using W C 0 .

  (a) W C 2 full view. (b) CTPS C 1 full view. (c) TPS full view. (d) W C 2 close view. (e) CTPS C 1 close view. (f) TPS close view.

Figure 7 . 2 :

 72 Figure 7. Problem 2: Mesh deformation using RBFs (only every 5 th grid point along each direction shown).

  S. SEN ET AL. (a) Full view. (b) Close view.

Figure 8 .

 8 Figure 8. Problem 2: Mesh deformation using CTPS C 2 a without and with the hybrid greedy algorithm (Since the results in both cases are nearly identical, only one set is shown here. Note that only every 5 th grid point along each direction is depicted).

Figure 9 .

 9 Figure 9. Problem 2: Mesh deformation using CTPS C 2 b without and with the hybrid greedy algorithm (only every 5 th grid point along each direction shown).

  (a) Full view. (b) Close view.

Figure 10 .

 10 Figure 10. Problem 2: Mesh deformation using IDW and α calculated by Eq. (18) (only every 5 th grid point along each direction shown).

Figure 11 .

 11 Figure 11. Sketch of the topology of the grid structure with the splitting into ten blocks.

Figure 13 .Figure 14 .

 1314 Figure 13. Air inflated structure in Gentofte (Denmark).

  S. SEN ET AL. (a) Full view of the rigid hemisphere: Slices (y/D = 0, z/D = 0.014) with the vorticity magnitude. (b) Full view of the deformed hemisphere: Slices (y/D = 0, z/D = 0.014) with the vorticity magnitude. (c) Full view of the rigid hemisphere: Slice (y/D = 0) with the grid skew quality metric. (d) Full view of the deformed hemisphere: Slice (y/D = 0) with the grid skew quality metric. (e) Close view of the deformed grid at the top: Slice (y/D = 0) with the grid skew quality metric. (f) Close view of the deformed grid at the bottom: Slice (y/D = 0) with the grid skew quality metric.

Figure 15 .

 15 Figure 15. Snapshot of the flexible hemisphere in a turbulent flow predicted by LES.

  (d). Again 25 steps similar to the time-marching procedure are used. Obviously the results obtained by the timemarching strategy and the transformation division algorithm are indistinguishable from each other. A particularly interesting observation that requires attention is that neither Fig.

17

 17

  (c) nor Fig.17

  (d) matches with Fig.17

  (

  IDW with origin located at the center of the mesh.

  IDW with origin located away from the center of the mesh.

  IDW with origin located away from the center of the mesh but based on an unsteady computation. IDW with origin located away from the center of the mesh but using the transformation division procedure.

Figure 17 .

 17 Figure 17. Deformed mesh of the beam and skew quality metric computed by IDW with or without the transformation division procedure.

Table I .

 I Radial basis functions.

  and Table II reveal that all four RBFs are overall leading to a reasonable mesh deformation but none is able to maintain the orthogonality particularly at the challenging bottom side. The grid clustering is best preserved by W C 2 and TPS. W C 0 and CTPS C 1 tend to shrink the region of clustered points under severe deformation.

	1.2							1.2
	1						1	1						1
							0.9							0.9
	0.8						0.8 0.7	0.8						0.8 0.7
							0.6							0.6
	0.6						0.5	0.6						0.5
							0.4							0.4
	0.4						0.3 0.2	0.4						0.3 0.2
							0.1							0.1
	0.2						0	0.2						0
	0	1	1.2	1.4	1.6	1.8	2	0	1	1.2	1.4	1.6	1.8	2
			(a) 21 × 21 grid.			(b) 301 × 301 grid.

Table IV .

 IV Comparison of relative computational costs and global measure of the mesh quality f skew for elementary and hybrid mesh deformation strategies applied to a block-structured grid.

	Method TFI-TFI RBF-RBF IDW-IDW RBF-TFI IDW-RBF IDW-TFI	Rel. comp. costs Remark 1 Block interfaces strongly non-smooth 2013.71 Overall high grid quality 19.71 High grid quality near moving boundaries 2000.78 Overall high grid quality 418.00 Overall good grid quality 5.46 Best compromise of quality & CPU-time	f skew 0.96 0.96 0.96 0.96 0.95 0.95

Table V .

 V Comparison of the floating-point operations (GFLOP for one time step) and of the relative computational costs for elementary and hybrid mesh deformation strategies applied to a real-life FSI problem.

	Method TFI-TFI RBF-RBF IDW-IDW IDW-TFI	GFLOP Rel. comp. costs ----9675 6.68 1328 1	Remark Cross-overs System too large to be solved High grid quality near moving boundaries Best compromise of quality & CPU-time

International Journal for Numerical Methods in Engineering, vol. 111 (3), pp. 273-300, (2017).

International Journal for Numerical Methods in Engineering, vol. (3), pp. 273-300, (2017).

Copyright c 2017 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2017) Prepared using nmeauth.cls DOI: 10.1002/nme.5465 International Journal for Numerical Methods in Engineering, vol. 111 (3), pp. 273-300, (2017).

which preserves the orthogonality in the vicinity of the deformed splitter plate to a great extent. However, in the region downstream of the splitter plate some cells with a relatively high skewness are observed in Fig. 12(c). This is consistent with earlier observations that the IDW interpolation pushes the skewness away from the deforming body.

The two hybrid versions involving IDW for the block interfaces (see Fig. 12(e) and 12(f)) lead to similar grid qualities as obtained by the pure IDW method. Similarly the hybrid RBF-TFI method (see Fig. 12(d)) delivers a grid of almost the same quality as the pure RBF technique. Thus, the lesson which could be learned from the present results is that the method applied for the deformation of the block interfaces is of major importance for the grid quality achieved, whereas the method used for the displacements of the inner grid points plays a minor role. The corresponding relative computational efforts related to all six computations are given in Table IV. Similar to the previous cases huge differences are already found between the pure elementary methods. As expected, pure TFI produces deformed meshes with the lowest computational effort. TFI is about 2000 times faster than pure RBF and still about 20 times faster than pure IDW. As visible in Table IV the major disadvantage of any RBF implementation is the large computational effort required consumed. As anticipated a combination of RBF and TFI only slightly reduces the computational costs in comparison with the pure RBF algorithm, since the time-consuming determination of the coefficients of the basis functions are required in both variants. Hence, this hybrid technique is also not suitable. Combining IDW with RBF reduces the computational costs of the pure RBF method significantly. However, it is quite slow compared to the IDW-TFI method. The hybrid IDW-TFI method is very promising since it only takes about one quarter of the computational effort of the pure IDW technique while a grid of comparably high quality is achieved.

APPLICATION OF THE HYBRID METHOD TO A REAL-LIFE FSI PROBLEM

In the previous section the hybrid IDW-TFI method proves its ability to deform the mesh conserving a high quality with reasonable computational effort. This hybrid method is now evaluated based S. SEN ET AL.

ACKNOWLEDGMENTS

The first author gratefully acknowledges the Indian National Science Academy (INSA), India and the Deutsche Forschungsgemeinschaft (DFG), Germany for their patronage and financial support under the INSA-DFG bilateral collaboration program. Financial support received from DFG under contract number BR 1847/15-1 allowed the first author to stay at the Helmut-Schmidt Universität, Hamburg for three months. Furthermore, the project is financially supported by the DFG under the contract number BR 1847/12-2.

A. DECOMPOSITION OF THE DISPLACEMENTS INTO ROTATION (MODELED WITH QUATERNIONS) AND TRANSLATION

The displacement at each boundary node can be decomposed into a rotation (quaternion Q) and a translation T . Q and T can be determined by using the steps outlined by Maruyama et al. [START_REF] Maruyama | High-quality mesh deformation using quaternions for orthogonality preservation[END_REF]:

1. Determine unit normal vectors nu and nd at each surface grid point corresponding to the undeformed and deformed mesh, respectively. This can be accomplished in two steps. The first step is to represent the coordinates of the neighboring nodes with respect to the node under consideration (green node in Fig. 16). These vertices are now denoted R u = [0, r u ] or R d = [0, r d ] for the undeformed and deformed mesh, respectively. A schematic representation for one point can be found in Fig. 16(a) and 16(b), respectively. The second step is to predict the unit normals to the face corner position vectors and to take the average leading to the computation of nu and nd , respectively. 2. The first quaternion Q 1 is then defined as

The quaternion Q 1 represents the change of the surface orientation during deformation and aligns nu to nd as depicted in Fig. 16(c). 3. Now the position vectors of the undeformed face are updated to R u = [0, r u ] using the relation

4. Next, it is necessary to rotate the updated undeformed surface in order to align it as good as possible with the corners of the deformed face. The corresponding rotation axis is nd and the corresponding angle θ 2 can be calculated as follows:

where N vrtx is the number of adjacent vertices. This process can be seen in Fig. 16(d). Thus, a second quaternion Q 2 which represents the rotation due to the cell shape modification during the mesh deformation procedure can be defined as

5. Finally, the quaternion Q is attained by using the composition of Q 1 and Q 2 ,

Once the rotation quaternion Q is determined at a boundary node, the translation quaternion T = R d -QR u Q * can be calculated following Eq. [START_REF] Rendall | Unified fluid-structure interpolation and mesh motion using radial basis functions[END_REF]. However, this translation vector is coordinate dependent since the rotation has been calculated with respect to a local coordinate system. To

16. Representation of the procedure to obtain quaternions at a surface node having four neighboring nodes.

eliminate this coordinate dependence, a strategy called transformation division algorithm has been proposed by Maruyama et al. [32]. This algorithm consists of repeating the identical elementary process m times. Each iteration of this elementary transformation is composed of a rotation about the same axis but with an incremental angle θ m, where θ is the angle of rotation prescribed by the quaternion Q. This is followed by the corresponding calculation of the translation vector. The entire procedure can be summarized as

The corresponding calculation of the translation vector is based on the relation

where N is the 3 × 3 element matrix corresponding to the quaternion Q and

B. COORDINATE DEPENDENCY OF IDW AND TRANSFORMATION DIVISION STEPS

In order to illustrate the coordinate dependency of the interpolated quaternions used for IDW leading to different grid deformations for an identical structural movement, the simple test case of a bending beam in a rectangular region is considered. Furthermore, this setup allows to discuss the subsequent strategy developed to avoid this phenomenon. Assuming a two-dimensional problem, a beam of unit length is immersed in a rectangle of dimension 4 × 2. Computations are carried out using two different grids, a coarse mesh of size