# A GPU-Accelerated 3D Mesh Deformation Method Based on Radial Basis Function Interpolation

Authors: Jiandong He; Chong Wu; Yining Jia
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, we developed a GPU parallelized 3D mesh deformation based on the radial basis function (RBF) interpolation using NVIDIA CUDA C++. The RBF mesh deformation method interpolates displacements of the boundary nodes to the whole mesh, which can handle large mesh deformations caused by translations, rotations, and deformations. However, the computational performance of RBF mesh deformation depends on the quantity of grids. For 3D mesh deformation, especially for mesh with large number of boundary nodes, RBF interpolation has been verified computationally intensive. By introducing GPU acceleration, the computational performance of RBF mesh deformation code improved significantly. Several benchmark test cases were carried out to show the feasibility and efficiency of GPU parallelization. In summary, the GPU parallelized RBF interpolation shows the potential to become an alternative way to deal with 3D mesh deformation problems in an efficient way.


---

## Body

## 1. Introduction

In recent decades, a numerical simulation method such as Computational Fluid Dynamics (CFD) has become an efficient way to solve complex scientific and engineering problems. For simulations such as fluid-structure interaction and aerodynamic shape optimization which typically involve moving boundaries of fluid domain, mesh deformation is needed to adapting the computational grid to the new domain.Mesh deformation based on the radial basis function (RBF) interpolation method was first developed in 2007 [Boer et al. [1]]. Compared to traditional mesh deformation methods such as the spring analogy [Batina (1989); Farhat et al. [2]; Degand and Farhat [3]] and transfinite interpolation (TFI) method [Dubuc et al. [4]], RBF interpolation has many advantages: (1) As a point-to-point interpolation method, RBF interpolation requires no grid connective information and thus is much easier to implement; (2) RBF interpolation is more flexible and can be applied to both structured and unstructured grids; and (3) RBF interpolation can handle large deformations of a mesh while maintaining an acceptable geometrical accuracy and mesh quality [Boer et al. [1]]. In the past few decades, RBF interpolation has been applied in many numerical simulation procedures such as aerodynamic shape optimization [Allen and Rendall [5]; Morris et al. [6]], CFD simulation of flapping wings [Bos et al. [7]], and fluid-structure interaction problems [Rendall and Allen [8]]. RBF interpolation has been stated as more computationally efficient than other mesh deformation method in literature [Boer et al. [1]; Morris et al. [6]]; however, this conclusion is established in the case of mesh with small number of grid nodes. As the grid quantity grows, the computational cost increased drastically as presented in Section 4, especially for 3D mesh deformation with a large number of grids. Due to this reason, RBF mesh deformation has not been popularized in actual numerical simulations.Since the introduction of RBF interpolation mesh deformation method, a lot of efforts has been taken to further improve its efficiency. These works can be divided into two categories: algorithm and hardware perspective. Rendall and Allen developed an RBF interpolation with the data reduction algorithm [Rendall and Allen [9]]. In a study [Rendall and Allen [8]], a data reduction RBF mesh deformation method based on the greedy algorithm was proposed. As demonstrated in Section 2, a system of linear equations consists of nb equations which need to be solved in the process of RBF interpolation, where nb is the number of boundary nodes. For data reduction RBF mesh deformation, only a small subset of boundary nodes is selected as control points, so the computational cost reduced drastically. However, this method brought in interpolation errors and geometry precision declines [Gillebaart et al. [10]; Niu et al. [11]]. There are also some works to improve the computational efficiency of RBF mesh deformation by the parallel algorithm with message passing interface (MPI) and run RBF codes on multiple CPU cores [Chao et al. [12]; Fang et al. [13]].In recent decades, the computational performance of graphic processing units (GPU) has reached to a great extent [Owens et al. [14]], coupled with the rapid development of its computing platform and application programming interface such as CUDA [Shuai et al. [15]]. CUDA is a general purpose parallel computing platform introduced by Nvidia company in 2006, which employs GPUs to solve various computationally intensive problems more efficiently than multiple CPU cores [Nvidia [16]]. Nowadays, GPU computing has extended its application from computer vision to many fields like numerical simulation of fluid flow [Harris [17]; Rogers et al. [18]; Liu et al. [19]] and structure deformation [He et al. [20]; He and Lei [21] problems. However, to our knowledge, there are few works that have attempted to accelerate the mesh deformation process on the GPU platform.The motivation of this paper is to develop a GPU-accelerated 3D mesh deformation algorithm based on RBF interpolation. The paper is organized as follows: In Section2, we introduced RBF formulation used in this paper, including the schemes and equations. In Section 3, the details of GPU parallelization with CUDA programming are presented. Test cases are carried out in Section 4 to validate the reliability of used RBF interpolation and the efficiency of GPU parallelization. Conclusions and future works are given in Section 5.


## 2. RBF Interpolation Formulation

Radial basis functions (RBF) refer to a series of functions in the form:(1)Φx=φx,where x is Euclidean distance between two points, which refers to the spatial distance between two grid nodes of the mesh. There are various forms of radial basis functions which can be used to interpolating the data, and they can be divided into two groups: (1) functions with compact support and (2) functions with global support. Here, we chose the functions with compact support, and the distance x is normalized by support radius r: ξ=x/r. Then (1) can be transformed into the form as follows:(2)Φr=φξ=fξ,0≤ξ≤1,0,ξ>1.When a radial basis function with compact support is applied, the influence of control points will be constricted inside a circle (for 2D) or sphere (for 3D) of compact radiusr. The commonly used forms of fξ are listed in Table 1. In this paper, we choose the function with form of CP C4 [21].Table 1
Radial basis functions with compact support.
NoNamefξ1CPC01−ξ22CPC21−ξ44ξ+13CPC41−ξ635/3ξ2+6ξ+14CPC61−ξ832ξ3+25ξ2+8ξ+15CTPSC01−ξ56CTPSC11+80/3ξ2−40ξ3+15ξ4−8/3ξ5+20ξ2logξ7CTPSCa21−30ξ2−10ξ3+45ξ4−6ξ5−60ξ3logξ8CTPSCb21−20ξ2+80ξ3−45ξ4−16ξ5+60ξ4logξIn the process of mesh deformation based on RBF interpolation, the displacements of internal grid nodessx can be approximated by the sum of basis functions as (3)sx=∑j=1nbαjφx−xbj+px,where xbj=xbj,ybj,zbj are the coordinates of control points, which refers to boundary nodes. nb is the total number of boundary nodes. x−xbj is the distance between an internal grid point and j-th boundary node. The coefficients αj polynomial px are determined by solving the linear equations as (4)MPbPbT0αβ=db0,where M is a square matrix of dimension nb×nb with the form as(5)M=φb1b1φb1b2⋯φb1bnb⋮⋮⋱⋮φbnbb1φbnbb2⋯φbnbbnb,and φbibj=φxbi−xbj, which is radical basis funtion of distance between the i-th and j-th boundary node. Pb is a matrix of dimension nb×4 which has the form as follows:(6)Pb=1xb1yb1zb1⋮⋮⋮⋮1xbnbybnbzbnb.In vectorα, there are interpolation coefficients of boundary nodes presented as(7)α=αb1⋮αbnb,and vector β contains the coefficients of polynomial px. Vector db is the prescribed displacement of boundary nodes:(8)db=db1⋮dbnb.


## 3. CUDA Implementation of RBF Mesh Deformation

In Section2, we introduced the RBF formulation for mesh deformation. Here, we present the GPU parallelization of RBF interpolation using CUDA C++. The serial version CPU code is also developed as a reference. The procedure of RBF mesh deformation can be decomposed to following subroutines: (1) importing grid points from mesh file; (2) constructing the matrix M, as described in Section 2; (3) solving the linear system of (4); (4) calculating the desired moving vector of inner grid points; and (5). updating the coordinates of boundary and inner grid nodes.For a CUDA program, the CPU is regarded as host to perform the task schedule, and the GPU is regarded as a device to deal with the parallel computing part of the program. Usually a CUDA program starts from the host side and copies data from the host side to the device side, then the host calls the kernel function executed on the device. By analyzing the algorithm of RBF interpolation and the computational cost of subroutines measured from the serial CPU program, we found that most of the subroutines in RBF interpolation can be parallelized and executed on GPU. The flow chart of CUDA implementation of RBF deformation in this work is shown in Figure1.Figure 1
Flow chart of CUDA implementation of RBF mesh deformation.In this work, we adopted four GPU kernel functions to accelerate the RBF interpolation. The first one is for constructing the matrixM, the second one is for solving the linear equation, the third one is for computing the displacement of inner grid nodes, and the fourth one is for updating coordinates of the inner and boundary nodes. The pseudo code of the kernel functions are illustrated in Algorithms 1–4.Algorithm 1: Kernel function for constructing matrix M.
function MatrixConstructionCbInput: MatrixCb containing the coordinates of boundary grid nodesOutput: MatrixM described in Section 2/∗ thread indextid, number of inner grid nodes nb ∗/iftid<nb thenfori=1 to nb dox=xtid−xi; //x: distance between two boundary nodes/∗r: support radius ∗/ifx<r thenMtid,i=φx; //φ: Radial basis functionelseMtid,i=0endendendAlgorithm 2: Kernel function for solving the linear equation system of equation (4).
function LinearSolverA,b/∗ Solve the linear equation system with the form ofAα=b ∗/Input: MatrixA of dimension nb+4×nb+4, vector bOutput: vectorαcusolverDnDgetrfA; //LU decomposition of matrix AcusolverDnDgetrsL,U,b; //solving linear equation system: LUα=bAlgorithm 3: Kernel function for computing the displacement of inner grid nodes.
/∗ Matrixxi containing the coordinates of inner grid nodes ∗//∗ Matrixxb containing the coordinates of boundary grid nodes ∗//∗ vectorα and β containing the coefficients from solving the linear equation system ∗/function ComputeDispxi,xb,α,βInput: Matrixxi, xb and vector α,βOutput: Displacementssx described in Section 2/∗ thread indextid, index of inner nodes ni ∗/iftid<ni thenforb=1 to nb dox=xtid−xb; /‖x‖: distance between a inner grid node and a boundary grid node/∗r: support radius ∗/ifx<r thencalculateφx; //φ: Radial basis functionelsesetφx=0endsum=sum+αiφx+pxendcalculate the displacements of inner grid nodesstid according to (3)endAlgorithm 4: Kernel function for updating the coordinates of inner and boundary grid nodes.
function MatrixConstructionC,s,δInput: MatrixC containing the coordinates of inner and boundary grid nodes, vector s and δ containing the displacements of inner grid nodes and prescribed displacements of boundary nodesOutput: MatrixCu: updated coordinates of inner and nboundary nodes/∗ thread indextid, number of total grid nodes nt ∗/iftid<nb thenCutid=Ctid+δtid; //update coordinates of boundary grid nodesendifnb<=tid<nt thenCutid=Ctid+stid; //update coordinates of inner grid nodesend


## 4. Results and Discussion

In this section, we will study two test cases for 2D and 3D mesh deformation based on the RBF interpolation algorithm described in Section2. Firstly, mesh deformation around an 2D airfoil with relatively smaller number of grid nodes is carried out. Then, we presented the mesh deformation of a 3D wing with large number of grid nodes. The computational cost of serial CPU and GPU parallelized code are compared for these test cases with different quantities of inner grid nodes and boundary nodes to verify the capability of GPU acceleration. In this work, the serial C++ code runs on Intel core(M) i7-8700K while the CUDA C++ code runs on NVIDIA GeForce RTX 2080Ti graphic processing units.


### 4.1. Mesh Deformation of Flow Field around an Airfoil

In this subsection, we perform a series of test cases for mesh deformation of flow field around an NACA0012 airfoil with chord lengthc=1.0, and the boundary of far field is a circle of radius R=5.0, as shown in Figure 2. During the process of mesh deformation, grid nodes on the upper and lower chamber of the airfoil move according to (9)Δy=αx2,where Δy is the displacement along the y direction of boundary grid nodes on the airfoil, x is the x coordinate of these nodes, and α=0.8 is a prescribed coefficient. The deformation is carried out for 10 steps which leads to maximum displacement Δymax=8.0 at the trailing edge of the airfoil. There are four meshes of different resolutions which are investigated, as listed in Table 2. The initial and deformed mesh with number of grid nodes nt=1.0×104 and boundary nodes nb=200 are shown in Figure 2. As shown in the figure, the mesh is deformed smoothly without negative volume cells. The change of mesh quality after deformation is not examined in this paper since it has been studied in many literatures [Boer et al. [1]; [22]].Figure 2
Deformation of the mesh around an airfoil. (a) initial mesh. (b) Final mesh using RBF interpolation. (c) Close up of the initial mesh. (d) Close up of the final mesh using RBF interpolation.
(a)(b)(c)(d)Table 2
Efficiency comparison of the serial CPU code and GPU parallelized code for 2D RBF mesh deformation.
Mesh nonbntΔtCsΔtGsSpeed up1501.0×1030.01020.19000.053×21001.0×1030.00990.19500.051×32001.0×1040.27940.19701.51×44001.0×1040.28500.20451.39×54001.0×10518.74162.01729.29×68001.0×10519.06511.98349.61×In Table2, the execution time of GPU-accelerated code and the corresponding CPU code for different movement types and resolutions are presented, where nb is the number of boundary nodes, nt is the total number of grid nodes, and ΔtC and ΔtG are the measured execution time spent on CPU and GPU for single deformation step, repectively. As shown in Table 2, in the case of small quantity of grid nodes (around 1.0×103), the execution time of CUDA code is longer than the CPU code. We analyse the reason that GPUs cannot be fully exploited in this situation, and the execution time of memory copy shares a large proportion. When the total number of grid nodes increases to the magnitude of 1.0×104, the GPU code becomes faster than the CPU code. As the total number of grid nodes grows to the level of 1.0×105, more than nine times acceleration was attained for the GPU code. Also, we found that both the CPU and GPU code are not sensitive to the number of boundary nodes in these cases.In Table3, the efficiencies of CPU subroutines and corresponding GPU kernel functions for 2D mesh deformation (nb=800,nt=1.0×105) are compared, in which subroutine No. 1 to No. 4 corresponding to the algorithms are presented in Section 3.Table 3
Efficiency comparison of CPU subroutines and GPU kernel functions for 2D RBF mesh deformation (nb=800,nt=1.0×105).
Subroutine noCPU (s)GPU (s)Speed up11.40320.19827.080×23.85771.20353.205×313.58580.521926.031×40.21540.05983.602×


### 4.2. Deformation of the Mesh around a 3D Wing

In this subsection, we investigate test cases for mesh deformation of flow field around a 3D wing, the configuration of the wing is extruded from a NACA1410 airfoil with chord lengthc=1.0, and the boundary of far field is a cylinder of radius R=5.0 and length L=10.0. As presented in Table 3, a series of test cases with different movement types and resolutions are performed. Figure 3 shows the initial and final deformed meshes.Figure 3
Deformation of the mesh around a 3D wing. (a) Cross-section of the initial mesh. (b) 3D view of the initial mesh. (c) Final mesh with deflection. (d) Final mesh with translation and rotation.
(a)(b)(c)(d)In Table4, we compared the efficiencies of GPU-accelerated code and the corresponding CPU code for different movement types and resolutions, where nb is the number of boundary nodes, nt is the total number of grid nodes, and ΔtC and ΔtG are the measured execution time spent on CPU and GPU for the single deformation step. As shown in the table, for these 3D mesh deformation cases, tens to hundreds of times acceleration is achieved via GPU parallelization. We also found that for these 3D mesh deformation cases, as the number of grid nodes grows, the execution time of CPU codes becomes sensitive to the number of boundary nodes while the corresponding GPU codes do not.Table 4
Efficiency comparison of the serial CPU code and GPU parallelized code for 3D RBF mesh deformation.
Mesh nonbntΔtCsΔtGsSpeed up15005.0×1049.61041.78625.37×210005.0×10419.60971.769911.08×315002.5×105150.83806.663522.64×430002.5×105432.68666.573865.82×530005.0×105877.683912.328471.19×660005.0×1052789.444012.3134222.54×In Table5, the efficiencies of CPU subroutines and the corresponding GPU kernel functions for 3D mesh deformation (nb=5000,nt=5.0×105) are compared, in which subroutine No. 1 to No. 4 corresponding to the algorithms are presented in Section 3. As presented in Table 5, the CPU computational cost of subroutine-2 increased drastically, and the GPU kernel function can accelerate this subroutine by more than 500 times for this case.Table 5
Efficiency comparison of CPU subroutines and GPU kernel functions for 3D RBF mesh deformation (nb=6000,nt=5.0×105).
Subroutine noCPU (s)GPU (s)Speed up114.71671.403210.488×22287.954.2662536.297×3474.16056.820469.521×412.61730.411930.632×


## 4.1. Mesh Deformation of Flow Field around an Airfoil

In this subsection, we perform a series of test cases for mesh deformation of flow field around an NACA0012 airfoil with chord lengthc=1.0, and the boundary of far field is a circle of radius R=5.0, as shown in Figure 2. During the process of mesh deformation, grid nodes on the upper and lower chamber of the airfoil move according to (9)Δy=αx2,where Δy is the displacement along the y direction of boundary grid nodes on the airfoil, x is the x coordinate of these nodes, and α=0.8 is a prescribed coefficient. The deformation is carried out for 10 steps which leads to maximum displacement Δymax=8.0 at the trailing edge of the airfoil. There are four meshes of different resolutions which are investigated, as listed in Table 2. The initial and deformed mesh with number of grid nodes nt=1.0×104 and boundary nodes nb=200 are shown in Figure 2. As shown in the figure, the mesh is deformed smoothly without negative volume cells. The change of mesh quality after deformation is not examined in this paper since it has been studied in many literatures [Boer et al. [1]; [22]].Figure 2
Deformation of the mesh around an airfoil. (a) initial mesh. (b) Final mesh using RBF interpolation. (c) Close up of the initial mesh. (d) Close up of the final mesh using RBF interpolation.
(a)(b)(c)(d)Table 2
Efficiency comparison of the serial CPU code and GPU parallelized code for 2D RBF mesh deformation.
Mesh nonbntΔtCsΔtGsSpeed up1501.0×1030.01020.19000.053×21001.0×1030.00990.19500.051×32001.0×1040.27940.19701.51×44001.0×1040.28500.20451.39×54001.0×10518.74162.01729.29×68001.0×10519.06511.98349.61×In Table2, the execution time of GPU-accelerated code and the corresponding CPU code for different movement types and resolutions are presented, where nb is the number of boundary nodes, nt is the total number of grid nodes, and ΔtC and ΔtG are the measured execution time spent on CPU and GPU for single deformation step, repectively. As shown in Table 2, in the case of small quantity of grid nodes (around 1.0×103), the execution time of CUDA code is longer than the CPU code. We analyse the reason that GPUs cannot be fully exploited in this situation, and the execution time of memory copy shares a large proportion. When the total number of grid nodes increases to the magnitude of 1.0×104, the GPU code becomes faster than the CPU code. As the total number of grid nodes grows to the level of 1.0×105, more than nine times acceleration was attained for the GPU code. Also, we found that both the CPU and GPU code are not sensitive to the number of boundary nodes in these cases.In Table3, the efficiencies of CPU subroutines and corresponding GPU kernel functions for 2D mesh deformation (nb=800,nt=1.0×105) are compared, in which subroutine No. 1 to No. 4 corresponding to the algorithms are presented in Section 3.Table 3
Efficiency comparison of CPU subroutines and GPU kernel functions for 2D RBF mesh deformation (nb=800,nt=1.0×105).
Subroutine noCPU (s)GPU (s)Speed up11.40320.19827.080×23.85771.20353.205×313.58580.521926.031×40.21540.05983.602×


## 4.2. Deformation of the Mesh around a 3D Wing

In this subsection, we investigate test cases for mesh deformation of flow field around a 3D wing, the configuration of the wing is extruded from a NACA1410 airfoil with chord lengthc=1.0, and the boundary of far field is a cylinder of radius R=5.0 and length L=10.0. As presented in Table 3, a series of test cases with different movement types and resolutions are performed. Figure 3 shows the initial and final deformed meshes.Figure 3
Deformation of the mesh around a 3D wing. (a) Cross-section of the initial mesh. (b) 3D view of the initial mesh. (c) Final mesh with deflection. (d) Final mesh with translation and rotation.
(a)(b)(c)(d)In Table4, we compared the efficiencies of GPU-accelerated code and the corresponding CPU code for different movement types and resolutions, where nb is the number of boundary nodes, nt is the total number of grid nodes, and ΔtC and ΔtG are the measured execution time spent on CPU and GPU for the single deformation step. As shown in the table, for these 3D mesh deformation cases, tens to hundreds of times acceleration is achieved via GPU parallelization. We also found that for these 3D mesh deformation cases, as the number of grid nodes grows, the execution time of CPU codes becomes sensitive to the number of boundary nodes while the corresponding GPU codes do not.Table 4
Efficiency comparison of the serial CPU code and GPU parallelized code for 3D RBF mesh deformation.
Mesh nonbntΔtCsΔtGsSpeed up15005.0×1049.61041.78625.37×210005.0×10419.60971.769911.08×315002.5×105150.83806.663522.64×430002.5×105432.68666.573865.82×530005.0×105877.683912.328471.19×660005.0×1052789.444012.3134222.54×In Table5, the efficiencies of CPU subroutines and the corresponding GPU kernel functions for 3D mesh deformation (nb=5000,nt=5.0×105) are compared, in which subroutine No. 1 to No. 4 corresponding to the algorithms are presented in Section 3. As presented in Table 5, the CPU computational cost of subroutine-2 increased drastically, and the GPU kernel function can accelerate this subroutine by more than 500 times for this case.Table 5
Efficiency comparison of CPU subroutines and GPU kernel functions for 3D RBF mesh deformation (nb=6000,nt=5.0×105).
Subroutine noCPU (s)GPU (s)Speed up114.71671.403210.488×22287.954.2662536.297×3474.16056.820469.521×412.61730.411930.632×


## 5. Conclusion and Future Works

In this paper, we developed a GPU-accelerated mesh deformation method based on radial basis function interpolation. The GPU parallelized code is designed with CUDA C++ provided by NVIDIA, and the implementation of GPU parallelization is described in detail. Two test cases are carried out to validate the reliability of prescribed radial basis funtion interpolation and the performance of GPU-accelerated mesh deformation code. Over hundreds of times speed up has been attained for the GPU parallelized code corresponding to the serial CPU code. The code is proved to be efficient and robust. In future, the authors will integrate data reduction radial basis function in the GPU parallelized code to further improve the performance of mesh deformation.



---

*Source: 6018008-2022-10-13.xml*