A Survey on Isogeometric Collocation Methods with Applications

: Isogeometric analysis (IGA) is an effective numerical method for connecting computer-aided design and engineering, which has been widely applied in various aspects of computational mechanics. IGA involves Galerkin and collocation formulations. Exploiting the same high-order non-uniform rational B-spline (NURBS) bases that span the physical domain and the solution space leads to increased accuracy and fast computation. Although IGA Galerkin provides optimal convergence, IGA collocation performs better in terms of the ratio of accuracy to computational time. Without numerical integration, by working directly with the strong form of the partial differential equation over the physical domain deﬁned by NURBS geometry, the derivatives of the NURBS-expressed numerical solution at some chosen collocation points can be calculated. In this study, we survey the methodological framework and the research prospects of IGA. The collocation schemes in the IGA collocation method that affect the convergence performance are addressed in this paper. Recent studies and application developments are reviewed as well.


Introduction
Isogeometric analysis (IGA) can be traced back to 2005, when Hughes introduced computational mechanics technology that advanced the seamless integration between computer-aided design (CAD) and computer-aided engineering (CAE) [1]. In CAD systems, curves and surfaces constructed in geometric design are generally expressed in the mathematical form of a non-uniform rational B-spline (NURBS) functions [2]. However, in CAE for classical finite element analysis (FEA) [3], a linear basis is often adopted as the basis function of the approximation space. Therefore, to perform finite element calculation of CAD models represented by NURBS curves and surfaces, these models need to be discretized into mesh representations; this process is called mesh generation [4]. Mesh generation often involves tedious manual interaction, which commonly accounts for about 80% of the entire finite element calculation process according to statistics [1].
To avoid mesh generation, Hughes et al. proposed IGA [1,5], which directly employs NURBS as the basis function of approximation space and the physical domain space, thus providing an effective tool for the finite element calculation of CAD models represented by NURBS curves and surfaces. Although FEA has widespread applications in simulation, the commonly used linear bases result in C 0 continuity between elements, and the approximate error caused by mesh generation can lead to low efficiency of the numerical results [5]. NURBS representations possess many benefits, which mainly include high-order continuity, a small number of degrees of freedom, local support, etc. First, as a tool in the CAD system, it can exactly express the geometry, such as conic curves and quadrics. Second, with NURBS expression, the geometry consists of control points and the corresponding NURBS bases defined in the parameter space whose number is much smaller than the number of nodes in a linear expression. Third, the local support property of NURBS can reduce the computational burden of assembly and facilitate parallel computing. Therefore, on the basis of these NURBS basis functions, IGA avoids the tedious and timeconsuming mesh generation that exists in classical FEA when CAD model simulation is performed. On the other hand, IGA considerably reduces the amount of computation and greatly improves the numerical accuracy of the solution. In conclusion, the merits of IGA compared with classical FEA are as follows: • Higher-order continuity across elements; • Higher accuracy without mesh generation; • Fewer degrees of freedom (DoFs).
As a general and powerful tool for analysis, IGA comprises two forms: IGA Galerkin and IGA collocation. This paper explores the more efficient and versatile IGA collocation methods for general PDEs.

Methods and Comparisons
IGA is implemented by two frameworks: IGA Galerkin [1] and IGA collocation [6]. In this section, we introduce these two methods and perform a comparison with related work.
The non-decreasing and non-repeated sequence of these knots is called the breakpoint sequence. Then, the ith NURBS basis function, N i,p (u), is recursively defined by the de Boor-Cox algorithm [2] as where B p i (u) is the B-spline basis function and ω i is the corresponding weight. In particular, when all the weights are 1, {N i,p } i=0,...,N are equal to B-spline basis functions {B p i } i=0,...,N . Figure 1 shows the B-spline basis functions under different orders with the same breakpoint sequence. Each parameter interval [u i , u i+1 ] has p + 1 related local bases, for which the property of local support is essential to the analysis. For high-dimensional NURBS surfaces and solids, their blending basis functions are calculated by their tensor product [2].
Because of the tensor product structure of NURBS, local modification is inconvenient. The construction of NURBS representation is no picnic, since it always involves a good parameterization, smooth joining and trimming, etc. However, the NURBS representation of the geometry has the privilege of high order continuity, which possesses better smoothness. And the property that the geometry is always lied within the convex hull of the control mesh results in reduced degrees of freedom for computation. Local support makes it also cheap and easy to modify the geometry. In addition, NURBS can represent quadrics and free-form shapes.

IGA Galerkin Method
When Hughes et al. proposed IGA, the calculation was implemented by the Galerkin method [3] through the weak form of the differential equation. Here, we call it the IGA Galerkin method [1]. In analysis, NURBS geometry Ω is regarded as the physical domain. Let the following Laplace equation over the physical domain be [5] ∆T + f = 0 in Ω, where the Dirichlet boundary Γ D , Neumann boundary Γ N , and Robin boundary Γ R form the boundary ∂Ω. The functions f , g, h, and r are all given, and T is the unknown function to be solved. Without loss of generality, suppose that they are all scalar functions. We call Equation (4) the strong form of the partial differential Equation (PDE).
Galerkin's method constructs the finite-dimensional approximations of infinite spaces S and V, i.e., S h ∈ S, V h ∈ V. These approximate spaces are spanned by the same NURBS bases as the physical domain Ω. Let the unknown numerical solution T r be a linear combination of these basis functions with N + 1 unknown coefficients, that is, Then, the weak form can be transformed into a(ω r , T r ) = L(ω r ).
By assembling the system of the weak form (8), we can solve the following linear equations and obtain the solution T r [5]: where q refers to the unknown coefficients, K = [K i,j ] is the stiffness matrix, and b = [b j ] is the load vector given by For elastodynamic problems [3], the weak form often leads to solving the following Lagrangian dynamics Equation [5,[7][8][9]: where q(t),q(t), andq(t) are the displacement, velocity, and acceleration vectors, respectively. Stiffness matrix K = [K i,j ] is related to partial derivative matrix B of the displacement and elastic coefficient matrix D. M = [M i,j ] is the mass matrix given by mass density ρ.
C is the damping matrix, which can be represented by the mass and stiffness matrices according to Rayleigh damping coefficients γ 1 and γ 2 . b(t) = [b i (t)] is the force vector given by external force f(t).
Given the local support of the NURBS basis, the assembly can be treated as the sum of the counterparts of elements according to their contribution. Thus, the calculation in the whole physical domain is transferred to the calculation on each element (Figure 2), which reduces the computational cost. The assembly involves integration on the physical domain, so some numerical integration technologies on the reference space arise to improve the calculation. Another issue is that the condition number of the assembled stiffness matrix can be large, especially in the presence of large DoFs. This issue can be alleviated by pre-conditioning operations [10][11][12][13] and high-quality parameterization schemes [14][15][16][17][18][19][20][21][22][23][24].

Isogeometric Collocation Method
When the degree of the NURBS basis used in IGA is higher than the order of the differential equation, the form of the equation still exists if the NURBS-expressed numerical solution is introduced into the differential equation. Then, we can select a number of points in the computational domain, introduce these points into the equation form above, and obtain a linear system. The unknowns in the equations are the unknown coefficients of NURBS functions. This is the isogeometric collocation (IGA-C) method [6] proposed by Auricchio et al. in 2010, where the points selected in the calculation domain are called collocation points. The convergence and consistency of some theoretical results of IGA-C have been examined and proven by Lin et al. [25,26]. Distinguished from IGA Galerkin, which focuses on the weak form, IGA-C directly deals with the strong form of PDE. More generally, the boundary value problem (BVP) can be summarized as follows [6]: where L represents a scalar differential operator, T is the unknown solution, and the given right-hand side f is called the load function. The boundary vector operator GT represents the boundary conditions. IGA-C solves the equations by introducing a set of collocation points {ū k } k=0,...,N int into parametric domainΩ and {ū k } k=N int +1,...,N int +N bd on the boundary of parametric domain ∂Ω. The corresponding points of the physical domain arex k = F(ū k ), k = 0, . . . , N int + N bd , where F is the NURBS mapping defined in Section 2.1. Thus, replacing the unknown T of BVP with the expression of T r at {x k } i=0,...,N int +N bd leads to a linear system.
The coefficient matrix of the linear system is called the collocation matrix. In the system obtained from the IGA-C framework, the number of equations is equal to the number of unknowns, i.e., N int + N bd = N. Direct solving and solving with iterative methods can be adopted to solve the linear system. Given that the partial differential operator acts on the physical coordinates, and the NURBS solution is a function of the parametric coordinates, the derivatives can be implemented by the chain rule. For instance, the derivatives of the 3D case involve the following calculations: where x i , x j refer to any one component of the physical coordinate x = [x, y, z] T calculated by the NURBS mapping where When the number of selected collocation points is larger than the number of unknown coefficients, i.e., N int + N bd > N, which indicates that the collocation matrix in the system is over-determined, the linear system can be solved in a least-squares manner [6,27,28]. On the basis of this idea, Lin et al. [27] proposed the isogeometric least-squares collocation method (IGA-L) and presented the conditions and theoretical proof for convergence and consistency. The numerical results of the IGA-L method validate the increased accuracy obtained by adopting numerous collocation points.
It is worth mentioning that the original application of IGA-C was on the physical domain constructed by NURBS in simulation, and later there were advanced IGA-C corresponding to more splines. Due to the exploit of NURBS, the privileges of IGA-C can be concluded as a high-order continuity of the solution and high efficiency for computation. Nevertheless, if the collocation points are not properly selected, the condition number of the coefficient matrix of the generated linear system may be very large, or even singular. Thus, the collocation schemes become crucial to the effectiveness and accuracy of IGA-C.

Comparison of IGA Galerkin and IGA-C
IGA Galerkin has stable calculation, a high convergence rate, and perfect theoretical analysis [1]. However, it requires a large amount of integration calculation, and the assembly of the stiffness matrix is highly complex. The integration efficiency of IGA was first considered by Reference [29]. Then, Gauss-Galerkin quadrature rules [30] and optimal and reduced quadrature rules [31] were developed. A promising approach for reusing the basis function evaluations in numerical integration has been recently presented to improve the accuracy and time cost [32].
The commonly used tensor product NURBS basis must be reparameterized globally when performing refinements, which adds a number of unnecessary DoFs to the system. To reduce the computational burden, Bazilevs et al. [33] explored T-splines with IGA for local refinement and provided the alternatives of the adoption of the bases. Several alternative spline technologies have been proposed as a basis for IGA. In addition to T-splines [34], hierarchical B-splines (HB-splines) [35][36][37][38], truncated hierarchical B-splines (THB-splines) [39], polynomial splines over hierarchical T-meshes (PHT-splines) [40][41][42], locally refined splines (LR-splines) [43], and hierarchical extension of analysis-suitable T-splines (HASTS) [44] are also capable of adaptive local refinement in IGA. By adopting non-polynomial functions, such as trigonometric or hyperbolic functions, generalized B-splines [45,46] and generalized T-splines [47] enhance the piecewise polynomial spline bases for analysis. For discontinuous problems, such as cracks, the extended finite element method (XFEM) [48] incorporates a discontinuous displacement field and asymptotic crack tip displacement fields into the displacement approximation through the partition of unity method. Developed from the concept of XFEM, IGA has been enhanced to extended IGA (XIGA) [49,50] to locally enrich the bases for simulating discontinuities and singular fields. In accordance with the crack location, the original bases in the area can be selected to split into several new bases multiplied by an enrichment function and then contribute to the overall approximation. In addition to crack problems, enrichment functions can be custom designed for other specific problems, such as imitating local refinement when using the splitting basis techniques in XIGA [51]. An overview of the computer implementation aspects of IGA can be found in Reference [52].
In particular, IGA Galerkin requires high-quality parameterization of the computational domain. If singular points exist, the system may suffer from collapse unless a special treatment is provided. Cohen et al. [53] proposed the concept of analysis-aware modelling to show that the quality of domain parameterization has a great impact on the analysis results and efficiency [54], and that the condition number of the assembled stiffness matrix affects the stability of the linear system to be solved [55]. The research on analysis-suitable planar [14][15][16][17][18][19] and volumetric [20][21][22][23][24] parameterization for IGA has gradually flourished. Moreover, Xu et al. attempted to decouple geometry and simulation and obtained satisfactory results [56]. They revealed that the bases for the solution space are not necessarily the same with the bases for the physical domain. Overall, constructing high-quality analysis-suitable parameterization from a complex CAD boundary is still an open problem.
In conclusion, the interests of IGA Galerkin are focused on the following aspects: • Numerical quadrature rules in the assembly; • High-quality parameterization for IGA-suitable analysis; • Construction of alternative basis functions for particular cases.
IGA-C was proposed for accelerating the solving of PDEs by avoiding the efficiency issue of numerical integration in IGA Galerkin. It focuses on solving the strong form of PDEs rather than the weak form in IGA Galerkin to increase the accuracy-to-computational cost ratio [6,36]. IGA-C considers the superior accuracy and smoothness of NURBS basis functions and the low computational cost, which is a compromise between accuracy and computational time. The assembly of the stiffness matrix becomes simple and easy, and the bandwidth of the matrix becomes narrow. Substituting the derivatives of the NURBS numerical solution at some collocation points into PDE leads directly to a linear system. The number of equations (i.e., the collocation points) in the system can be equal to or greater than the number of unknown coefficients of NURBS functions. The solution (or the least-squares solution) of the system determines the coefficients of the NURBS function as the numerical solution of PDE.
Compared with the IGA Galerkin method, the IGA-C method has many merits. First, the calculation amount is small, and the calculation speed is fast. Second, the form is simple and flexible, and the linear equations of the IGA-C are easy to construct. Lastly, even in the presence of singular points, IGA-C calculation can still be carried out as long as appropriate collocation points are selected.
To date, the most popular numerical methods in solving PDEs are finite difference method (FDM) and finite element method (FEM). FDM uses the difference to replace the exact derivative, which is fast but loses some precision in derivative approximation. FEM solves the variational form (weak form) equivalent to the original problem (strong form), in which the quadrature in the generated linear system is treated numerically. Both methods require mesh generation of the computational domain. The biggest difference between IGA-C methods and classical FDM and FEM is that IGA-C methods avoid mesh generation, because the computational field and the physical field are the same and are both represented by NURBS. Moreover, the degrees of freedom of IGA-C refer to the number of control points, which is much less than the degrees of freedom in FDM and FEM, which refer to the number of mesh nodes. Because IGA Galerkin also needs complex quadrature computation in weak form similar to FEM, whereas IGA-C methods select some appropriate collocation points by collocation scheme and substitute them into the strong form directly for solving, IGA-C is much easier and faster to implement than other numerical methods. With the high accuracy-to-computational-time ratio, the research on IGA-C methods and their applications is significant.
However, the theoretical system of the IGA-C method is far from being established. Lin et al. revealed the conditions for the consistency and convergence of IGA-C and IGA-L methods with the stable or strongly monotonic differential operator and provided the proof in their work [25,27]. In existing work, the convergence rate was mostly verified by numerical examples. Lin et al. also deduced the convergence rate and developed the necessary and sufficient conditions for the consistency of the IGA-C method [26]. The IGA-C method is consistent if and only if the differential operator in BVP and the interpolation operator on the load function are uniformly bounded when the knot grid size tends to zero, which advances the numerical analysis of the IGA-C method.
Moreover, the calculation accuracy and convergence rate of IGA-C are not high. These defects seriously restrict the further application and development of the IGA-C method. A few studies have focused on how to select ideal collocation points for the parameterization of a given domain to improve the convergence rate of IGA-C [28,[57][58][59][60]. In fact, many factors, such as the parameterization quality of the computational domain [61,62], the collocation selection scheme [28,[57][58][59][60], and the differential equation [63,64], affect the convergence rate and calculation accuracy. The current work does not consider the comprehensive impact of multiple factors on the IGA-C method. Therefore, the urgent tasks for current studies on IGA-C can be summarized as follows: • Clarification of the consistency and convergence conditions; • Design of efficient collocation schemes; • Improvement of calculation accuracy and speed.

IGA-CL and IGA-GL Methods
When solving an IGA problem, the exact solution is usually unknown, and the approximation error is difficult to estimate. The problem lies in the determination of how many degrees of freedom are needed to realize acceptable accuracy and when to stop the refinement. Motivated by this problem, the PDE LT = f with the strongly stable differential operator L can be solved with IGA methods by fitting the load function f [25,27,65]. The exploitation of the load function can control the error estimation. The approximation error of numerical solution T r no longer requires the exact solution T * of PDE as a prior. Without knowing T * when calculating the error T * − T r 2 , the approximation error takes good advantage of load function f , which is defined as According to the necessary and sufficient condition for the consistency of the IGA-C method [26], suppose L is the differential operator in BVP (13), I ρ is an interpolation operator satisfying I ρ f = LT r , T r is the IGA-C numerical solution, the IGA-C method is consistent if and only if L and I ρ are both uniformly bounded when knot grid size ρ → 0, i.e., LT * − LT r 2 → 0, when ρ → 0.
Because IGA-CL is actually a knot insertion of the initial parameter space for constructing numerical solutions in IGA-C, the essence of IGA-CL solution is the IGA-C solution. Then, with LT * = f holding, we can obtain Thus, the necessary and sufficient condition for the consistency of the IGA-CL method is derived as well. Specifically, in a typical IGA, the knot vector of the NURBS numerical solution is determined only by the physical domain. On the one hand, solving PDE with IGA is regarded as fitting load function f defined on the NURBS geometry Ω with the differential operator on the NURBS numerical solution LT r . On the other hand, in the area of data fitting in geometric design, the design of the knot vector has a close relationship with the NURBS functions to be fitted. Therefore, the detected feature points of the load function are integrated into the initial knot vector of the physical domain to construct the new knot vector of the numerical solution in the IGA framework. The concept of combining the information from the load function and physical domain to improve accuracy applies not only to the IGA-C framework but also to the IGA Galerkin framework, which are referred to as IGA-CL and IGA-GL methods, respectively.
Thus, given a PDE defined on the NURBS physical domain, we can always apply the feature detection strategy to the load function and generate the improved knot vector for the NURBS solution with high accuracy using IGA-GL and IGA-CL methods. The error estimation becomes easy to control without the exact solution.

Collocation Schemes and Convergence
The selection of collocation points influences the efficiency of IGA-C greatly, and it is the crucial part of the procedure because it directly affects the accuracy and convergence of the numerical solution. Some popular options, such as Demko and Greville abscissae [6,57], superconvergent (SC) points [28], Cauchy-Galerkin (CG) points [58], and alternating/clustered superconvergent (ASC/CSC) points [59], have been widely utilized in IGA-C problems since they were proposed. In this section, we introduce commonly used collocation schemes.

Demko and Greville Abscissae
In IGA-C, the commonly used collocation points are Demko [57] and Greville [6] abscissae. Demko abscissae are obtained from the extreme points of Chebyshev splines and can be calculated by the function "chbpnt" in MATLAB or by referring to the book by deBoor [66]. They are proven to be stable for any mesh and degree. Greville abscissae are acquired from the averaging operation of the knot vector. However, the convergence orders of both choices are O(h p−1 ) for the odd degree p of the NURBS bases and O(h p ) for the even degree p.
Assuming that the NURBS basis for the numerical solution is associated with an open knot vector U = {u 0 , u 1 , . . . , u N+p+1 } with multiplicity p + 1 at both ends, Greville collocation points are generated via an averaging operation as follows [6]: Figures 3 and 4 display the distribution of Demko and Greville points under odd and even degrees. Each knot interval includes at least one collocation point, and the number of generated collocation points is equal to the number of unknown coefficients N + 1. Except for points near the boundary, Greville points locate exactly on the knot positions for an odd degree, and Greville points locate exactly on the centre of the knot intervals for an even degree. Demko points are very close to Greville points but not quite equal.

Superconvergent Points
On the basis of superconvergent theory, Anitescu et al. developed superconvergent points (SC points) [28] to increase the convergence rate.
They attempted to find collocation points in which the collocation solution behaves similarly to the standard Galerkin solution. On the basis of the observation that the error in the second derivative for collocation approximation is the smallest at the collocation points, they sought for the collocation points in which the error in the second derivatives of Galerkin approximation is also small. The computed numerical solution is expected to be close to the Galerkin numerical solution up to higher-order terms as well. The problems are that the accuracy of SC points are not high enough and that there are nearly twice as many points as DoFs. Given that the SC points are all adopted as collocation points, a solution of the overdetermined linear system needs to be computed by least-squares approximation. Anitescu et al. improved the convergence order of the IGA-C solution to the optimal O(h p+1 ) for odd degree p and kept one-order suboptimal O(h p ) for the even degree. We present the distributions of SC points under odd and even degrees (e.g., p = 3 and p = 4) in Figure 5. Except for points near the boundary, at least two symmetric SC collocation points exist in each knot interval for the odd degree, whereas SC points are located at the centre and both ends for the even degree. The numbers of SC points are 2(N + 1 − p) + 2 and 2(N + 1 − p) + 1 for odd and even cases, respectively, which are greater than the number of unknown coefficients N + 1 and are solved in a least-squares manner.

Degree
Second Derivative Superconvergent Points

Cauchy-Galerkin Points
Cauchy-Galerkin points (CG points) have been proven to exist such that the IGA-C solution at these points can coincide with the IGA Galerkin solution [58]. CG points are selected as a certain subset of the zeros of the Galerkin residual, i.e., SC points (Table 1). Under some hypotheses, these approximated CG points achieve superconvergence of the second derivatives of the Galerkin solution. IGA-C with CG points has been verified to have the convergence order of O(h p ) for both odd and even degrees. Figure 6 depicts the distribution of CG points, which preserves one SC point in each knot interval while maintaining the global symmetry as well. Thus, the number of CG points is equal to DoFs N + 1.

Alternating/Clustered Superconvergent Points
Montardini et al. found a subset of SC points and obtained the optimal convergence order O(h p+1 ) for odd degree p [59]. They were inspired by the selection way of CG points and devised a scheme of alternating SC point selection (ASC), which chooses collocation points from every other SC point. The collocation points in the elements near the boundary can be chosen flexibly without affecting the convergence rate. Examples of the distribution of ASC points are shown in Figure 7. Furthermore, Montardini et al. attempted to acquire a globally periodic distribution of collocation points while maintaining the symmetry at the local element scale. The collocation points are selected as two SC points for every other element, which is shown in a clustered way. We call these clustered SC points (CSC points). Figure 8 illustrates an instance of the clustered scheme under odd and even numbers of elements for the odd degree. The collocation points near the boundary are determined by considering the DoFs, symmetry, and periodicity. To compare these collocation schemes, we present an example of numerical experimental results under different collocation schemes in Figure 9. Here, we solve the source problem in reference [26] defined on a quarter of an annulus with degree p = 3, and we display the absolute errors of the IGA-C numerical solutions with different collocation schemes after some h−refinements. Specifically, the relative errors in L 2 norm are shown in Table 2. Finally, we compare the convergence rates of the above-mentioned collocation schemes under odd and even cases and present the results in Table 3. Table 2. Relative errors in L 2 norm of IGA-C with different collocation schemes after twice h−refinements (p = 3).  Table 3. Comparisons of orders of convergence in L 2 norm [59].

Other Collocation Schemes and Open Problems
Wang et al. established a superconvergent IGA-C method collocating at Greville points that utilizes recursive gradients to replace the conventional gradients of the bases [60]. Marino et al. investigated the influence of different parameterization and knot placement techniques on the accuracy of collocation-based formulations [61]. A discussion of the convergence rates of the IGA-C method with various collocation schemes can be found in [59]. The open questions regarding the selection of collocation points can be summarized as follows: • How can a new collocation scheme be designed to improve the computational accuracy and convergence rate of the IGA-C method even for even-degree bases? • How can boundary conditions (Dirichlet, Neumann, mixed boundary conditions) be dealt with in the new collocation scheme? • How can a collocation selection method with the same convergence rate as the IGA Galerkin method be designed?
• How can a collocation scheme in the presence of singular points in the parameterization of the computational domain be designed such that the IGA-C method can perform an effective calculation with the highest possible calculation accuracy?

Applications
In the past 20 years, IGA-C has been widely applied in some scientific fields of physical simulations such as structural mechanics, elastic mechanics, fluid mechanics, acoustics, magnetics, shape and topology optimization, and so on. Due to the high efficiency of computation of IGA-C methods for solving PDEs, they have been verified to have good application prospects in these various fields, and we are convinced of their great potential in solving physical problems in other unexplored scientific fields in the future. In this section, we categorize and summarize the related applications with IGA-C methods. Readers may refer to the related references for more details.

Structural Mechanics
IGA-C was originally introduced to solve the thin structural problems described by Bernoulli-Euler beam and Kirchhoff plate models [67], and the cases under different boundary conditions have been discussed. For mechanical simulation of slender, elastic, and spatial rods and rod structures with multiple rods subject to large deformation and rotation, IGA-C solves the strong form of the equilibrium equations of forces and moments [68]. The centreline of a rod is represented by a NURBS curve, and the cross-section orientation is expressed as 4D NURBS curves based on unit quaternions. A kinematic model with a mixed isogeometric collocation formulation for shear-deformable beams has been presented and determined to be locking-free for any choice of approximation degree [69]. IGA-C provides a highly robust, efficient, and accurate results, and it is powerful in handling geometrically complex beam problems. The locking issue existing in Timoshenko beams and rods has also been addressed with mixed collocation schemes [70,71]. Structural dynamics with IGA-C consisting of explicit and implicit integration strategies were considered in [72,73]. The IGA-C method has now been more widely applied in structural mechanics [74,75].

Elastic Mechanics
For static elasticity and elastodynamic problems [76][77][78][79], IGA-C provides a combination of high-order accuracy, stability without hourglass modes, and high efficiency, because it requires a minimum number of quadrature evaluations. For nearly incompressible elasticity problems, Morganti et al. [80] employed IGA-C and achieved the same accuracy results as those for compressible elasticity via the mixed displacement-pressure method. Fracture problems with contact have also been considered in elasticity [81].

Fluid Mechanics and Mixed Problems
Paired with plane wave (PW) and generalized harmonic polynomial (GHP) enrichment functions based on XIGA [49,50], IGA-C has been utilized in acoustic problems to solve 2D problems for the Helmholtz equation [82]. A comparison of Greville, CG, and SC points has shown that the CG scheme is the best choice in terms of overall error, convergence rate, and ease of solving the linear system. The application of IGA-C in acoustic problems can also be found in References [83,84]. IGA-C has also been applied in mixed fluidstructure interaction problems [85]. For porous media, IGA-C simulates poromechanics and models the linear elastic response of a fully fluid-saturated porous medium [86], and for PDEs involving the fractional Laplacian operator, it solves the fractional porous media equation [87].

Shape and Topology Optimization
For the IGA Galerkin framework, the idea can be extended to the isogeometric boundary element method (BEM) to solve the potential flow in the shape optimization of hydrofoils [88]. It uses generic B-splines to generate hydrofoil shapes and the same basis of geometric representation to represent the velocity potential. For 2D low-fidelity aeroelastic analysis, isogeometric BEM is involved in the optimization framework, which is a closely coupled isogeometric potential flow model and isogeometric curved Timoshenko beam model [89]. Isogeometric topology optimization (ITO) in the Galerkin formulation was proposed to achieve an effective and efficient solution process [90][91][92]. An energy-based homogenization method (EBHM) to evaluate material effective properties is numerically implemented by IGA for the systematic design of 2D and 3D auxetic metamaterials [90]. For topology optimization of mixed problems, such as vibrating structures interacting with acoustic waves to minimize the radiated sound power level or multi-frequency acoustic topology optimization, IGA has shown its validity and effectiveness as well [93,94].
Shape and topology optimization in engineering can also be extended to the IGA-C framework. Ship hull optimization is an application of IGA-C [95], which represents the ship hull model with T-splines and employs the IGA-C hydrodynamic solver to calculate the ship wave resistance. The resulting boundary integral equation is numerically solved by a higher-order collocated BEM with IGA. The hydrodynamic solver and the ship model are then integrated in an optimization for local and global ship hull optimizations against the criterion of minimum resistance.
Given the high efficiency of IGA-C, researchers have proposed a new topology optimization method using IGA-C based on moving morphable components [96]. They provided the Young's modulus of the material and the density field by employing three new Ersatz material models based on uniform, Gauss, and Greville abscissae collocation schemes and obtained good performance.
IGA-C has likewise been utilized for shape optimization of non-linear curved 3D beams and beam structures subject to large deformations and rotations [97]. The centrelines of the initial and current beams are expressed by the NURBS curve, and IGA-C discretizes the strong form of the balance equations. Then, the positions of the control points of the NURBS curve are optimized, and the performance is validated.

Discussion and Conclusions
This study provides a systematic introduction to the isogeometric collocation method. Starting with IGA in Galerkin formulation and its pros and cons, we present the development of the efficient isogeometric collocation method in detail. First, we introduce and compare the two methods. Second, a new perspective to the IGA method is developed by fitting the load function of PDE. Current mainstream collocation schemes are also investigated, and the corresponding convergence rates are considered. Lastly, we review the previous work on various applications with the isogeometric collocation method. The research directions for future work are summarized as follows: • Quality optimization for computational domain parameterization.
Similar to how mesh quality in classical FEA has an important influence on calculation accuracy, the parameterization quality of the computational domain represented by NURBS is also one of the important factors that affect the calculation accuracy of the IGA-C method. Optimizing the quality and robustness of computational domain parameterization is crucial.

•
Design of the collocation scheme.
Currently, the number of collocation points can be equal to or greater than the number of unknown coefficients of the NURBS function, and the positions of collocation points are the most decisive factors that affects the calculation accuracy and convergence rate. The convergence rate of existing collocation schemes does not exceed the convergence rate of IGA Galerkin. Although some studies have proven that the collocation points with the same convergence rate as the IGA Galerkin method exist, its construction method is still unclear. In some special cases, singular points must be included in the computational domain parameterization. In this situation, the collocation scheme determines whether the IGA-C method can be effectively performed or not.
• Consistency, convergence, and convergence rate of the IGA-C method.
At present, the consistency and convergence of most collocation schemes have not been theoretically proven, and the convergence rate can only be verified numerically. Obviously, clarifying the consistency and convergence conditions and deriving their theoretical convergence rate formulas are fundamental for establishing the theoretical system of the IGA-C method. Numerical verification has shown that the convergence rates of most collocation schemes for odd and even degrees of NURBS functions are different. Notably, the convergence, consistency, and convergence rate of the IGA-C method are closely related to various factors, such as differential operators and right-hand side load functions, and these factors need to be considered comprehensively.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

IGA
Isogeometric analysis IGA Galerkin Isogeometric analysis with Galerkin IGA-C Isogeometric analysis with collocation IGA-GL Isogeometric analysis with Galerkin by fitting load function IGA-CL Isogeometric analysis with collocation by fitting load function SC Superconvergent CG Cauchy-Galerkin ASC Alternating superconvergent CSC Clustered superconvergent