Reusing the Evaluations of Basis Functions in the Integration for Isogeometric Analysis

: We propose a new approach to reuse the basis function evaluations in the numerical integration of isogeometric analysis. The concept of reusability of the basis functions is introduced according to their symmetrical, translational and proportional features on both the coarse and refined levels. Based on these features and the parametric domain regularity of each basis, we classify the bases on the original level and then reuse them on the refined level, which can reduce the time for basis calculations at integration nodes. By using the sum factorization method and the mean value theorem for the integrals, a new integration method with high integral efficiency is proposed. We validate the proposed method by some structural analysis problems in domains with different dimensionality. Comparing the numerical result accuracy and the time cost of the proposed integration method with the full Gauss integration quadrature, it turns out to be very promising.


Introduction
The efficient quadrature for isogeometric analysis, which was initially considered in Hughes et al. [Hughes, Reali and Sangalli (2010)], is a hot research topic in recent years [Wang, Xu and Pasini (2017)]. However, most of the works so far in isogeometric analysis still use the "full" Gauss quadrature on each element that is carried over from standard finite element analysis. It is no doubt that this integration method can give high-precision numerical results, but there still exists a non-ignorable problem: the number of Gauss quadrature points grows exponentially on each refined level and thus much more time is spent on the integral computation for the stiffness or mass matrix generation. To minimize the integration points for high efficiency, many integration methods are proposed. Hughes et al. [Hughes, Reali and Sangalli (2010)] proposed the "half-point rule" based on the precise smoothness of basis function across element boundaries, which recalculate the quadrature points on the macro-elements. And the number of quadrature points in this proposed method is less than the full Gauss method [Wang, Arabnejad, Tanzer et al. (2018)]. This integration rule is verified by some structural mechanics and fluid dynamics problems. Compared with the standard Gauss method, this method has higher integral efficiency. Takacs et al. [Takacs and Jüttler (2011)] discussed the special case of singularly parameterized NURBS surfaces and volumes represented by non-quadrangular or non-hexahedral domains without splitting, and presented an approach for verifying the integrability of these singular parameterizations and a method for modifying the test functions when a diverging integral is tested. Auricchio et al. [Auricchio, Calabro, Hughes et al. (2012)] developed several new quadrature rules for isogeometric analysis taking into account the features of basis functions, which are tested by some boundary value problems. Rypl et al. [Rypl and Patzák (2012)] have compared the computational efficiency of four available numerical quadrature in IGA: GSR (Gauss Standard Rule), GBE (Gauss rule on Bezier Elements), GPS (Gauss rule with basis functions Precomputed for all knot Spans), and HPR (Half-Point Rule), and they have got the conclusion that the HPR scheme is more appropriate than others. Schillinger et al. [Schillinger, Hossain and Hughes (2014)] proposed some non-tensor-product structure monomial quadrature rules based on the tensor-product Gauss and Gauss-Lobatto rules. It enjoys the same accuracy and stability behavior as the full Gauss quadrature with fewer integration points. Adam et al. [Adam, Hughes, Bouabdallah et al. (2015)] proposed a new approach to construct selective and reduced integration rules and developed a robust and computationally efficient local algorithm for computing the quadrature points and weights element by element. Calabro et al. [Calabro, Sangalli and Tani (2017)] proposed a new algorithm using no-element-wise assembling loop. Instead, this proposed method loops over the matrix rows based on a specifically designed weighted quadrature (WQ) rule. According to the tensor product structure of the B-spline basis functions and the sum-factorization implementation, this quadrature method also has high integral efficiency. Mantzaflaris et al. [Mantzaflaris, Jüttler, Khoromskij et al. (2017)] introduced a new quadrature method which reduces demanding multi-dimensional quadrature operations on tensor-product B-splines to inexpensive one-dimensional operations on univariate B-splines. For the refinement or the hierarchical refinement, each basis function must be recalculated on each level in traditional integration methods. Their low quadrature efficiency is a bottleneck for IGA, especially when a hierarchical refinement is adopted [Wu, Huang, Liu et al. (2015)]. Actually, the basis functions on the refined level have a shape similar to those on the original level, except for the width of their spans. And the values of basis functions on higher levels can be obtained from those at lower levels or the original level through transformation. According to this transformation concept, we could reuse lower-level basis functions without recalculating the basis functions on each refinement level. In Auricchio et al. [Auricchio, Calabro, Hughes et al. (2012)], the entire basis functions have been divided into three categories: the boundary functions, bubble functions and transmission functions. However, few works discuss the quadrature based on these transform features. Actually, our verification shows that the integration efficiency can also be improved by exploiting these basis function features. The structure and content of this paper is organized as follows. Section 2 begins with a brief review of the basis functions for splines and follows up with the discussions on their reusability and the descriptions of the basis-function classification based on the reusable rules formed according to the multiplicity of the nodes in the definitions of individual basis functions. The definition of reusable basis function is also discussed in this section. In Section 3, the Gauss integration method is introduced first, and then the extended Gauss integration method based on three features of basis functions is proposed. In addition, the integration points and weights under given knot vector is given there. Section 4 presents three examples, to which the proposed method is applied, and the comparative study that discusses about the computational efficiency and precision of the proposed method versus the standard Gauss method and the weighted quadrature method. Finally, we finalize the paper with a summary of some open problems related to this research in Section 5.

The redefinition of basis functions
We start with the discussion on three features of basis functions. Then, we propose a classification method for basis functions on entire refinement levels based on the features and their node multiplicity as well. At last we define the calculation method for basis function of NURBS according to the proposed classification method. For further details on the features of basis function and refinement, we refer to the fundamental works of Hughes and his co-workers [Hughes, Cottrell and Bazilevs (2005)]; Schillinger and Rank (2011);Schillinger, Evans, Reali et al. (2013); Bornemann and Cirak (2013)]; Scott, Thomas and Evans (2014); Auricchio, Calabro, Hughes et al. (2012); Wang, Wang, Xia et al. (2018)].

Three features of basis functions
The basis functions of B-spline curves with degree + ∈ Z p are defined from a non-decreasing knot vector The basis functions are piecewise polynomials on their parametric spans. At the two end spans, the corresponding basis functions comply with the symmetry. This feature is called symmetry rule in this paper. In addition, most of basis functions with the same multiplicity in their parametric domain have the same shape though they are related to different nodes with various spans, and these basis functions can be gotten from one function through transformation. This characteristic is called translation rule here. For the basis functions of univariate B-splines with degree p=3, knot vector {0 0 0 0 1 2 3 4 5 6 6 6 6}  , , , , , , , , , , , , Ξ = shown in Fig. 1(a), the two corresponding rules are illustrated in Fig. 1 For the basis functions from the refinement at the midpoints of intervals, there exists a ratio of basis functions between the coarse level and refined level. This scaling relation between basis functions at two levels is called scaling rule in this paper, which is shown in Fig. 1

The classification of basis functions
The sharpness of a basis function depends on the knot multiplicity of the support knot spans supp i N  . In this paper, the knot multiplicity is the total number of multiple knots in a specific support knot span. It should be noted that the multiplicity of each knot is different in the support spans of different basis functions. We suppose that the multiplicity of each support span is identical or symmetric, the corresponding basis function is the same. For reusing the basis functions based on the three rules in Section 2.1, we propose a reclassification of the basis functions according to the multiplicity of support knot span. We assume that the similar multiplicity of support spans supp k i N  of the refined basis functions can be found out on the spans 0 supp i N  of the original level. We only classify the original knot vector, and the refined basis function can be classified according to the multiplicity of support spans. Here, we take the original knot vector {0, 0, 0,1, 2, 3, 4, 5, 5, 5} Ξ = and its refinement knot vector {0, 0, 0, 0.5,1,1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5, 5} Ξ = for p=2 as an example to illustrate our classification method. The classification results are shown in Tab. 1. From Tab. 1, all the basis functions are divided into three groups based on the multiplicity of the knots on their support spans. For the basis functions on the same level, they only meet the translational or symmetry rule. But for the basis functions on the different levels, they satisfy the scaling rule in addition to the above two rules. Here, we propose the following classification principles: (1) Choose each support span on the refined levels and calculate the multiplicity of each knot node. (2) Compare them with those for the bases on the initial level one by one; if the multiplicity vectors are equal or symmetrical, their corresponding bases belong to the same class. For the spans [0 0 0 1] and [4 5 5 5] on the original level and the span [4.5 5 5 5] on the refined level, the corresponding multiplicity vectors are [3 3 3 1] , [1 3 3 3] and [1 3 3 3] respectively, which meet the requirements of symmetric rule, and so they belong to the same group. However, for span [0 1 2 3] on the original level and [0 0 0.5 1] on the refined level, their multiplicity vectors are different ( [1 1 1 1] vs. [2 2 1 1] ), and it is necessary to classify them into different groups (see Tab. 1). Here, for simplification, we suppose that all level refinements are conducted by subdividing all the knot span elements of the knot vector at their midpoints. The knot interval length on the original level is the same and the length among different refined levels is proportional to each other.

The transformation of basis function according to the classification rules
Based on the classification method discussed above, we select one basis function in each group to redefine the others. Here, this selected basis function is called reference basis R N , and its corresponding support domain is called reference span supp R N  . Generally, each group has only one reference basis and one reference span. Actually, for the basis functions spans in a same group, there exists an affine transformation between the defined reference span and others: where the coefficient C is a constant, and the variables µ and β are the first node value of the reference span and the mapped span respectively. Choosing the support span of . Note that it is wrong to confirm the variable β in Eq.
(2) using the above proposed method directly for the symmetric spans and they must be inverted at first and then calculate µ and β according to the corresponding reference span.
The coefficient for the proportional spans where k is the refinement level and the minus indicate that the selected span is symmetry with the reference span. Now, all the basis functions can be expressed by their corresponding reference bases through the affine transformation of their spans. And the value of basis function at a point is equal to that of reference basis functions. As for the evaluation of the compared basis function, it is important that the right reference basis is chosen. Here, we follow a "left" principle to select them: trying to select the basis function and its span located on the left end of the original level as the reference basis and span. Usually, these selected spans always contain node 0, and it is convenient to confirm the affine transformation for them. Based on affine transformation, all the basis functions ( ) N u can be re-expressed in terms of the reference basis ( ) This formula defines the translation relation of basis functions from the reference span to the other spans. However, in the process of calculating, the calculated span has always been known beforehand, and it is necessary to confirm its inverse transformation that projects the selected span to its corresponding reference span. This inverse transformation can be written in this way: So the basis function on the selected span can be expressed by the reference basis function: and the corresponding derivative has the form as follows  According to the knot multiplicity of their reference spans, the basis functions on these two levels are classified into four groups shown in Fig. 3(b) and Fig. 3(c). The variables in Eq. (4) are given in Tab. 2. Due to the chosen reference basis functions are distributed on the left of the original parametric domain, which contains node 0, so the variable µ is equal to 0 in general. In addition, 0 C < indicates that the span is symmetrical to the reference span.

The basis functions of NURBS
The basis function can be re-expressed by a reference basis via Eq. (5). The value of basis function at a point can be calculated by evaluating its corresponding reference and performing the inverse affine transformation Eq. (4). And the bivariate basis functions of NURBS can be written with the reference bases: where 1 2 ij ϖ is the control points weights.
If 1 ij ϖ ≡ , this formula can be simplified as follows: The partial derivatives for Eq. (7) are shown as: In the above equations, the values of basis functions are expressed with the reference basis functions. It is possible to reuse the evaluations for reference bases in the calculation of the basis functions according to this feature.

The numerical integration method
Nowadays, the integration efficiency is still a bottleneck for IGA. In this section, we propose a numerical integration method according to the concept of reference basis functions defined in Section 2, which takes less time in the integration compared with the "full" Gauss quadrature.

The quadrature based on the reference bases
The tensor product structure of B-splines basis functions can be expressed as [Antolin (2015)]: And its gradient can be written as follows: where In Antolin et al. [Antolin, Buffa, Calabro et al. (2015)], the integration formulation for the mass matrix and the stiffness matrix can be obtained according to the univariate quadrature rules in nested integrals: u N u c u du N u N u N u N u N u N u c u ,u ,u where ( ) c u is a factor about the coefficients of the investigated partial differential equation and the problem geometry, i.e., the Jacobian of the geometry. From the extended basis function defined in Section 2.3, each basis function in Eq. (12) and Eq. (13) has its own reference basis. They can be expressed using the reference basis as follows:  N u D N u D N u D N u D N u D N u c u u ,u u are the univariate reference basis. According to Eq. (11), an integral on a multi-dimensional domain has become a product of serval integral parts on one-dimensional domain. The format of each part is similar to each other.

The integrand in quadrature 3.2.1 The quadrature derived from Gauss method
From the principle of Gauss integration, where i Q and i W are the integration points and weights respectively, n is the number of Gauss points. According to Eq. (14), there are similar integral forms based on the reference basis functions defined in Section 2.3, and hence one reference basis might have been calculated several times in an integral operation. However, this repeated calculation has two disadvantages in the integral process. One is the low efficiency of integral calculation; another is a waste of a lot of computational recourses. To overcome this dilemma, we extended the Gauss method to a weighted Gauss integration.
Here, we define a new integration method similar to the general Gauss integration: where ( ) g x is a polynomial weight function. However, in this new method this function is associated with reference bases in the way ( )= ( ) ( ) . Note that ( ) g x does meet the following condition: This proposed integration method, which is called the extended Gauss integration method (EGM), is different from the traditional Gauss method shown in Eq. (15). The first is the integration polynomial that is defined by two basis functions. The second is the integration domain that is alien to the normal span [-1,1] of the traditional Gauss quadrature. In this proposed integration method, this domain is related to the function ( ) g x , that is, the integration domain is determined by the intersection of spans of the two basis functions in ( ) g x . Although it can be easily converted to a certain domain using the three features defined in Section 2.1, it needs more transformation parameters.

Combinations of the reference bases
In the previous section, the integral polynomial is a compound of two arbitrary reference bases. It is known that the number of reference bases is limited. It is possible to calculate the integration on all defined domain through this characteristic. However, due to the two differences mentioned above, it is not easy to deal with ( ) g x in the same way as traditional Gauss. For reusing the reference bases in our new integration method, the combinations of the reference bases are studied. In the integration process of IGA, there are four situations for each combination (such as When we calculate the integration of mass matrix, only 1 ( ) g x needs to be chosen. For the computational mechanics or the fluid problems, it is probably to select the 4 ( ) g x to calculate the stiffness matrixes or fluid equations. It is known that the basis function has its own influence domain, that is, a basis function may cover several elements. Although one can calculate the integration on several elements (macro-element in half point rule), it is hard to deal with the integration domain when the combinations of different reference basis functions are taken into account. In the process of integration calculation, our proposed method will continue to perform the integration on each element with the Gauss method which can be written as follows: where ( ) And the integration on the domain [0, 2] can be shown as follows: Due to ( ) g x is a piecewise polynomial on [0, 1] and [1, 2], this integration can be decomposed into the following form: According to Eq. (16) and the traditional Gauss method, each integration domain in Eq.
(23) has its own integration points i x and weights i ω , which can be obtained through the following equations: Here, we can obtain equations as follows (p=2):  u N u N u u g u u u N u N u u u x x g u u u N u N u u u x x g u u u N u N u  Using the same method, we can get their corresponding integration points and the corresponding weights for the all reference bases. With the experimental verification, for higher order of basis functions (less than 5), it can obtain the accurate integral results using only 2 integration points in one direction for each element. Actually, it is necessary to take three or more integration points in each direction according to the Eq. (23) to reduce the calculation error. Here, we consider a d-dimensional model problem on a single-patch domain, the degree of tensor-product space p and total dimension . NDOF For the standard Gauss method, each local stiffness matrix has dimension (p+1) 2d and each entry is calculated by quadrature on (p+1) d integration points. The total cost is about

The analysis of the integration condition for integral polynomial
According to the integration condition Eq. (17), the integral polynomial should be always great than or less than 0 on the whole integration domain [a, b]. However, this condition cannot be satisfied by all the situations of ( ) k g x defined in Eq. (18) at the same time. To overcome this problem, it is necessary to divide the integral domain into several sub-domains and keep the value of ( ) k g x always being negative or positive on each sub-domain. These sub-domains can be confirmed by solving the equations: In this equation, the integration points on some domains may be doubled. Due to the non-negative characteristic of basis functions, this situation does not appear when the integral polynomial is 1 ( ) g x .It is obvious that each calculated point is still located in its corresponding integration domain. However, there exists a large integral error using these integration points and weights without considering the condition of Eq. (17). So it should check the integrand carefully when calculating the integration points and weights.

The standard integration domain for the proposed integration method
In the above section, the calculated integration points and weights are present on the whole spans of reference basis. It is inconvenient to convert these reference basis spans to the actual integration domain, especially for the longer spans of the high-order basis functions. So it is necessary to define a standard domain for simplifying the calculation process and keeping accordance with the traditional Gauss integration method. Considering the features of knot vectors, this standard domain is set to [0, 1]. According to the three characteristics illustrated in Section 2.1, any span of basis function can be converted to this standard domain. This transformation is similar to the affine transformation Eq. (2) between coarse and refined basis spans that can be extended as follows: whereθ is a scaling ratio between the integration and the standard domain, η is a shift factor that equal to the lower bound of integration domain. Illustrating the integration points and weights on the standard domain, we take the knot vector in Section 3.2.2 as an example and give the points and weights for combinations of all reference basis functions, which are shown in Tabs. 3-5. In these tables, g2(x)/g3 (x) indicates that the integration points and weights are equal to each other under this combination situation.   17) is not considered . It is necessary to divide their corresponding domain into two subdomains according to Eq. (25), which means that there are four integration points on the domain of [0, 1]. According to the above table, for any interval integration [a, b], the integration points and weights can be expressed as follows: where i is the index of integration point/weight and C is the scaling factor in Eq. (5).

The integration method for multivariable bases
Substituting Eq. (5) and Eq. (23) into Eq. (16), one obtains 1 2 Note that C θ is equivalent for all reference bases in the same refined level. Eq. (29) can be expressed as follows according to the transformation of the upper and lower bound: Here, we assume that the knot vector is{0,0,0,0.5,1,1,1}, This numerical result is equal to the exact integral result. Whether the integration is for mass matrix or the stiffness matrix, there is the same calculation form shown in Eq. (28). Due to the reusability of the basis functions, it is not necessary to calculate every basis function on the whole integration domain. Similar to the traditional Gauss method, the stiffness and mass matrix in Eq. (14), can be expressed as  N u D N u c u du D N u D N u D N u D N u D N u D N u c u u ,u u DF u u ,u DF u u ,u DF u u ,u  c u u ,u DF u u ,u γ ⋅ respectively. The variable γ is the transformation coefficient between the current basis function and the reference basis function. For NURBS, due to the weights of control points are not all equal to 1, so the stiffness and mass matrix can be rewritten as follows (for 2-dimensional problem): ,u ,v ,v ,u ,v ,u ,u ,u ,v ,v A u,v A u,v s c u,v dudv B u,v B u,v B c u,v c u,v A u,v A u,v dudv A u,v A u,v dudv B B B c u,v B B c u,v A u,v A u,v dudv A u,v A u,v u,v A u,v m c u,v dudv B u,v B u,v c u v N u N u N v N v dudv B u,v x x where A and B are defined in Section 2.4. Compared with the weighted quadrature method, there are many differences in our proposed method. Firstly, the integration principle is different. The Gauss-lobatto rule based weighted quadrature is a fixed-point quadrature method, that the integration points and the corresponding weights are not changeable for the same test and trial functions. However, this proposed integration method is not a fixed point quadrature method, that is, the integration points and the corresponding weights are different for the different combinations of the reference bases. Secondly, the number of integration points in one element is different. For the weighted quadrature, only 2 integration points are needed in each direction far away from the domain boundary, while p+1 points are taken on boundary elements. For the proposed method, the number of integration points are the same for all elements on the whole domain if the corresponding basis functions satisfy the defined condition in Eq. (17). Thirdly, the computation complexity is different. It needs to confirm the value of trial functions on all integration points when calculating the integrand value defined in Eq. (4.1) in Calabro et al. [Calabro, Sangalli and Tani (2017)] for weighted quadrature. However, it does not need to calculate the value of basis functions for the integrand in the proposed method. Because the value of integration points and the corresponding weights already contains the results of the integrand.

Numerical examples
In this section, we present three application examples for the proposed integration method. Upon the results obtained for these examples, the accuracy, efficiency and convergence of the method are discussed. The results are produced with the programs developed in the environment of MATLAB, which implement the proposed integration method as well as other related methods for comparison. It is worthwhile to note that these programs utilize the same linear equation solver provided by MATLAB, which implements the algorithm of Gauss elimination with partial pivoting for the case of our problems.

A Poisson equation in 1D
Considering the following one-dimensional Poisson equation [Nguyen, Anitescu, Bordas et al. (2015)]: u u .The exact solution of this problem is We refine the definition domain to obtain the numerical solution of Eq. (31) using the standard Galerkin method. Here, we choose two integration points in each element in Tabs. 3-5 to obtain the approximate solution. For the higher accuracy of numerical solution, we refine the whole parameter domain nine times using h-refinement and estimate the relative error (between the exact solution and the proposed method, between the exact solution and the standard Gauss method) in In the proposed method, it is unnecessary to calculate the basis functions on the refined level when the stiffness matrix is generated. And the experiment results show that the computation efficiency of proposed method is higher than the Gauss (at least three times in Fig. 6). Here, the consumption time contains two parts: the pre-processing times and the assembling the global stiffness matrix time. The pre-processing time means that the consumption time of establishing the domains of reference basis functions and obtaining the integration points and the corresponding weighting coefficients. We also observe that it needs longer computation time as the degree of B-spline increases. For the B-splines with the same degrees, the consumption time of the proposed method is far less than the Gauss method.

A simple cantilever beam example in 2D
We test the proposed integration method with a simple cantilever beam problem in 2D, which is described in Fig. 7. The analytical solution of this problem is ( We analyze this problem using the above integration method to get the numerical displacement results on the whole domain. In this example, the initial bases degrees are set to p=2, the initial parameter domains are defined by the knot vectors , and the corresponding control point weights are all set to be 1. Fig. 8 shows the convergence of the solutions obtained with B-splines of 2 polynomial degree in the 1 H semi-norm and 2 L norm, respectively. This figure presents the error curves of three integration methods: Gauss integration method, the weighted quadrature method, and the proposed integration method. These methods are respectively denoted with Gauss, WQ and R-Gauss in this paper. It can be observed that the convergence rate of the R-Gauss is faster than the WQ method. WQ and R-Gauss methods are slightly more accurate than the Gauss method. It is noticed that the proposed extended Gauss integration method shows no advantage on first two refinement levels. Since the new reference basis functions are created in the process of subdivision refinement, the corresponding integration points and weights should be recalculated.

Hole within an infinite plate
The problem of a hole with an infinite plate subject to an x-direction traction x T at infinity has an exact solution is: This problem is described in Fig. 10(a), where symmetry is applied to the right and bottom edges and the traction boundary conditions are applied the left edges. The corresponding mesh and control points are shown in Fig. 10(b). R is the radius of the hole, L is the length of the finite quarter plate, E is Young's modulus, and v is Poisson's ratio. This example is analyzed using the Gauss integration method, weighted quadrature method and the proposed integration method. Based on the close-form solutions, the relative error norm of domain displacements is considered for convergence study: Here, the h-refinement is used to carry out the convergence study. Fig. 11 also shows the convergence comparison between the proposed integration method and the weighted quadrature method. From this figure, it is observed that the proposed integration method can converge at the same rate as the weighted quadrature method. It also has higher convergence rate as the degree of B-spline increases. Efficiency is another important measure for evaluating a solution method. The time consumption curves vs. the degrees of freedom under Gauss, weighted quadrature and R-Gauss integration method with B-splines of even polynomial degree (quadratic and quartic) and odd polynomial degree (cubic and quintic) are shown in Figs. 12(a)-12(b). It can be seen that the R-Gauss spends much less time in assembling the stiffness matrix than the other two integration methods. In order to keep the same integration points of the finite element method, we take 3 integration points to calculate this problem when the B-splines polynomial degree greater than 3.

Conclusions
In this article, we have discussed the basis functions classification method based on the concept of middle subdivision for the parameter domain, introduced the reuses of the basis functions among coarse and refined levels, and defined the reference basis functions as well as their spans, which can be utilized to decrease the calculation time in iterations. We have also proposed an extended Gauss integration method based on the Gauss method and the defined reference bases, and applied the combination of two reference bases to calculate the integration points and weights. The classification for basis functions turns out to be suitable for analysis model, because it is unnecessary to calculate every basis function on refined levels. The extended Gauss integration method appears to be applicable for IGA. It has the same accuracy as the traditional Gauss method while it bears less computational cost for the same DOFs. This proposed integration method performs well on several test problems, which shows that it is very promising for efficiency improvement in IGA. Nevertheless, there are some topics that will be investigated in the future. First, the classification for the basis functions at different levels should account for the multiplicity of knot nodes. An extension to complex situations such as multi-variable basis function is another possible topic requiring future research. Currently, the reference basis used is deduced from the relations between the corresponding spans and other spans. Further adjustments and more convenient formulas that are directly deduced from parametric knot vectors are under investigation. Second, each combination of reference basis function has its own integration points and weights, which is inconvenient for the calculation in integration. It would be helpful to carry out an elaborated study on the calculated integration points and weights for achieving higher integration efficiency. Third, due to the complexity of the parameter domain of the multi-patch and trimmed models, it is necessary to further study the integration method of the proposed method for the boundary basis functions. Here, the proposed integration method is only tested with some elliptic problems and we have not examined its effectiveness on hyperbolic and parabolic differential equations. Despite this, we still believe that the present developments have provided inspiration for studies on those more complex engineering design and analysis problems.