Meshfree analysis with the aid of NURBS boundary

We present a meshfree analysis framework to integrate the geometric exactness of non-uniform rational B-splines (NURBS)-based isogeometric analysis (IGA) (Hughes et al., Comput Methods Appl Mech Eng 194:4135–4195, 2005) with the flexibility of meshfree approximations. In the framework, only the NURBS boundary surface immediately available from CAD tools is used to describe the exact problem domain, and meshfree particles are inserted inside the boundary surface, in a flexible manner, for construction of the approximation for analysis. Nitche’s method is employed for imposing essential boundary conditions and the domain integration in the Galerkin formulation is performed based on variationally consistent integration (VCI) to recover integration exactness. The NURBS boundary surface from CAD serves as an aid in selecting particle distributions and as the integration net for the boundary integration required both for the Nitche’s method and the VCI. As shown in numerical studies, the VCI is essential for the solution accuracy of the method. Several benchmarks are tested to examine the effectiveness of the proposed framework and numerical results are compared with those obtained by the IGA.


Introduction
In conventional finite element (FE) analysis procedures, the geometry of the problem domain is first represented using computer-aided design (CAD) tools, the mesh is then generated from CAD data, and lastly the analysis is performed using FE methods.It has been found that the effort for converting CAD data to a quality mesh is greater than that for analysis.To minimize the conversion effort, Hughes et al. pioneered isogeometric analysis (IGA), where the non-uniform rational B-splines (NURBS), a standard basis function in CAD, is adopted to match the CAD geometry and the same NURBS bases are used for analysis [1].The distinct characteristics whereby the geometry exactly matches with the CAD data and stays unchanged during mesh refinement together with other desirable properties of the IGA have been intensively explored for advanced engineering and scientific computing [1][2][3][4][5][6][7].
However, in computer graphics and CAD tools, the geometric representation of a three-dimensional object is described merely by its boundary (the enclosed surface) using the NURBS surface.As to analysis, it is necessary to create a well defined approximation function over the entire problem domain, including the boundary and the interior (volume).Since multivariate NURBS basis functions are constructed by tensor product operations, a trivariate NURBS solid typically is augmented by adding the third parametric coordinate pointing inwards from the boundary surface, e.g., [8,9].This way of constructing a NURBS solid often introduces overly redundant control points and poses challenges for the hadaptivity.Furthermore, when the geometry is topologically complex, subdivision into patches is needed [9,10], due to the fact that a trivariate NURBS mesh is uniform on a rectangular domain in the parametric space.Many advanced B-spline-based techniques, such as T-splines [11,12], have been introduced to remedy the abovementioned issues in the NURBS-based IGA.The local hierarchical refinement technique of NURBS has been also introduced to enhance the flexibility of IGA [13][14][15][16][17]. Nonetheless, the construction of T-splines or the employment of local refinement still relies on a structured gird in the parametric domain and the h-adaptive refinement remains restricted.Thus, many efforts have been made to blend NURBS-based IGA with other numerical methods that are versatile in regard to the h-adaptivity, such as meshfree methods.
On the other hand, meshfree methods do not require structured nodal connectivity information, such as meshes, to construct the approximation function and are less sensitive to the nodal distribution, and thus highly versatile for discretization.In addition, the meshfree shape functions that are based on polynomial reproduction, e.g., moving least squares (MLS) [18] or reproducing kernel (RK) [19,20], can achieve essentially arbitrary order of continuity and completeness.However, due to the nature of "meshfree," the boundary description of a domain discretized by meshfree particles is less precise in comparison to a mesh-based approximation.A common means to describe the boundary of a group of meshfree particles is through linear connections between contiguously outermost particles, which is similar to the FE discretization.Consequently, meshfree methods usually do not preserve the geometry when the discretization is refined or coarsened.Furthermore, since a meshfree approximation is inherently rational with complex support overlapping, leading to quadrature inaccuracy, extra efforts need to be made when it is employed in the Galerkin formulation.Many attempts to resolve the quadrature inaccuracy require certain types of mesh, such as a background mesh conforming with support domains, a Voronoi diagram, a triangular tessellation, or a octree mesh [18][19][20][21][22][23][24].Inevitably, the need of a structured mesh for domain integration not only deviates from the spirit of meshfree but also complicates the implementation.
To take advantage of the accuracy of geometric representations in IGA and the flexibility and adaptivity in meshfree approximations, many attempts have been made to couple NURBS-based IGA with meshfree methods.Rosolen and Arroyo presented a method to couple the local maximum entropy meshfree approximation with NURBS based isogeometric analysis through the imposition of reproducing conditions with a maximum entropy optimization [25].In this method, a thin layer of NURBS volume elements adhering to the boundary is generated for isogeometric boundary description and the interior of the problem domain is discretized by local maximum entropy basis functions.Wang and Zhang introduced a consistently coupled isometricmeshfree method based upon the reproducing conditions that are enforced at newly derived reproducing points in the NURBS parametric domain to ensure consistency [26].Similar to [25], a thin layer of NURBS elements extruded from the boundary is introduced for boundary description whereas the interior region is discretized by RK meshfree functions in [26,27].Valizadeh et al. coupled the IGA with meshfree method by enforcing reproducing conditions in the physical domain [28].In this work, a narrow band of NURBS volume near the boundary is introduced in order to couple with the interior meshfree approximation.Triangular B-splines have also been introduced to achieve both smooth geometric representation and the flexibility of triangular tessellation [29], although the optimal convergence is not guaranteed due to lack of completeness in a triangular Bsplines based approximation.Jia et al. adopted the triangular B-splines as the kernel function in the reproducing kernel framework to recover the completeness of the triangular Bspline [30].The abovementioned coupling methods can be viewed as specific applications under the general framework of partition of unity [31,32].Although the NURBS functions or B-spline functions are popularly used for IGA, Simkins et al. [33] proposed using the Reproducing Kernel Element Method [34] for isogeometric analysis.The RKEM approximation was first constructed based on quasi-uniform meshes that are loosely discretized the domain and then transformed, through a least-squares projection, to smoothly interpolate the boundary of domain.
Compound B-splines based IGA methods have also been proposed to alleviate the complexity of constructing analysissuitable NURBS volume elements.In these methods, the approximation function for analysis is constructed by trimming or applying Boolean operations to a primitive background NURBS mesh [35,36].Kim and Young proposed the spline-based meshfree method, where the approximation function is defined based on a background Eulerian NURBS mesh and then trimmed by the boundary that is represented by spline surfaces and follows the material deformation [37].Although trimmed NURBS patches well represent the geometry, in the Galerkin based formulation the domain integration for the trimmed patches requires special treatments to ensure accuracy.Moreover, the adaptive refinement in this type of methods is performed on the background mesh and hence is inflexible.
The idea of the proposed meshfree framework is that only the boundary surface that is immediately available from CAD tools is used to describe the exact problem domain, without the need of construction of NURBS volume, and then either structured or unstructured meshfree particles can be inserted inside the boundary surface for construction of the approximation function for analysis.Evidently, it enjoys both exact geometry that is directly obtained from CAD and the flexibility in discretization.Nonetheless, numerical issues associated with meshfree approximation, such as the imposition of the essential boundary condition and the domain integration, need to be adequately treated to ensure the desirable properties of the proposed method.In the proposed methods, we adopt the Nitsche's methods for imposing the essential boundary condition [38] and employ domain integration schemes suited for achieving "isogeometric meshfree" under the framework of variationally consistent integration (VCI) [39].
The remainder of the paper is organized as follows.An introduction of NURBS functions to construct the boundary of a body is given in Sect. 2. The algorithm for inserting nodes into the domain is also illustrated in this section.Section 3 is to review the Reproducing Kernel Particle method including the construction of shape function, the imposition of the essential boundary conditions and how to perform the domain integration to achieve the variational consistency.In Sect.4, several numerical examples of the benchmark problems are conducted to demonstrate the performance of the proposed isogeometric mesh-free method.Finally, the summary and the conclusions of the proposed method will be given in Sect. 5.

NURBS-based geometric representation and meshfree discretization
NURBS have been commonly used for representing curves and surfaces in geometric modeling [40].In the proposed method, the NURBS entities will be used to describe the boundary of the problem domain and the meshfree discretization, a set of points inside the boundary, will be created for construction of the RK approximation function.We review in this section the basic properties of NURBS functions and sketch the procedure of isogeometric meshfree discretization.

B-Spline basis functions
A B-spline basis function of order p is defined by the knot vector in the parametric space written as 1 and n is the total number of B-spline basis functions that form the B-spline entity.If knots in a knot vector are equally spaced, they are called uniform; otherwise, they are called non-uniform.The ith B-spline basis function in one dimension can be defined using the Cox-de Boor recursion formula [40], starting with piecewise constants ( p = 0): The support of N i, p (ξ ) is compact with the interval [ξ i , ξ i+ p+1 ), which is termed as the knot span.Therefore, the support of N i, p (ξ ) contains p + 1 knot spans.B-spline functions of order p is C p−1 -continuous if the knots are distinct.If a knot has multiplicity k, the basis functions at this knot are C p−k -continuous.Similarly, the order of continuous derivative decreases ktimes for knots or control points with multiplicity of k.
A B-spline curve can be generated by the linear combination of the B-spline basis functions: The coefficients P i for each basis function are referred to as control points.Piecewise linear interpolation of the control points gives the control polygon.B-spline basis functions do not possess Kronecker delta properties and, therefore, the control points, in general, are not interpolated by B-spline curves.
Multi-dimensional B-spline basis functions, B α (ξ), can be constructed by using the tensor product of one dimensional ones.
Adopting the multi-index notation, here α = (α A rational B-spline entity is formed from a projective transformation of a piecewise B-spline entity in the onehigher-dimensional space.Therefore, rational B-splines can be used to exactly represent any conic section, such as circles and ellipses.The basis functions of the rational B-spline curves can be expressed as where w i is the weight of the ith control point, which can be treated as the value of the control point in the (d + 1) dimension.A multivariate NURBS entity is analogously expressed as where is the rational B-spline basis in the parametric direction d.B-spline based functions possess the property of partition of unity.Moreover, an affine transformation is applied to a B-spline or NURBS entity by applying it to the control points.
Knots can be inserted without changing the geometry or parametric coordinates of a B-spline curve.This is one of the advantageous properties in B-spline applications.When a knot is inserted, one associated control point should be inserted simultaneously and the number of the basis functions is increased by one.Thus, knot insertion is analogous to hrefinement in FEM.Consider a B-spline curve, as in Eq. ( 2), defined by a knot vector, ξ = ξ 1 , ξ 2 , ...., ξ n+ p+1 .Let ξ ∈ ξ m , ξ m+1 ) be a new knot to be inserted.The corresponding new knot vector is ξ = ξ 1 , ξ 2 , .., ξ m , ξ, ξ m+1 , .., ξ n+ p+1 and the corresponding new basis function, Ni,p ξ , is generated by Eq. ( 1).The new control points, Pi , can be obtained by setting After solving a system of linear equations raised from Equation ( 7), the formula for computing all the new control points, Pi , can be shown as where It is noted that from Eqs. ( 8) and ( 9), only p new control points must be computed.Figure 1 shows an example of a B-spline curve in two dimensions, which is constructed from a uni- Fig. 1 B-spline curve, control polygons, and knot insertion form knot vector ξ = {0, 0, 0, 0, 1 4 , 2 4 , 3 4 , 1, 1, 1, 1}.Two levels of knot insertion are performed by inserting new knots at the middle of each knot span in the coarser knot vectors.The associated control points and control polygons obtained from ( 8) and ( 9) are shown in Fig. 1.It is worthwhile to note that control points are not interpolated by the B-spline curve; however, when more knots are inserted, the new control points get closer to the B-spline curve.

Isogeometric meshfree discretization
In computer graphics and CAD, a 3-D solid is typically described by its enclosed surface using NURBS.Without increasing the dimensionality, the NURBS-represented boundary cannot be used for numerical analysis.In contrast, in the proposed method, the NURBS boundary is adopted merely for exact boundary description and the analysissuitable approximation is constructed by point-based RK basis functions, which are located inside the NURBS boundary.The point distribution in the NURBS boundary is flexible since construction of RK shape functions does not require a structured grid.
The general procedure of creating the isogeometric meshfree discretization is described as follows: (1) Create points on the NURBS boundary in the physical domain (Fig. 2a).This can be simply achieved by adopting uniformly distributed nodes in the parametric space of the NURBS boundary.Note that this step is optional since the NURBS boundary is merely for boundary integration and determination of the physical domain in the proposed method, which will be further elucidated in the following section.Nonetheless, it helps to capture the variation near a curvy boundary.(2) Insert a uniform or non-uniform set of background points that covers the NURBS-represented domain (Fig. 2b).Note that advanced background-cell structures, such as the octree cells [24], can be used to efficiently select nodal distribution and identify neighboring information.The reader is referred to [24] for the details.(3) Deactivate background points that are outside or too close to the NURBS boundary (Fig. 2c).Many algorithms have been introduced to determine whether a point is inside/outside a NURBS-represented boundary [40][41][42].
In this study, an algorithm based on the closest-point projection is used and the detail is given in Appendix 1, from which the vector from a background point to the closest point on the NURBS boundary (denoted by V B N in equation ( 65)) can be obtained.For a point that has been determined to be inside (from (65)), the following distance criterion is adopted to determine whether the point is too near the boundary and should be excluded.
where h B is the minimum nodal distance in the background points.

Reproducing Kernel approximation
The reproducing kernel (RK) approximation was first introduced by Liu et al. [19] as the correction to the kernel approximation, which controls the error of discrete kernel approximation in the finite domain up to a desired order of accuracy.The discrete RK approximation can be also obtained by enforcing the polynomial reproducing condition and hence the consistency is achieved [19,20].The discrete RK approximation for a function is expressed as: where N P is the number of nodes (or basis functions), d I is the generalized coefficient associated with node I .The shape function I (x) consists of two parts: a (x−x I ) is the kernel function, which introduces the locality and defines the smoothness of the approximation function with a compact support measured by 'a'.C(x; x − x I ) is the correction function to be determined by satisfying the reproducing conditions.
The general correction function in m-dimension can be expressed by a set of nth order complete monomials with corresponding coefficients as: Here α = (α 1 , α 2 , . . ., α m ) denotes the m-dimensional index, with its length defined as |α| ≡ The coefficients b(x) are obtained by satisfying the nth order reproducing condition that states as: After applying (12), it is equivalent to and, substituting (13) into (15), the b(x) is determined as follows: where M(x) is the moment matrix of the kernel function a (x − x I ): The RK shape function is therefore obtained: 123 It has been shown [43][44][45] that the RK shape functions also satisfy the following gradient reproducing conditions: Unless stated otherwise, the C 2 continuous kernel function that is formed by the tensor product operation from the onedimensional cubic B-spline function is employed. and where a i is the support size in each dimension.
Remark 3.1 It is manifest that the RK approximation in Equation ( 18) can be constructed purely based on a set of points without a structured mesh, as depicted in Fig. 3. To be admissible for the approximation solution of a boundary value problem in ¯ ≡ ∪ ∂ , where is the problem domain and ∂ is its enclosed boundary, the RK approximation requires to be defined in the entire problem domain, that is, the moment matrix M(x) must be invertible for all x in ¯ .The necessary condition to have an invertible M(x) is to have sufficient kernels cover at x.The detailed discussions on the invertibility of M(x) can be found in [43,46].Unlike mesh based methods where the boundary of the mesh that is used to construct the approximation function also describes the boundary of the problem domain, the RK approximation typically does not possess a precise boundary description, unless pre-defined otherwise.Though not necessary the boundary of an RK approximation is commonly considered as straight

Piecewise linear approximation of boundary NURBS represented boundary
Fig. 3 RK discretization in two-dimension and NURBS represented boundary line segments that connect adjacent points on the set of predefined "boundary" points in 2D (or flat surface patches in 3D) (Fig. 3).The reason of selecting straight line segments or flat surface patches for boundary description is primarily for the simplicity of domain integration arising in the Galerkin method.

Statement of the linear elasticity problem
For purpose of demonstration, consider the elasticity problems as follows: and σ = C : ∇ s u, ∇ s u ≡ n is the outward unit normal vector on ∂ ; t is the surface traction on ∂ h , g is the prescribed displacement on ∂ g ; b is the body force.Note that the enclosed boundary ∂ of the domain is described by a NURBS or B-spline entity, ∂ = (ξ), as in Equations ( 2) and ( 6), which can be directly acquired from CAD tools.

Imposition of Dirichlet boundary condition
Due to lack of the Kronecker delta property, i.e., I (x J ) = δ I J , many techniques for the imposition of essential boundary conditions in the Galerkin meshfree formulation have been developed, for instance [38,47].Nitsche's method has been proved to be effective and variationally consistent among those techniques, and is particularly well-suited in the proposed isogeometric meshfree framework, which will become apparent as we proceed.
The weak form of ( 22)-( 24) with Nitsche's method for imposing the Dirichlet boundary condition states: find u ∈ U, for all v ∈ V , such that where the spaces of trial functions and test functions are U = u|u ∈ H 1 ( ) and V = v|v ∈ H 1 ( ) , respectively; (• , •) denotes the L 2 inner product of its two arguments; α is a large constant.Nitsche's method is often viewed as a consistent improvement of the penalty method.
The additional boundary integrals, the third term in the l.h.s. and the third term in the r.h.s in (25), are added to maintain the symmetry of the bilinear form and to recover consistency of the weak form.They can also be derived from the Lagrange multiplier method with substitution of the Lagrange multiplier by the traction on the Dirichlet boundary.The penalty terms are provided to ensure the coercivity of the bilinear form provided α is large enough.It has been proved in Nitsche's original work [48] that if α = β/ h is chosen, where β is a large enough constant and h is a discretization measure, the solution of ( 25) converges to the exact solution with the optimal rate in L 2 norm [38].The choice of parameter α depends both on the PDE and the Dirichlet boundary conditions.A tight estimation of β to ensure convergence can be obtained from the maximum eigenvalue of a generalized eigenvalue problem [49].

Variationally consistent integration (VCI)
The approximation solution of the weak form in Eq. ( 25) can be obtained by the Galerkin procedure and the optimal rate of convergence in L 2 norm can be achieved provided the numerical integration involved in ( 25) is sufficiently accurate.Nevertheless, since the mesh free approximation functions are rational with overlapping supports, conventional integration techniques cannot provide exact integration for these functions, leading to sub-optimal convergence in the Galerkin meshfree methods.Chen et al. [22] introduced the concept of fulfillment of integration constraints to recover the variational consistency and proposed a smoothed gradient over each conforming nodal representative domain to achieve the linear exactness.Similar ideas have been generalized in [39] where arbitrary order of integration exactness can be recovered for virtually any numerical integration techniques via modification of the test function to enforce variational consistency.As such, a desired order of exactness can be obtained in the proposed isogeometric meshfree method.The Galerkin formulation of (25) seeks to find Let the trial and test functions be approximated by: where I and ˆ I are proper basis functions, such as RK functions in (11), and d I and c I are associated coefficients.
To demonstrate that nth order exactness is achieved in the Galerkin formulation, consider the model problem in ( 22)-( 24) with an exact solution as nth order polynomials, denoted as u n : where x α and A α are nth order complete monomials and associated arbitrary coefficients in the multi-dimensional index notation, respectively.H (x) and A denote the vector forms of x α and A α , respectively.The body force and boundary conditions are consistent with the solution, i.e., Adopting the RK basis function as I in ( 27), the trial function can exactly recover the nth order solution u n , as which follows immediately from (14).Applying ( 30)-( 32), the weak form in (26) becomes Employing numerical integration in Eq. ( 34) and applying the approximation in (28) leads to the linear integration constraints to nth order of Galerkin exactness [39]: where , and , ∂ are quadrature versions of ( , ) and ( , ) ∂ .Considering an arbitrary nth order solution and arbitrary elasticity constants, the integration constraint can be further reduced to where H(x) is the vector form of {x α } |α|≤n−1 .The detailed derivation is provided in Appendix 1.
In order to achieve computational efficiency, an assumed gradient method is employed in (26) as in [39]; therefore, the integration constraint in Eq. ( 36) with a shifted center can be rewritten as: The assumed gradient is selected as: where R I (x − x I ) is a window function which can be taken in the simplest form as and bIi is the coefficients of H (x − x I ) to be determined.After substituting ( 38) into (37), the coefficients bIi is determined: with and Here r I i is the violation of the integration constraint.If a numerical integration scheme satisfies the integration constraint (35), no correction is needed and the method is the same as the unmodified one.

Total Lagrangian formulation
The total Lagrangian formulation is employed for hyperelasticity problems.In the following, X denotes the material coordinates in the initial configuration and x denotes the position in the current configuration.The domain of interest that initially occupies X is transformed to x with the enclosed boundary ∂ x in the current configuration.The variational form of the Lagrangian equation of motion is expressed as follows: where τ i j is the Cauchy stress.The Dirichlet boundary condition is imposed by the penalty method to avoid the linearization of traction on ∂ g x .Equation ( 43) can be transformed to where δu i, j = ∂u i ∂ X j , b 0 i and t 0 i are the body force and surface traction in the initial configuration, respectively.P i j is the first Piola-Kirchhoff stress: where F i k is the deformation gradient and S k j is the second Piola-Kirchhoff (P-K) stress.The second P-K stress can be obtained from the strain energy density function W : where E i j is the Green-Lagrange strain tensor defined as The linearized incremental form of the weak formulation can be written as where Similar to the linear case described in previous section, the integration constraint for the total Lagrangian formulation can be obtained by assuming an nth order solution for the test function in equation ( 44) as: where ∇ X ≡ ∂/∂ X and N is the unit outward normal vector in the initial configuration.Assuming that the first P-K stress P(u n ) can be described by polynomials as where D is an arbitrary coefficient vector, the integration constraints to be satisfied is The assumed gradient is employed as in (38) and the correction can be obtained following the procedure given in ( 37)- (42).Subsequently, the variational equation ( 44) can be solved by Newton-Raphson method.

Linear and Quadratic path tests for linear elasticity
The linear and quadratic path tests are conducted over a quarter of a circular ring as shown in Fig. 4a, where the boundary is exactly represented by NURBS.The inner and outer radii R 1 and R 2 of the quarter ring are 1 and 2, respectively.The boundary conditions are designed to be consistent with a linear and a quadratic solution for linear elasticity when the body force is absent, as given in Table 1.
In the tests, Neumann boundaries are prescribed on 1 and 3 whereas Dirichlet boundaries on 2 and 4 .The plane strain formulation is used with Young's modulus E = 10 5 and Poisson's ratio v = 0.3.The linear and quadratic basis RK approximations with corresponding order of VCI correction are adopted for linear and quadratic patch tests, respectively.Two different discretizations are considered while only one uniform background integration is used for the domain integration as depicted in Fig. 4b-d.
For demonstration purpose, the numerical results that are labeled by RK-VCI are obtained using 1x1 point Gauss integration within each background integration cell with VCI correction whereas those labeled by RK-Std.GI are obtained using standard Gauss integration methods with 2x2 (for linear) and 3x3 (for quadratic) integration points within each background cell.Table 2 gives the comparison of numerical results for RK-VCI and RK-Std.GI.It can be seen that numerical solutions with RK-VCI recover both linear   and quadratic exact solutions, up to the machine precision, regardless the nodal distribution.

A circular ring subjected to an internal pressure
A circular ring is subjected to a pressure on the internal surface while the outer surface is traction-free.The analytical solution of this problem for small deformation under plane strain assumption is given as follows [50]: A quarter domain of the ring is modeled with proper symmetric boundary conditions.The inner and outer radii R 1 and R 2 are 1 and 2, respectively.The internal pressure P = 10, the Young's modulus E = 10 5 , and Poisson's ratio v = 0.3 are adopted.The geometry of the problem domain is identical to example 4.1 as in Fig. 4a, where the boundaries 1 and 3 are Neumann boundaries and 2 and 4 are Dirichlet boundaries.For RK-approximations, normalized support sizes 2.0 and 2.8 are used for 1st order and 2nd order basis, respectively.The total support size of RK shape functions is the product of the nodal distance and the normalized support size.The exact geometry of the problem domain is represented by the NURBS functions.The numerical results obtained from RK-VCI are compared with those from RK-Std.GI and Isogeometric analysis (IGA).For RK discretization, the nodal distributions and the back ground integration cells are illustrated in Figs. 5 and 6, respectively.The integration cell size is uniform and proportional to associated nodal spacing.Within each integration cell, the 2 × 2 and 3 × 3 Gauss integration schemes are used for 1st order and 2nd order RK approximations, respectively.The refinement of IGA discretization is made by uniform knot insertion in the parametric space as shown in Fig. 7.The 2nd order NURBS is employed in IGA and 3 × 3 Gauss integration scheme within each knot cell is performed.Figures 8 and 9 show the convergence study in L 2 error norm for the displacement and stress fields.Figure 8 compares the convergence results from RK-VCI to the ones from RK-Std.The convergence rate, which is obtained by using linear regression of the convergence curve, is given in the legend of Fig. 8, indicating that with VCI corrections the numerical results from proposed method converge at the optimal rate in average.Figure 9 shows the comparison of convergence between RK-VCI and IGA.Considering different discretization refinement strategies for RK-VCI and IGA, the logarithm of total degrees of freedom (DOF) is used for convergence study in Fig. 9. Figure 10 shows the radial stress contours obtained by different methods using comparable DOFs.
As can be seen from Fig. 8, RK-VCI shows better accuracy than RK-Std.GI since RK-VCI has been corrected to be variationally consistent.Although the RK-VCI does not achieve the optimal convergence for the displacement error due to the refinement strategy, it yields nearly optimal converged results as the model is refined, while the numerical results of RK-Std.GI suffers from poor convergence.On the other hand, the comparison between RK-VCI and IGA is also made in terms of total DOFs.Under the same number of DOFs, the displacement solution of RK-VCI is less accu-123 rate than that of IGA whereas the stress solutions of the two methods are comparable.

Plate with a circular hole under uniaxial tension
A plate with a circular hole subjected to remote traction as shown in Fig. 11a is studied.The analytical solution under plane strain condition is given as [50]: A finite quarter model with symmetric boundary conditions is modeled as depicted in Fig. 11b, where the boundaries 1 , 3 and 4 are Neumann boundaries and 2 and 5 are Dirichlet boundaries.The plane strain formulation is adopted with Young's modulus E = 10 5 and Poisson's ratio v = 0.3.As shown in Fig. 12, four levels of refinement are employed for RK approximation in the convergence study.Normalized support sizes 2.0 and 2.6 are used for linear and quadratic RK approximations, respectively.The boundary of the discretization in Fig. 12 is exactly represented by NURBS.The integration cells are proportionally refined with the nodal distributions by the same manner as in Sect.4.2.The discretization of IGA is given in Fig. 13, where the refinement is implemented by uniform knot insertion.It is noted that the refined IGA mesh is not well balanced globally due to isoparametric mapping.Many knot insertion schemes have been introduced to generate balanced IGA mesh, e.g.[1,51], which, however, is not within the scope of this study.We only adopt the uniform knot insertion scheme when refining IGA mesh in the numerical examples.
The numerical results for the comparison between RK-VCI and RK-Std.GI are shown in Fig. 14.Similar to the previous example, the RK-VCI yields better accuracy than RK-Std.GI.The numerical result of RK-Std.GI has slower convergence than RK-VCI or even fails to converge in fine discretization.RK-VCI also shows better accuracy in comparison with IGA in this problem as shown in Fig. 15.The Fig. 12 Refinements of isogeometric meshfree method for a plate with a circular hole problem Fig. 13 Refinement of Isogeometric analysis for a plate with a circular hole problem radial stress contours of different methods under about the same number of DOFs are given in Fig. 16, in which the result from RK-VCI shows a smoother stress solution than those from RK-Std.GI and IGA.It is noticed that the IGA results (Fig. 16c) show a strong pattern that mirrors the IGA mesh shown in Fig. 13 (the third finest level).This indicates that the IGA shape function becomes distorted and produces high derivatives through isoparametric mapping, which eventually causes oscillation of stress solution.

A plate with an elliptical hole under bi-axial tension
This case is chosen to examine the proposed isogeometric meshfree method for problems with high stress concentration.The plate with an elliptical hole is subjected to remote bi-axial tension.This geometry and the loading conditions of this problem is given in Fig. 17a.The uniform surface traction σ 0 is applied in the infinite distance from the elliptical hole.The analytical solution under plane stress condition for this problem is given as [50]: A finite quarter model with symmetric boundary conditions is adopted and the dimension of this problem is illustrated in Fig. 17b, where Neumann boundaries are prescribed on 1 , 3 and 4 and Dirichlet ones are on 2 and 5 .The length of the major axis a = 2 and the minor axis b = 0.4, the Young's modulus E = 10 5 , and Poisson's ratio v = 0.3 are implemented in the numerical analyses.Four levels of discretization are used for RK approximation with normalized support sizes 2.0 and 2.8 for linear and quadratic RK approx-Fig.19 Refinements of Isogeometric analyses for a plate with an elliptical hole problem imations, respectively (Fig. 18).The comparable refinement of IGA discretization, which is implemented by uniform knot insertion, is also given in Fig. 19.
The solution comparison between RK-VCI and RK-Std.GI in terms of displacement and stress L 2 error norms is given in Fig. 20, where the RK-VCI solution shows better accuracy.However, due to high stress concentration, the RK-VCI with quasi-uniform refinement yields lower convergence rates in both displacement and stress solutions.Comparable low convergence rate is also observed in IGA as shown in Fig. 21.It is noted that the accuracy of IGA is marginally better than that of RK-VCI even though much denser mesh (or DOFs) appears near the apex of major axis in the IGA discretization.The stress contours are also given in Fig. 22, where the RK result shows slightly oscillatory stress distribution near the apex of the ellipse due to quasi-uniform nodal distribution.Nevertheless, the proposed method (RK-VCI) allows local refinement to be implemented in the region of interest in a flexible manner.To demonstrate, additional RK nodes, about 20 % more nodes, are inserted into the quasi-uniform mesh (Fig. 23a), in the vicinity of the apex of the major axis as shown in Fig. 23b.The numerical results listed in Table 3 shows that with the local refinement the displacement solu- tion accuracy has near an one-order improvement while stress solution accuracy has an about 25 % improvement.It is also noted that without VCI correction, RK std.GI, no improvement is observed when RK std.GI is used with local refine-ment.Figure 24 gives the stress contours for locally refined solutions.It is worthwhile to note that while the numerical oscillation persists when RK std.GI is used, it is noticeably reduced when RK VCI is used with local refinement.

Hyperelasticity problem: an infinite long cylinder subjected to an internal pressure
The performance of the RK-VCI for nonlinear nearly incompressible problems is examined in this section.The geometry and loading conditions are the same as in example 4.2 whereas Mooney-Rivlin hyperelastic model is adopted, in which the strain energy density function is: where I 1 , I 2 and I 3 are the first, second, and third invariants of Cauchy-Green deformation tensor and C 10 = 10, C 01 = 10, and k = 10000 are the material parameters.The analytical solution of the pressure-inner radius relationship under plane strain assumption is given in [41] as: where R 1 = 1 and R 2 = 2 are the inner and outer radii of the undeformed cylinder, respectively.r is the inner radius of the deformed cylinder and P is the internal pressure.A quarter model is used for numerical simulation with symmetric boundary conditions.For numerical analysis, linear RK-VCI and linear RK-Std.GI with 2 × 2 Gauss integration method are selected for this problem.Two different levels of discretization, with 86 and 296 DOFs, are considered.The numerical results of pressure-inner radius relationship are shown in Fig. 25, where the horizontal axis is the averaged inner radial displacement.With the coarse discretization, the pressure field has a error about 6.6 %, while with the finer discretization, the error reduces to 1.1 %.Both integration methods (RK Std.GI and RK VCI) provide quality results when finer discretization is used; nonetheless, RK-VCI yields better accuracy.

Conclusions
This paper presents a meshfree method that combines the benefit of exact geometry of IGA with the flexibility and adaptivity of meshfree approximations.Unlike coupled IGA and meshfree methods, only the boundary (surface) of a three-dimensional object that is immediately available from CAD is utilized to describe the exact problem domain, and meshfree nodes are inserted inside the boundary, either in a structured or unstructured manner, for analysis without the difficulty of creating an analysis-suitable NURBS volume mesh.In addition, the geometric exactness of the proposed method can be achieved regardless of the order of continuity and consistency of the approximation, i.e. even the linear order approximation in the present method can be used to perform isogeometric analysis.The Nitche's method is employed for imposition of essential boundary conditions.As such, identifying boundary nodes and recovering Kronecker Delta properties can be avoided and the required boundary integration can also be performed readily through the NURBS boundary description.The domain integration in the proposed method is performed based on the background grid, and more importantly, the variationally consistent integration (VCI) is employed to recover integration exactness and thus ensure solution accuracy.A generalized form of integration constraints for multi-dimensional elasticity is provided.As noted in the numerical studies, the variationally consistent domain integration is essential for the solution accuracy of the present method.With the correction of VCI, the proposed methods offer compatible solution accuracy to that of IGA with full integration while retaining the advantages of being a meshfree method.Although we demonstrate the proposed method using 2D problems, it can be extended to three-dimensional problems with ease.The method has potential to narrow the gap between design and analysis processes and accelerate the construction of computer models with complex geometry, for example, patient-specific models in the biomedical applications.order polynomials of the displacement: where C is the matrix form of elastic modulus tensor and ε(u n ) is the matrix form of strain tensor as given below.
In (67), H(x) is the vector form of {(x) α } |α|≤n−1 and Ā is the corresponding coefficients that are consistent with Eq. ( 29).Collectively, the Cauchy stress tensor can be rewritten as and consequently, the integration constraint, Eq. ( 35), can be expressed in a matrix form as:

Fig. 4
Fig. 4 Patch tests on a quarter of a circular ring.a geometry, b nodal distribution 1, c nodal distribution 2, and d background integration cells

Fig. 14 Fig. 15 Fig. 16 Fig. 17 A
Fig. 14 Convergence study of RK-VCI and RK-Std.GI for a plate with a circular hole problem.a displacement L 2 error norms and b stress L 2 error norms

Fig. 18 Fig. 20 Fig. 21 Fig. 22 Fig. 23
Fig.18 Refinements of nodal distributions for RK approximation for a plate with an elliptical hole problem

Fig. 24 Fig. 25
Fig. 24 Contour of stress σ xx for a locally refined plate with an elliptical hole problem.a 1st-order RK-Std.GI, b 1st-order RK-VCI, c 2nd-order RK-Std.GI, and d 2nd order RK-VCI

Table 1
Problem set up for linear and quadratic patch tests