Geometric Multigrid Method for Isogeometric Analysis

Houlin Yang,Bingquan Zuo,2,＊,Zhipeng Wei,2,Huixin Luo,2 and Jianguo Fei,2

1Key Laboratory of Metallurgical Equipment and Control Technology Ministry of Education,Wuhan University of Science and Technology,Wuhan,430081,China

2Hubei Key Laboratory of Mechanical Transmission and Manufacturing Engineering,Wuhan University of Science and Technology,Wuhan,430081,China

ABSTRACT The isogeometric analysis method(IGA)is a new type of numerical method solving partial differential equations.Compared with the traditional finite element method,IGA based on geometric spline can keep the model consistency between geometry and analysis,and provide higher precision with less freedom.However,huge stiffness matrix from the subdivision progress still leads to the solution efficiency problems.This paper presents a multigrid method based on geometric multigrid(GMG)to solve the matrix system of IGA.This method extracts the required computational data for multigrid method from the IGA process,which also can be used to improve the traditional algebraic multigrid method(AGM).Based on this,a full multigrid method(FMG)based on GMG is proposed.In order to verify the validity and reliability of these methods,this paper did some test on Poisson’s equation and Reynolds’equation and compared the methods on different subdivision methods,different grid degrees of freedom,different cyclic structure degrees,and studied the convergence rate under different subdivision strategies.The results show that the proposed method is superior to the conventional algebraic multigrid method,and for the standard relaxed V-cycle iteration,the method still has a convergence speed independent of the grid size at the same degrees.

KEYWORDS Isogeometric method;geometric multigrid method;reflecting matrix;subdivision strategy

1 Introduction

In recent years,as a numerical method for solving partial differential equations(PDEs),isogeometric analysis[1](IGA)has been widely concerned and studied.The basic idea of IGA:Based on the isoparametric element theory,geometric splines can not only describe the exact geometric modeling but also the solution space of numerical methods,which can unify the geometric model and the analytical model.Therefore,the spline information in CAD model is directly used as the initial mesh in IGA.Since it was proposed,most of the research on IGA has focused on parametric modeling and physical spatial discretization.However,when dealing with large-scale problems,it is still a problem to solve the system of linear equations which are obtained from the geometric discretization.The reason is that the direct solution is inefficient and usually requires targeted development of faster and more efficient iterative solver.

The multigrid method[2,3](MG)is often used as an efficient iterative algorithm to solve large-scale linear equations.The main idea is to eliminate high-frequency errors through fine mesh smoothing,coarse grid correction to eliminate low frequency errors,and an iterative strategy for fine grid and coarse grid cycles.According to the different mapping matrix generation methods,the MG can be divided into two types:geometric multigrid method(GMG)and algebraic multigrid method(AMG).Usually the AMG is more widely used then the other because it obtains different level meshes in an algebraic way instead of constructing them,which may be complicated in some cases[4].However,the IGA directly generates the hierarchical mesh through the refinement method(p-refinement,h-refinement,k-refinement),which encourages us to implement the GMG in IGA.

It seems natural to extend these methods to IGA,and several remarkable works has been made in recent years.MG for IGA based on classical concepts have been considered in[5,6],and a classical multilevel method in[7].An approach to local refinement in IGA based on a full multigrid method(FMG)can be found in[8].The related studies mainly include Gahalaut et al.[9]applied multigrid to IGA,and did not propose the concept of mapping matrix;Hofreither et al.[10]proposed a robust error approximation.These studies above mainly focus on the implement of MG in IGA.And the efficiency in the adaptation process has been studied.Other research mainly includes Guo[11]to improve efficiency through parallel computing and multigrid method,Liu et al.[12]combined conjugate gradient method with multi-grid method,and Chen et al.[13]based on geometric analysis such as least squares coupling multigrid method.Recently,in[14]some progress has been made by using a Richardson method preconditioned with the mass matrix as a smoother(mass-Richardson smoother).Similar techniques have been used to construct symbol-based multigrid methods for IGA,see[15,16].However,all the above research uses the coordinate mapping relationship in the process of geometric mesh construction to modify the traditional algebraic multigrid iterative process but ignore the stiffness information discretely obtained during isogeometric subdivision.In addition,the extension matrix transposition is directly used in the construction of restriction matrix,which will lead to the confusion and imbalance of mapping relationship.

Based on previous generations,this paper uses the geometric information discretized by isogeometric subdivision to couple the MG with IGA and its application.Especially for the numerical algorithms such as IGA,we use geometric method to optimize them in the shortcomings of AMG,and it has integrated into the process of FMG.

2 Isogeometric Analysis

2.1 B-Splines and NURBS Basis Functions

Like FEM,IGA also have shape functions,not the Lagrange interpolation functions but the spline basis functions.As the two common spline functions used in CAD,B-spline and NURBS are briefly introduced below.With a given node vector,the B-spline basis function[17]can be recursively defined as:

where the one-dimensional knot vectorξ=[ξ0,ξ1,...,ξn+k+1]is a set of non-decreasing sequences between 0 and 1.kis the degree of the basis function(p=k+1 is the order of the basis function).nis the number of basic functions,which is numerically equal to the number of control points.

If the basis function of B-spline is multiplied by an appropriate weightw,the non-uniform rational B-spline(NURBS)can be introduced.Its expression is:

2.2 The Framework of Multigrid Isogeometric Analysis

FEA and IGA have the same theoretical basis(i.e.,the weak form of partial differential equation).Take the common elliptic partial differential equation as an example:

whereΩ⊂R2is the open and bounded Lipschitz domain and the boundary is∂Ω,ΓD∪ΓN=Γ≡∂Ω,ΓD∩ΓN=∅;Given source termf:Ω⊂R2;Symmetric positive definite matrixk:Ω⊂R2×2;Dirichlet boundary valueg:∂Ω⊂ΓD;Neumann boundary valueh:∂Ω⊂ΓN.

It is assumed that the physical domainΩis represented by a parameter domainΩ0mappingUsing the control pointcij,Fcan be defined by NURBS as:

The weak form of the differential equation is obtained by multiplying with the trial function and integrating in the domainΩ.Letandthen the equation becomesa(p,v)= 〈f,v〉.According to Green’s method,the linear equations are obtained by discrete assembly:

The detailed process can be seen in[18,19].

2.3 Isogeometric Subdivision

Based on the requirements of light weight,modern CAD systems follow the principle of sufficient enough when characterizing models[20].So,the spline order in CAD is low usually,just meeting the basic characterization needs,and the node of the parameter domain is relatively sparse.But the order of the shape function(spline basis function)and the grid size(knot vector interpolation space scale)have a great impact on the analysis accuracy.As a result,the first step in IGA is subdivision,which works like the refinement in FEM,but with a different method.Therefore,it is necessary to introduce the subdivision process before constructing the isogeometric multiple meshes.

There are mainly three subdivision methods[21]in IGA,they arep-refinement,h-refinement andk-refinement.All these subdivision methods can keep the original shape the same.The changes are mainly at the spline order and knot vector.In particular,thep-refinement retains the original parameter domain grid and only increases the order of the basis function;thehrefinement does not change the basis function order but reduces the grid size to construct a finer parameter domain grid.As a combination ofp-refinement andh-refinement,k-refinement increase the order of basis function firstly,and then refines the knot vector.

Taken the one-dimensionalh-refinement as an example for introduction.Inserting the knott=[t1,t2,...,tm]⊂[τk,τn+1]into the knot vectorT=[τ0,τ1,...,τn+k+1]at one time can get the new knot vectorthen the corresponding coordinate calculation formula of the new spline control point can be obtained[22]:

Ifk=3(the cubic B-spline),Formula(6)can be written as follow:

where,the coefficients can be calculated as follows:

Duringp-refinement process,the degree of the spline curve is raised once,and the shape and continuity of the curve remain unchanged.Cohen[23]offered a solution to calculate the new control points:

K-refinement is a combination ofh-refinement andp-refinement.In order to ensure the continuity of elements,k-refinement is often defined as “first upgrade,then insert nodes” in IGA.This means,for the B-spline curve of orderp,the order is elevated toqbyp-refinement firstly,and then knots inserting is taken byh-refinement.

3 Multigrid Method

For FEM,most of MG method uses the algebraic multigrid method to construct interpolation operator and coarse the initial grid to get different hierarchical grid.In IGA,the geometric model is constructed from coarse to fine.There are different levels of grid naturally,and the mapping matrix between each level can be obtained from the subdivision process.In addition,each level has its own stiffness matrix and righthand vector.If this advantage is used to combine the multigrid method with IGA,the efficiency of solution would be improved.

3.1 Multigrid Mapping Matrix

3.1.1 Classical Algebraic Multigrid Mapping Matrix

In AMG,the main work is focused on constructing the interpolation matrix by coarsen the fine mesh,coarsening refined netΩmto get rough gridΩm−1(m=2,...,n).And the coarse gridΩm−1is a subset of the fine grid defined withCm,usingFmto define the remaining setΩm−Cm,whenthe pointiis strongly connected to pointj.Letdefine all the strong connection points of pointi,then the strong connection rough set isand the interpolation matrixrefers to the interpolation formula given by Chang et al.[24].

For MG,the mapping matrix is important.And its generation can affect the computational efficiency of MG directly.So,it is necessary to compare the generation rate of mapping matrix under different degrees of freedom.the time-consuming of mapping matrix generation(using Chang’s method)in Tab.1 shows that the time of generating the mapping matrix of the algebra multigrid rises sharply with the increase of degrees of freedom.Obviously,a faster constructing method for mapping matrix will be helpful for the implement of MG.

Table 1:The generation rate of mapping matrix

3.1.2 Suitable for Isogeometric Analysis of Multigrid Mapping Matrix

As discussed in Chapter 2.3,in the subdivision of IGA can generate the mapping relationship between new and old control pointsDmandDn.It is assumed that the relationship between(on fine grid)andu(coarse grid)also satisfies this relationship.According to the formula(7),the extension matrix from coarse to fine corresponding toh-refinement can be expressed as:

The corresponding linear relation is that the non-zero element in row vectorris numbered as the old control point number related to the new control point,which is usually a continuous interval segment.

P-refinement extended matrix can be obtained according to the calculation of Formula(8):

Therefore,the extended matrix ofk-refinement is:

In theory,the restriction matrixPfrom fine grid to coarse grid is the generalized inverse matrix of extension matrixR,but it is difficult to calculate the inverse matrix of extension matrixR,and there are stability problems and limitations in inversion calculation.Therefore,in this paper,transposed matrix is chosen as generalized inverse,P=RT.But in fact,the mapping matrixPobtained by transpose cannot accurately represent the mapping relationship from coarse grid to fine grid,so the mapping matrixαiis modified by parameterPhere:

Since the control pointDis a coordinate vector,the difference between the positive and negative coordinates will affect the effect of the correction.To ensure the positive correction of the parameterαi,the analysis model can be translated to the first quadrant to ensure that the control point coordinates are greater than zero.

In IGA,the mapping matrix obtained by different refinement strategies are shown in Tab.2.The detailed process can be referred to[25].

Table 2:Mapping matrix

Based on one-dimensional subdivision,two-dimensional subdivision only needs to subdivide inu,vdirections at the same time to get mapping matrixRu,iandRv,j,and then get the mapping matrixRξ,ηby global numbering of coarse and fine meshes.n,n′is the number of control points of the grid.Rξ,ηcan be represented by Eq.(13).

A time-consuming compare of mapping matrix generation between this method and AMG is given in Tab.3.With the freedom increase,more time can be saved in pre-processing stage of MG.

Table 3:The generation rate of mapping matrix

3.2 Multigrid Algorithm

3.2.1 Geometric Modification of Algebraic Multigrid Method

The multigrid method is divided into initial period and calculation stage.In the initial period,GMG is based on the subdivision of the solution domain to get the grid of different scales.AMG only uses the information of coefficient matrix in algebraic equations to construct multigrid,and the method is consistent in the iterative solution stage.In order to improve the implement of AMG in IGA,geometric information will be introduced to AMG,what is called Geometric modification of algebraic multigrid method(GMG).

1.Initial period

1).The nested grid layerΩ1⊂Ω2⊂...⊂Ωmis constructed by subdivision,and the extension matrixfrom layerm−1 grid to layermis stored.The mapping matrixfrom layermgrid to layerm−1 is modified according to Eq.(12).

2.Calculation stage

Algorithm 1:um ←GMG(Am,fm,um,d)1.If m=1,solve directly or precisely iteratively on Ω1,otherwise proceed to the next step.2.Pre-relaxation:um ←Sλ1(Am,fm,um)3.Coarse-grid correction:(a)Algebraic residuals rm−1=Pm−1 m(fm −AmAMG um);Geometric residual rm−1=Pm−1 m(fm −AmGMG um)(b)Take the minimum residual rm−1 from(a)(c) em−1 ←GMG(Am−1,rm−1,0,d)(d)Correction of residuals um=um+Rmm−1em−1 4.Post relaxation:um ←Sλ2(Am,fm,um)

The iterative process of multigrid method(GMG)is given above.For a given parameter spaceΩ1⊂(0,1)2,a rough gird for IGA can be built.Though refinement,multilevel grids can be constructed.,so the limit matrix of layermto layerm−1and the extension matrix fromm−1 tom.Parameterdin GAMG determines how to access each layer,as shown in Fig.1,d=1 is a V-type cycle andd=2D is aW-type cycle.The traditional AMG,the left term of the coarse grid is calculated using the mapping method,.In GMG,the correspondingAmof each geometric level can be calculated directly.That means,GMG in this paper introduces geometric information as a more selection of residual and selects the smaller residual as the next relaxation iteration by comparing the residual between geometric grid and algebra grid.

Figure 1:Multi-grid cycle

As for the iteration of relaxation operatorSv1(Am,fm,um),v1=v2=2 represents the number of times of relaxation.In order to get the iteration structure,Amis divided intoAm=Mm−Nm,and its iterative format can be written as:

Mmis selected to obtain various iterative methods for relaxation.Jacobi,Gauss–Seidel and SOR iterations are commonly used.For rough grid correction,both the residualand errorsatisfy the residual equationAmem=rm,soem−1can be solved by the same relaxation operatorS.

3.2.2 Full Multigrid Method

As known to all,Amandfmare discretized on each level grid in IGA,in order to make full use of these information,a full multigrid method(FMG)method suitable for IGA is established.

In the loop process of FMG,both AMG and GMG can be selected as solving algorithm on each layer,hence the accuracy of each layer can be guaranteed.Besides,and with the help of the extension matrix,the results can be interpolation to the previous layer as the initial solution which can be useful to reduce the iteration steps.As shown in Fig.2,the FMG process is as follows:

Figure 2:Full multigrid method

Algorithm 2:FMG algorithm 1.Using exact solution in the thickest grid u1=A1/f1 2.The initial solution of the second layer grid u2=R21u1 3.For the level m,the initial solution of the third layer um=Rmm−1um−1 4.Calculate um ←AMG(Am,fm,um) or um ←GMG(Am,fm,um)5.If m is the final level “exit”,else go to Step 3.

4 Numerical Examples

In order to verify the validity of the algorithm,this paper has carried out the verification on the Poisson equation and Reynolds equation.

According to the relationship between residual norm and error norm:‖rk‖ = ‖Aek‖ ≤‖A‖‖ek‖,the convergence of the residual and the error is consistent.

The relative error can be expressed as follows:

The convergence rate is defined as:

The algorithm is written in OCTAVE,the computer processor is Core i5-4570,with 3.20 GHz frequency.We contrast the convergence and calculation time of different subdivision strategies.The unified cycle format is V cycle,the number of grid layers m = 3,and the subdivision mode of coarse grid to fine grid isp1h1p2h2p3h3,corresponding to the 6-bit coding(such as 102026)shown in the figure below,respectively.Select the corresponding subdivision strategy from Tab.5(Reynolds)and Tab.8(Poisson),and the relaxation algorithm selects Gauss–Seidel iteration,and the number of internal iterations is controlled to bev1=v2=2.

4.1 Reynolds Equation Model

The piston-cylinder liner is used as a lubrication system,and the oil film pressure is derived based on the Reynolds equation.Fig.3a shows the piston liner system,and the area of lubricating film is shown in Fig.3b.Based on the Navier–Stokes equations and the mass continuity equation,the Reynolds equation[26,27]can be established as follows:

whereρis the lubricant density,pis the oil film pressure,ηis the dynamic viscosity,u0anduhare the oil film circumferential speeds on the piston and the cylinder liner surface,v0andvhare the film shaft speed on the piston and the cylinder liner surface,his the film thickness,tis time.

Since the lubrication area is symmetrical relative to the planeφ=0,only the right half needs to be considered.In general,the piston moves periodically,and the values of density and viscosity will change with the change of load value.However,in order to simplify the numerical model,the values of density and viscosity are taken as constant.The calculation parameters of oil film pressure are shown in Tab.4.

Figure 3:Piston-cylinder liner lubrication system(a)piston-cylinder system(b)oil film lubrication area

Table 4:Calculate parameters

The motion parameters of the piston-cylinder liner model have been given above,and the relevant detailed derivation process is referred to[28,29].

Table 5:Number of control points from MG based on different refinement strategies

When the Reynolds equation is discretized,the knot vector of the initial parameter domain isv=u=[0 0 1 1].In the initial state,given differentp,h,k-refinement strategies will get different levels of grid control points.

In the case ofh-refinement,n is the thinnest grid,and the appropriate coarse grid is selected to make it have the fastest convergence rate.It can be seen from Tab.6 that whenkis constant,the grid size is independent of the convergence rate.Since the Reynolds equation is an elliptic equation of second order,so it has the fastest convergence rate atk=2.Whenkincreases,the convergence rate decreases sharply.

Table 6:Convergence rate

Fig.4 only shows the comparison betweenh-refinement andk-refinement in order to explore the influence of different refinement strategies on the convergence rate.The broken lines in the figure indicate the number of subdivisions of each layer of the grid.For example,252627 represents the second-order fifth-level grids,the second-order sixth-level grids and the second-order seventh-level grids.The figure below is similar.The traditional AMG and GMG are used to compare the convergence of different degrees of freedom and spline degreek.The convergence rate of AMG and GMG iterations are shown in Fig.5.Tab.7 shows the CPU calculation time of the AMG algorithm and the GMG algorithm.It can be seen that:when the number k is constant,the convergence rate of AMG and GMG is almost the same under different grid sizes.Due to the different ways in which the coarse grid matrix is generated,the non-zero elements of the matrixAare different.As shown in Fig.6,nz is the number of non-zero elements.The band length in the GMG matrixAis reduced globally and locally,and the number of non-zero elements is less than that in the AMG matrix.Therefore,when the degree of freedom is large,the GMG of the natural discrete matrixAcan save more time and computing resources.Meanwhile,as the right items of different levels of are also produced in the discrete process,FMG can be used to calculate the original right items in the current layer.It can be seen from Fig.7 that the full multiple cycles have superior initial values in the iterative calculation process,and the desired error accuracy can be obtained in less iterations.

Figure 4:Relative error of refinement strategy(a) H-refinement relative error(b) K-refinement relative error

Figure 5:Relative error between AMG and GMG(a) k=1(b) k=2(c) k=3(d) k=4

Table 7:CPU calculation time of AMG and GMG

Figure 6:Coarse grid matrix A(a)GMG(b)AMG

Figure 7:Relative error between FMG and GMG

4.2 Poisson Equation Model

The equation of the two-dimensional Poisson problem with Dirichlet boundary conditions in the domainΩwith the following expression:

The domainΩis the quarter-ring of the inner radiusRi=1 and the outer radiusRe=2 is shown in Fig.8.When the right term heat source issin(2arctan(y/x)))/(x2+y2),the analytical solution issin(2arctan(y/x)).

Figure 8:Two dimensional poisson heat transfer problem

For the parameter domain knot vector of two-dimensional Poisson equation isv=u=[0 0 0 1 1 1]in the initial state,the convergence and calculation time of subdivision strategy are compared under the working environment mentioned above.The related calculated data correspond to the Reynolds equation are listed below.

Table 8:Number of control points from MG based on different subdivision strategies

Table 9:Convergence rate

Figure 9:Relative error of refinement strategy(a) H-refinement relative error(b) K-refinement relative error

Tab.9 discusses the influence of the choice of grid level on convergence in the AMG algorithm under different subdivision situations.Fig.9 only shows the comparison betweenh-refinement andk-refinement in order to explore the influence of different refinement strategies on the convergence rate.It can be seen from the figure that the mapping matrix satisfies the mapping relationship between the control points when there is a large gap between the number of control points of the coarse and fine grid,whetherh-refinement ork-refinement.For the assumed scale series betweenanduare large,the linear distortion leads to the inaccurate mapping relationship,and the convergence rate is not converged as expected.Within a certain range of the number gap between the coarse and fine grid points,the convergence rate ofh-refinement is different in the early iteration,but the convergence curve tends to coincide in the later 10 iterations,while the trend ofk-refinement only convergence rate is the same.Therefore,the selection of the appropriate grid hierarchy is beneficial to the convergence of the multi-grid method.Fig.10 shows the relative error of the AMG algorithm and the GMG algorithm.Tab.10 shows the CPU calculation time of the AMG algorithm and the GMG algorithm.Although the convergence rate of the AMG algorithm and the GMG algorithm is almost the same in the verification of the Reynolds equation,the GMG algorithm still saves more time in the initial period calculation.

Figure 10:Relative error between AMG and GMG(a) k=2(b) k=3(c) k=4

Table 10:CPU calculation time of AMG and GMG

Fig.11 shows the convergence comparison between FMG and GMG under the Reynolds equation model.It can be seen that:at the start stage,with the help of the good initial value,the FMG algorithm has a better accurate result,but the convergence tends to be consistent with others in the later stage.

Figure 11:Relative error between FMG and GMG

In this paper,the algebraic multigrid method and the geometric analysis of the natural discrete stiffness matrixAare used instead of the transpose generation matrixand the geometric multi-grid method and the algebraic multigrid method are compared.Numerical results show that:convergence rates of the two methods keep consistent,but GMG reduces the time cost at the initial period.For the subdivision strategies of different levels(3 levels),the convergence rate ofh-refinement intermediate layer is faster when it is closer to the fine grid,and the convergence rate ofk-refinement intermediate layer is faster when it is in the middle of coarse and fine grid.When FMG is applied to IGA,its convergence rate is faster than that of GMG.

5 Conclusion

Based on the framework of IGA,the analytical models of Reynolds equation and Poisson equation are established,and multigrid is applied to the IGA.According to its solving process and method,a calculation program was developed in OCTAVE(FMG based on GMG modified the traditional AMG).The convergence rate of the iteration is independent of the size of the discrete mesh,so the multigrid method has the optimal computational complexity.When the degree of freedom of fine grid is determined,the degree of freedom of coarse grid can be reduced and the boundary information can be quickly spread to the whole,which can not only accelerate the convergence rate but also save the cost of calculation time.However,the selection of grid level also affects the convergence accuracy.In particular,the selection of the middle level has a greater impact on the results,and the performance of the middle grid close to the top grid is better.As for the combination of FMG and IGA,the results show that the algorithm can make full use of IGA to discretize the geometric information of each layer.Compared with the V-cycle,it has a faster initial convergence speed,but it is consistent with conventional algorithms in the later stage.How to make better use of the existing information to accelerate calculation efficiency and improve calculation accuracy is the direction we will continue to study in the future.

Funding Statement:This research was supported by the Natural Science Foundation of Hubei Province(CN)(Grant No.2019CFB693),the Research Foundation of the Education Department of Hubei Province(CN)(Grant No.B2019003)and the open Foundation of the Key Laboratory of Metallurgical Equipment and Control of Education Ministry(CN)(Grant No.2015B14).The support is gratefully acknowledged.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.