NURBS-Enhanced Meshfree Method with an Integration Subtraction Technique for Complex Topology

: In this paper, we present an integration subtraction technique to model holes interactively in a predesigned domain for adaptive problems. This technique involves two approaches, the normal subtraction technique and the moving subtraction technique. In the basic normal subtraction technique, the predesigned domain can be meshed using any methods as an initial integration background cell for meshfree analysis. Holes are described using closed non-uniform rational B-spline (NURBS) curves to preserve the exact computer-aided design (CAD) geometry and are meshed alone using the homotopic method, so they can easily be subtracted from the predesigned domain with no refinement. On the other hand, when the hole size is varying, the moving subtraction technique, in which only the changing part between the new and old boundaries needs to be integrated and subtracted, is more efficient. Compared with the standard radial point interpolation method (RPIM) and finite element method (FEM) in three linear elastic examples with different holes, the excellent accuracy and good efficiency of the proposed method are demonstrated, and its feasibility in complex topology problems is verified.


Introduction
As isogeometric analysis (IGA) has offered the possibility of integrating methods for analysis and computer-aided design (CAD) into a single, unified process, the time taken from design to analysis is greatly reduced [1]. IGA has been a topic of great interest in the area of recently developed technology in computational mechanics. When non-uniform rational B-spline (NURBS) has been used to describe complex geometry precisely and interactively in analysis, more physically accurate results and greater convenience in solving adaptive problems have been obtained.
In conventional isogeometric analysis, as proposed by Hughes, NURBS is used to integrate both of geometric and analysis information in one model [2], which leads to various advantages over the standard finite element method (FEM) [3][4][5][6][7]. However, the need for tensor product meshes has been a limitation in the analysis of arbitrary complex topology. To increase the efficiency of generating tensor product meshes for arbitrary complex topology, many methods have been developed, such as trimmed surfaces, patching operations, subdivision surfaces, and skeleton partitions [8][9][10][11]. In addition, other kinds of NURBS-enhanced FEM methods have been developed [12,13] at the same time. Furthermore, many other spline methods are promoted to overcome the drawback of tensor product requirement of NURBS, i.e., T-splines, PHT-splines, LHB-splines, etc. [14][15][16]. However, once the shape boundary or topology of an analysis domain is changed, the complicated mesh process needs to be performed again, and the global stiffness matrix needs to be partially corrected both in value and in size. These measurements are still somewhat tedious and time consuming in interactive design, topology optimization, crack growth, and other adaptive problems.
As the tensor product mesh has been a limitation in conventional isogeometric analysis, combinations of various meshfree methods and IGA have been developed by many researchers and applied successfully in many fields [17][18][19][20]. In meshfree methods, the domain of a problem is represented, ideally, only by a set of arbitrarily distributed nodes. The first meshfree method is invented by Gingold [21] and Lucy [22] for astrophysical problems without boundaries which is known as smooth particle hydrodynamics (SPH). To overcome the consistency and stability problem in SPH, reproducing kernel particle method (RKPM) and element free galerkin (EFG) are developed by Liu [23] and Belytschko [24]. As for the advantage of fast convergence ability and simplicity, radial basis function (RBF) is introduced to solve partial differential equations (PDEs) by Kansa in 1991 [25]. When a problem is solved in PDEs of integrated weak form in meshfree method, background cells or similar subdomains are always needed to achieve a consecutive integration. Triangular background cells, which can be easily created automatically for complicated domains, are sufficient for numerical integration operations in meshfree methods [26]. With the advantage of no special requirement for the mesh quality in meshfree methods, the mesh process in the integration region becomes easier. However, once the structure topology is changed, remeshing work for the new integration region or refinement along the inner boundary is needed. As new nodes and elements are added and the unnecessary ones are deleted in the process of refinement, the global stiffness matrix still needs to be partially corrected both in value and in size.
An attractive and competitive way to solve this problem is to couple a boundary element method (BEM) with IGA. A recent framework of the combination is proposed by Simpson [27], in which the NURBS functions were used to approximate the unknown fields in BEM. This IGABEM method circumvents the need to discretize the domain, as required by the IGAFEM [28]. However, the requirement for so-called fundamental solutions has been a limit for its application [29]. The scaled boundary finite element method (SBEM), introduced by Wolf and Song [30], shares the advantages of FEM and the BEM. Like in the FEM, no fundamental solution is required and like in the BEM, the spatial dimension is reduced by one, which would greatly decrease the freedom of degree. Furthermore, with the combination of IGA and SBEM, applications to problems in complex geometry have been made by Li [31] and Klinkel [32]. Moreover, some three dimension problems with singular [33] or nonsingular [34] boundary element method have also been investigated. The SBEM still requires a subdomain (having the same role as an element in the finite element method) to satisfy the scaling requirement (i.e., the whole boundary is directly visible from a point) in the discretization process [35], which is still a challenge work in the adaptive problems for complex two dimension domain and even three dimension domain.
The method we present here is the application of an integration subtraction technique to the meshfree method. The meshfree radial point interpolation method (RPIM) is chosen for its quick implementation. In the integration subtraction technique, the predesigned domain can be meshed as an initial integration background cell using one of many methods, while the awaiting designed hole is meshed using a homotopic method (such as the discretization idea in the scaled boundary method [32]) and subtracted from the initial integration background (short for RSUB), as illustrated in Figure  1. Unlike the scaling requirement in SBEM, there are no limitations for the shape of hole, which means convex and concave holes are all available, as the concave part would be cancelled out in the circle subtract integration process of the hole. Therefore, when the problems are expanded from two dimension to three dimension, the mesh for holes can still be easily generated, and no local refinement is needed with this homotopic meshing strategy. With this subtraction technique, the trimming curve used in CAD models to describe arbitrary topology can be directly applied in analysis procedures. As the integration of the new hole is calculated, the relevant part of the stiffness matrix can be directly subtracted from the global one to construct the new stiffness matrix with no change in size. Then, if the hole is changing, a more efficient moving subtraction technique (short for RMOV) is developed in which only the changing integration part between the old and new boundaries of the hole needs to be recalculated and subtracted. This will be discussed in detail later in Section 3.2.5. This paper is organized as follows. Section 2 provides the basic theory about the NURBS and RPIM meshfree methods. In Section 3, the NURBS-enhanced meshfree method and the integration subtraction technique are introduced. Section 4 provides numerical examples. In Section 5, final conclusions are presented.

NURBS
The ability of NURBS curves to efficiently describe arbitrary free-form shapes and all kinds of conic sections exactly has made NURBS the basis for all standard geometric exchange formats, e.g., Initial Graphics Exchange Specification (IGES) and Standard for the Exchange of Product Model Data (STEP) [36]. NURBS is the general form of the B-spline and is actually a linear combination of B-spline basis functions. B-Spline basic functions are defined over a parametric domain that is constructed by a so-called knot vector. A knot vector U is defined as: where all the i u are nondecreasing real numbers called knots, n + 1 is the number of B-spline functions, and p is the polynomial degree of the B-spline. In more detail, the knot vector U should be [37]: where i P are control points to determine the shape of the curve and , ( ) i p N u is the ith B-spline basic function of degree p, which is constructed by recurrence algorithms: The B-spline function An example of a B-spline basic function and control points to describe the arbitrary hole in Figure 1 with a polynomial degree of p = 3 is illustrated below in Figure 2. NURBS is the generalized form of a B-spline with the B-spline basic functions multiplied by their corresponding weights: Thus, with the tensor-production method, NURBS surfaces with two normalized parameters are defined as: where , ( )

RPIM Meshfree Method
Meshfree methods are developed to eliminate several key drawbacks of the finite element method [38]: 1. It is time consuming to generate a quality mesh in an arbitrary geometry with the desired accuracy. 2. It is difficult to construct approximations with an arbitrary order of continuity, making PDEs with higher-order differentiation or problems with discontinuities difficult to solve. 3. Performing h-or p-adaptive refinement is tedious. 4. The finite element method is ineffective in dealing with mesh entanglement-related difficulties (such as those in large deformation and fragment-impact problems).
As shown in Figure 3, the approximation function of unknowns in the meshfree method is constructed by scattered nodes, which are called field nodes. Thus, the information for approximation has been separated from the integration cells, which means that the mesh quality associated with approximation has been weakened, thereby giving us more freedom in discretizing the material domain. The RPIM approximation basic function is constructed by the radial basis function and point interpolation method on each field node within the supported domain. The known function values vector s U of inner nodes in the supported domain can be expressed as [26] where 0 R is the radial basis matrix, m P is the polynomial matrix, and a and b form the coefficient vector, which needs to be calculated. The polynomial term must satisfy the following conditional positive-definite constraint [39]: Combining the above equations and rewriting them in matrix form, we have Then, the function values of other unknown points (red triangles in Figure 3) in the supported domain can be approximated by where ( ) Φ x is the shape function of the RPIM, which represents each node's contribution to the unknown points in the supported domain.

Field nodes
The effects of different radial basis functions and shape parameters have been discussed in [26]. Here, in this presentation, the multiquadric radial basis function (MQ-RBF) is chosen for static elastic problems with shape parameters where c d is the average distance between field nodes and i r is the distance between the calculating point and the ith node. For a general 2D linear elastic problem in domain  , the static governing equation and boundary constraints can be described as: where T L is the gradient operator, which in 2D problems is generally defined as where σ is stress, b is body force, t and u are the traction vector and the displacement, respectively, t is the traction on the natural boundary t  , and u is the displacement on the essential boundary u  .
The general weak form for this problem is given by where D is the material constant matrix.
Assuming d u is the displacement of nodes in the support domain, h u and ε are the unknown displacement and strain of the integration point, respectively, which satisfy: After being integrated on the background meshes and assembled with all the field nodes, the weak form with no body force would become: It follows from the Kronecker  function property of the RPIM that boundary constraints can be directly applied to the field nodes. Then, Equation (16) can be rewritten as  KU F (17) which is the final discretization equation shaped by the meshfree RPIM. When the linear equation in Equation (17) is solved, the displacement U of all field nodes will be gained, and the stress σ can be recovered from:

NURBS-Enhanced Meshfree Method
In the NURBS-enhanced meshfree method, NURBS is used to describe the shape of the boundary and construct the background integration meshes as the foundation. As a result of the low requirement for integration mesh quality in the meshfree method, hybrid triangle elements can be chosen to accomplish an easy mesh for complex topology within the boundaries, as in Figure 4a. The integration elements and quadrature points are shown in Figure 4b. The integration points for the interior elements are generated with standard Gaussian quadrature rules. For the boundary elements, NURBS mapping quadrature rules for curved elements, such as the rules for high-order curved elements and NEFEM [13], are used. Even if the hybrid triangle mesh is relatively easy to implement, it is still time consuming to reconstruct the meshes and the stiffness matrix when the inner boundary is changed or a new hole is added. An alternative way to obtain the new stiffness matrix is to partially correct by refinement, but the measurement involved in deleting extra elements, repairing the gaps, reconstructing elements along the new boundary and redistributing nodes in the refinement region is relatively complex, as shown in Figure 5. Therefore, we recommend a subtraction technique in the integration process to avoid the complex measurement required for dealing with the nodes and elements in the refinement process, as illustrated in Figure 1.

Integration Subtraction Technique
In the integration subtraction technique, several steps are followed.

Meshing for a Predesigned Simply Connected Domain
Meshing for an arbitrary simply connected domain is relatively easy. Hybrid triangles are acceptable, while paver triangles with higher quality are still easy to apply, or a more convenient homotopic mesh method can be used, as illustrated in Figure 6. Homotopic meshing in 2D means the process of mapping the boundary curve 1 ( ) are closed NURBS curves, and the arrows in Figure 6c show the circumferential direction of 2 ( ) u C . Then, the physical locations of any points P within can be expressed as a function of the circumferential parameter u and radial parameter v by: collapses to a point, the homotopic mesh will be the same as the discretization idea in the scaled boundary method. What should be emphasized here is that the curve should be completely within the curve 1 ( ) u C . Otherwise, there will be elements (red slashed region in Figure 6d) outside the predesigned domain without field nodes to shape the approximation RPIM basic function. However, a point style 2 ( ) u C does not always exist for a nonconvex domain to ensure that all the elements are located inside the predesigned domain. However, there are still many other ways to accomplish the mesh for an arbitrary simply connected predesigned topology.

Deploying the Field Nodes
In conventional FEM or meshfree methods, field nodes are only deployed in the material domain within the boundaries. In the proposed method, field nodes are deployed all through the predesigned domain, including the region with holes. This is because the integration of the hole mesh requires the support of field nodes inside the hole.
To deploy the field nodes conveniently, nodes are first settled within a uniform interval in a region larger than the predesigned one. Then, the nodes outside the predesigned domain (blue in Figure 7) are removed using Gauss' flux theorem, the 2D form of which is where R is the Euclidean distance between the field node and the point on the boundary. Zero is just an ideal number to reach in the numerical solution, therefore we define any node with a result less than  to be outside. Finally, the remaining nodes are deployed along the boundary of the predesigned domain, as shown in Figure 7, in which purple labels of "12" means the fixed boundary condition in X and Y directions.

Meshing Homotopically for the Holes
Generally, the point-style homotopic mesh is available for holes of most shapes, unless they are too close to the boundary of the predesigned domain and are nonconvex, which would generate integration meshes outside the predesigned domain. With a point-style homotopic mesh, the integration of a hole can be executed in the parameter domain.
In the point-style homotopic mesh, the hole boundary 1 ( ) u C and reference point 2 ( ) u C are defined as: Then, the hole can be easily meshed in the parameter domain with n numbers in the circumferential direction u and m numbers in the radial direction v .

Numerical Integration and Normal Subtraction Technique
The integration on the physical domain for any function   , F x y can be transferred to the parameter domain in the form of Equation (23) is numerically integrated by standard Gaussian quadrature methods at each element. As the displacement on the Gaussian quadrature point can be expressed by the displacement of local supported field nodes in the RPIM interpolation, the stiffness matrix for a predesigned domain and h holes will be constructed on all the field nodes within the predesigned domain, as illustrated in Figure 7.
In a general 2D linear elastic problem, the stiffness of the jth element in the ith hole can be integrated over the parameter domain where g n is the total number of Gaussian quadrature points in an element, k w is the Gaussian integral weight of the kth point, and k F is the stiffness function expressed here as where k B is the strain matrix and D is the material constant matrix. In addition, k J in Equation (25) is the area integral Jacobian matrix where ( ) j du and ( ) j dv are the circumferential and radial lengths, respectively, of the jth element in the parameter domain. As a hole can be easily meshed in the parameter domain with n numbers in the circumferential direction and m numbers in the radial direction, the stiffness matrix of the jth element can be obtained by summing Assuming that the stiffness matrix of the predesigned domain is 0 K , the total stiffness matrix can be constructed by the integration subtraction technique as follows.
Then, the displacement of all the field nodes in Figure 7 is calculated as in the above equation. Even if holes are added or changed, as in Figure 5, no extra field nodes need to be added or deleted.

Moving Subtraction Technique for a Varying Hole
A hole is said to be varying if the control points of the hole are moving and overlapping regions still exist between the new and the old holes. In this case, the moving subtraction technique is more efficient than the normal subtraction technique. However, if no overlapping regions exist, as in faraway translation, implementing the deletion of old holes and the addition of new holes using previous normal subtraction techniques is more practical.
In the moving subtraction technique, when a hole varies from 1 ( ) , only the region between 1 ( ) u C and 2 ( ) u C needs to be integrated, as shown in Figure 8. This is mathematically equivalent to the direct operation of deleting hole 1  C and adding hole 2  C in Equation (29), but with less region to be integrated.
As the homotopic mesh is conveniently applied between 1 ( ) u C and 2 ( ) u C , as shown with the lines in Figure 8, the new global stiffness matrix can be easily obtained by adding the changing part of ( 2 1)  C C K .

Numerical Examples
In this section, three linear elasticity examples are presented to investigate the applicability of the proposed integration subtraction technique in a complex domain. The first example of the classic problem of a plate with a circular hole is used to discuss the convergence and efficiency of the presented method. The other two examples of bubbles with two holes are used to verify the feasibility of the complex topology of the normal subtraction technique and the moving subtraction technique.

Infinite Plate with a Circular Hole
The first example is an infinite plate with a circular hole in the center, subjected to a unidirectional tensile load of p = 1.0 in the x-direction, as shown in Figure 9a. Due to symmetry, only one quadrant of the plate is modeled in Figure 9b On the other hand, a plane stress case is discussed here, for which the exact displacement is: where the parameters  and k are given by Convergence and efficiency are first compared in the three tested methods, which are the presented integration subtraction method, the traditional RPIM meshfree method, and the bilinear isoparametric FEM, denoted hereafter as RSUB, RPIM and FEM, respectively. In the traditional meshfree RPIM, the mesh is the same as in the FEM, in which the nodes and elements are treated as field nodes and background cells, respectively.
The discretization of RSUB using integration meshes and the subtraction process is shown in Figure 10, while the discretization of RPIM and FEM is shown in Figure 11. A node interval of 0.2 h  is used for illustration.   T  q  e  I  I  I  I  disp  I  T  n  e  e  I  I (37) where the superscript e denotes the exact solutions, while the superscript q denotes the numerical solutions. The five different node intervals h = 0.25, 0.2, 0.125, 0.1, and 0.0625 are investigated.
In Figure 12, the two meshfree methods, RPIM and RSUB, demonstrate higher accuracy than FEM in displacement and energy, with RSUB appearing to be the most accurate among the three tested methods. The better accuracy of RSUB with coarse mesh may be obtained by the duplicate integration for subtraction, which may be able to improve the result around the boundary of hole. When the node interval gets smaller, this advantage in precision of RSUB becomes less significant and converges to RPIM. In addition to having greater displacement accuracy, RPIM has higher convergence rates. Moreover, the presented RSUB method achieves extremely high accuracy on a rough mesh. However, with respect to energy, the convergence processes of RPIM and RSUB are not good, although their accuracy is better than that of FEM. It has been shown by Liu and Gu [40] that a purely MQ-RBF cannot always ensure to exactly reproduce a linear field function. This could be one of the major reasons for the poor h-convergence in using MQ-RBF for field variable interpolations. As a result of the shortcomings of the MQ-RBF, its energy efficiency is still not satisfactory. However, with regard to displacement, the meshfree methods, RPIM and RSUB, achieve better efficiency than FEM, as shown in Figure 13. While RSUB seems to be a little more time consuming than RPIM, the duplicate integration in the hole region for subtraction yields greater accuracy. Thus, RSUB is actually more efficient than RPIM.

Bubble with Arbitrary Holes
As the problem of a single circular hole dug in the predesigned domain of a rectangle has been verified above, we now propose another example of an arbitrary domain: a bubble with arbitrary holes. In FEM, this bubble is discretized with an edge length of 0.03 with paver quad4 meshes, as illustrated in Figure 15. The material properties and constraints are as follows: Young's modulus 3 10 E  , Poisson's ratio 0.3 v  , the left and right side of the bubble are fixed in X and Y directions (shown with purple labels "12" in Figure 15), and a single force of 1.0 is applied to the upper side. As for the presented method, the predesigned domain and holes are all discretized in a point-style homotopic mesh, with n numbers in the circumferential direction and m numbers in the radial direction, as marked in Figure 16. Field nodes are deployed at a uniform interval of 0.04 within the predesigned domain with the same nodes as FEM along the boundary, as shown in Figure 16c.  In the presented RSUB method, the displacement and stress of the example are solved in the stiffness matrix constructed on the field nodes, as in Figure 16c. Then, the displacement and stress in the y-direction are compared for FEM, RPIM and RSUB, as shown in Figure 17. In addition, the detailed results of the displacement of nodes on the upper side of the boundary are listed in Figure  18 at the cusp where the force is loaded.  Comparing the solutions of displacement and stress in Figures 17 and 18, it is obvious that the RSUB results are consistent with those of RPIM and FEM. However, little oscillation is found near the boundary of the right-side irregular hole in the stress results of RSUB. The reason would be that the difference between the displacement of field nodes inner and outer the hole is too great. While in the stress recovery process, both of inner and outer field nodes are needed. In addition, stress is only accurately solved at the Gaussian interpolation point, so oscillation may appear in the recovery result on the boundary of hole.

Bubble with Varied Holes
The bubble with varied holes is considered here to verify the feasibility of the moving integration subtraction technique. Unlike the previous bubble, the hole with a complex boundary on the right side is changed from a red hole to a blue one in Figure 19b. Nevertheless, the essential boundary and natural boundary with a single traction force remain the same. There is no need to adjust the field nodes. The discretization in FEM and the mesh between varied holes (shown with a shadow line) are illustrated in Figure 19, in which purple labels of "12" means the fixed boundary condition in X and Y directions. The moving subtraction technique, which has been discussed in Section 3.2.5, is called RMOV for brevity. The shadow line in Figure 19b shows the discretization in the circumferential direction. The material property and constraints, the parameters of the point-style homotopic meshes of the circle hole, and the outer boundary remain unchanged. Then, the results of displacement and stress in the y-direction are compared among the four methods of FEM, RPIM, RSUB, and RMOV in Figure  20. The detailed displacement in the y-direction near the force is plotted in Figure 21.  Comparing the solution of displacement and stress given by the four methods in Figures 20 and  21, it can be seen that the result of RMOV is consistent with those of the other three methods. Therefore, RMOV is more efficient in dealing with hole-varying problems than RSUB, as RMOV requires less integration than RSUB. On the other hand, RMOV causes less oscillation than RSUB near a hole boundary with a complex shape in the stress field. Therefore, RMOV is more suitable for addressing hole-varying problems.

Conclusions
A NURBS-enhanced meshfree method with an integration subtraction technique is developed to model holes interactively in a predesigned domain for adaptive problems. With this method, the integration domain with complex topology can be easily and exactly described by its boundaries. By combining this method with the homotopic mesh method, any holes can be directly dug out from the predesigned domain with no refinement. In addition, a more efficient moving subtraction method is expanded to speed up hole-varying problems, in which only the region between the new and old holes boundaries must be integrated. The major discoveries arising from the comparison among FEM, RPIM, RSUB, and RMOV are as follows: 1. Both RPIM and RSUB are more accurate than FEM and show higher convergence in displacement but do not have good enough energy convergence. 2. The proposed RSUB method achieves superior accuracy on a coarse mesh and a smoother stress result with a regular hole. 3. The duplicated integration of holes for subtraction yields improvements in precision, making RSUB more efficient than the original RPIM. 4. The results of both the proposed RSUB method and RMOV in displacement and stress are in good agreement with those of FEM and RPIM but cause little oscillation in the stress field near the boundary of the irregular hole, with RMOV causing less than RSUB.
In conclusion, the proposed integration subtraction technique not only provides a convenient mesh for complex topology but also demonstrates excellent accuracy and good efficiency. Unfortunately, its performance in energy convergence is not good, and little oscillation around an irregular hole appears in the method. Other integration methods, such as nodal integration, may improve the performance in the area of energy convergence and greatly decrease the computation cost. Moreover, a new stress recovery method can be developed to eliminate the oscillation. These are the two major aspects of planned future work to perfect the proposed method.