A Meshless Local Petrov-Galerkin Shepard and Least-Squares Method Based on Duo Nodal Supports

The meshless Shepard and least-squares (MSLS) interpolation is a newly developed partition of unity(PU-) based method which removes the difficulties withmany othermeshlessmethods such as the lack of theKronecker delta property.TheMSLS interpolation is efficient to compute and retain compatibility for any basis function used. In this paper, we extend the MSLS interpolation to the local Petrov-Galerkin weak form and adopt the duo nodal support domain. In the new formulation, there is no need for employing singular weight functions as is required in the original MSLS and also no need for background mesh for integration. Numerical examples demonstrate the effectiveness and robustness of the present method.

Some currently popular meshless approximations are the moving least-squares (MLS) approximation, Shepard shape functions, partition of unity (PU), radial basis functions (RBF), reproducing kernel particle method (RKPM), point interpolation (PI), and Kriging interpolation (KI).Among them, the MLS approximation [30] is one of the most widely used approximations at present due to its global continuity, completeness, and robustness.However, the MLS approximation suffers from a number of problems that practically limit its application, namely, the high computational cost in obtaining the shape functions and their derivatives, difficulty in retaining accuracy with respect to nodal arrangement, and also the difficulty with which essential boundary conditions can be imposed due to the lack of the Kronecker delta property.Efforts have been made to address these problems by various means.Breitkopf et al. [31] developed the analytical forms for computing shape functions and diffuse derivatives of shape functions by assuming that some terms are constant and complete derivatives of shape functions.However, these formulations are dependent on the number of nodes and the formulation grows unwieldy when there are a large number of nodes in the support domain.Orthogonal basis functions pertaining to MLS for efficient and accurate computation of shape functions are investigated in [32,33].Singular weight functions are introduced by Kaljević and Saigal [34] to produce an interpolatory MLS approximation for the direct imposition of essential boundary conditions.Chen and Wang [35] developed a matrix transformation method for the imposition of boundary conditions; however, the formulation is complicated for implementation.The Shepard and least square (MSLS) interpolation developed by the authors [36] satisfactorily keeps the consistency of the Shepard shape functions up to the order of basis function and satisfies the delta property.However, singular weight functions have to be used to enforce the interpolating property of shape functions, which results in the loss of the smoothness of the interpolation as well as local oscillation.To eliminate this problem, an improved MSLS interpolation possessing the delta property without using singular weight functions as PU was proposed in [37].
Apart from the interpolation, the integration scheme of the weak form is also an important factor affecting the solution accuracy, which has been the problem for many meshless methods.Background cells have usually been used to integrate the weak form as has been the case with the EFGM, RKPM, and PIM [38].Due to the complexity of the shape functions, a large number of integration points are needed to avoid the underestimation of the weak form, which is computationally expensive but still not adequate to give accurate solutions.Furthermore when an irregular nodal arrangement is used, the background mesh has to be refined where nodes are densely distributed [39].To remove this difficulty, the nodal integration scheme which is free from the background cell was proposed in [40] and later used in [41,42].However, the performance of this scheme is unstable and also reduces the accuracy of the results.The MLPG-type meshless methods well solved this problem in a natural way by partitioning the local domain into a number of subdomains that may or may not overlap.In this way, integration of high order accuracy can be obtained for the global stiffness matrix [43] without background mesh.A similar but conceptually different approach was developed in [44] using the partition of unity quadrature scheme.However, the algorithms are complicated and more computationally demanding compared to the classical quadrature on subcells.Although many different approaches have been carried out for the weak form integration, finding a simple, efficient, and accurate integration scheme for meshless methods remains an open question.In this paper, the local weak form along with numerical integration over local subdomains is used to derive the discrete equations.
The content of the paper is outlined as follows.In Section 2, formulation of the improved MSLS interpolation is described in detail including the local approximation and support domain with dual definitions.The Kronecker delta property of the interpolation is also proved.In the Section 3, the discretised formulation of the present interpolation is derived using the local Petrov-Galerkin weak form.It is followed by numerical tests demonstrating the convergence characteristic and accuracy of the present interpolation in Section 4. Discussions are given at the end highlighting the features of the present method and suggestions on further studies.

Formulation of the MSLS Based on Duo Nodal Supports
In this section, the MSLS based on duo nodal supports is described in detail.We start the description of the formulation using a 2D problem domain of arbitrary shape as shown in Figure 1.The formulation is described in the context of elastostatics, with the fundamental field variable being a displacement.For an arbitrary node , its displacement vector in 2D is u  = (  , V  ) T , where   and V  are the nodal displacements in the  and  directions, respectively (the following formulation is derived only for   in the  direction but an identical process can be used for V  in the  direction).The interpolation at an arbitrary point x inside the domain in the  direction is expressed as where { 0  (x),  = 1, . . ., } is a set of shape functions that forms a partition of unity (PU); that is, ∑  =1  0  (x) ≡ 1,  is the node index, and  is the number of the nodes whose supports   include point x as shown in Figure 2.   (x) here is not the nodal displacement in the FEM or the "fictitious" nodal values in the EFGM but the local approximation of node  at x where the superscript  indicates local.Shepard shape functions used as PU are given by which is the same as in the original MSLS.The construction of the MSLS interpolation takes the following steps: firstly, construct the local approximation at each node and secondly apply the PU approximation to the local approximation to get the interpolation.The definition of the nodal support domain will be given in detail and the local approximations at a node will be described.

Local Approximation at a Node. The local approximation
(x) at an arbitrary node  is given by where   is the nodal displacement for the th node in support of node ,  is the total number of nodes in the local cover   of node  as shown in Figure 1, and    (x) is the modified least square shape function given by where    (x) are the modified least square shape functions of node  and are determined by the following equations: T is a polynomial basis, and  is the number of monomials in the basis.In the development of the MSLS interpolation, we use a bilinear basis throughout in 2D such that p T (x) = [1, , , ], and It can be seen from ( 4) that It has been proved in [36,37] that a singular weight function used as   (x) will enforce the interpolation, the equivalent equation ( 1) here, satisfying the delta property.A similar approach has been used by [45] to produce interpolatory MLS approximation.However, the use of singular weight function will bring some other problems such as the loss of smoothness of interpolation.

Duo Support Domain of a Node.
The support domain of a node is the area where a node exerts influence on the field variable.In this paper it is defined as a circle centered on that node although it may take other shapes such as a rectangle.
Here, two support domains are defined at each node, one is used in the construction of local approximation and the other in the PU approximation.In Figure 2, for example, node  has two support domains associated with it, namely,   with radius   and   with radius   .If a node, for example, node  in Figure 1 falls inside   , then node  will be used in constructing the local approximation at node .Similarly, if a point, point x in Figure 2, for example, is contained in   , then the local approximation   (x) will contribute to the PU approximation at x.For an arbitrary node, for example, node  in Figure 1, the size of   is defined by where  is a scale factor that ranges between 1.0 and 2.0,  is a coefficient such that  = 2 for a node lying on the boundary and  = 1 for all other nodes, and   is the distance between  and the fifth nearest neighbor node to .If there is a predefined triangular background mesh,   can be defined as the maximum distance between  and the nodes of triangles which are connected to the node .
For a node having its local support domain completely inside the domain, for example, the subdomain   of node  in Figure 2, the size of   is the same as   : For a node having its local support domain close to or intersecting the boundary, for example, node  shown in Figure 2, the definition of subdomain follows these steps.Firstly, find the nearest boundary node to  among the neighbor nodes which have been used as nodes in defining   , and secondly calculate the distance between the nearest boundary node and , denoted as   , and then the size of   is set by If there is an interior node where prescribed values needed to be applied, the procedure described above for setting the support domain of near boundary nodes can be repeated to that node using (10) by assigning this interior node as a boundary node.If we want all nodes to take nodal values at the nodes, the size of the   can be taken as the distance between the  and its nearest node for every node .The following quadratic spline function is used as the weight function over support domain in (2): where   = ‖x − x  ‖ is the distance between the point x and node  and x  is the coordinate of node .The aim of separately defining local domain and support domain is to produce MSLS interpolations having the delta property without using a singular weight, so that the difficulties associated with the use of singular weight function can be removed.Indeed this aim is achieved here if the domain for local approximation and domain for PU are defined by the method described above as will be proved later in Section 2.3.

2.3
. Delta Property at a Node.Consider a boundary node  to be applied with boundary conditions.If the support domain of the nodes is set according to ( 9) and ( 10), then the  will be the only node contained in   .Thus the MSLS interpolation equation ( 1) at x  becomes As there is only one node in the PU, then (2) becomes It is known by (7) that the local approximation u  (x  ) at node  satisfies Substituting ( 13) and ( 14) into (12) gives Hence, the present MSLS interpolation takes nodal value at boundary nodes and the essential boundary conditions or point load conditions can be directly imposed as in the FEM.Also, the present interpolation preserves the consistency up to the order of the basis function, which is a necessary requirement of accuracy.The proof is the same as has been presented in [36,37].

The Local Petrov-Galerkin Weak Form
Let  be the total number of nodes associated with the given point x; then (3) can be rewritten as where Φ 0 is the vector of Shepard shape function,  is a matrix comprising the modified least square point interpolation (LSPI) shape functions, and   (x) is the MSLS shape function.For domain Ω bounded by Γ (Figure 3), the equilibrium equations and boundary conditions of linear elasticity are given by where   is the stress tensor,   are the body forces,   are the unit normal to the domain, and Γ  and Γ  are the global boundaries with prescribed displacements and tractions, respectively.Similar to [36], the local polygonal subdomains are constructed for the purpose of simplifying the integration and the discrete equations.For example, based on Delaunay algorithm, a local polygon Ω  is constructed by using the  nodes in local cover   , as shown in Figure 3.A generalized local weak form of the equilibrium equation in ( 17) is written as where Ω  is the integration domain or subdomain for node  and V  is the test function.Using the divergence theorem in (18), we obtain the following local weak form: where   is the outward unit normal to the boundary Ω  .The boundary Ω  for the subdomain Ω  is usually composed of three parts: the internal boundary Γ  , the boundary Γ  , and Γ  , over which the essential and natural boundary conditions are specified.Substituting     =   in (19), the following is obtained: In order to simplify (20), we can deliberately select the threenode triangular FEM shape functions   , which correspond to the node  of the triangles constructing the polygonal subdomain Ω  , as test functions V  , such that they vanish over Γ  .Substituting shape functions   for V  in (20), we obtain the following local weak form: For a local polygonal subdomain Ω  located entirely in the global domain Ω, there is no intersection between Ω  and the global boundary Γ, and the integrals over Γ  and Γ  in ( 21) vanish.For a local polygonal subdomain Ω  near the boundary, the first item of ( 21) can also be omitted because of the properties of the test functions   .Substituting the MSLS approximation in (21) into the above equation leads to the following discretised system of linear equations: denoted as where D is the elasticity matrix: Equation ( 23) can be individually integrated over each triangle constructing the local subdomain Ω  , as shown in Figure 3.In the present work, seven Gaussian points are used in each triangle.

Numerical Examples
The proposed MSLSM interpolation has been coded in C++.
In this section, we show the performance of the present interpolation on a range of test problems.The results are compared with the exact solutions, the MLPG solutions, and the linear FEM solutions.The weight functions used in the MLPG for testing purpose are the Gaussian type weight functions given by where   is defined by (8) and   = 0.3  is used for all test examples.The scale factor  in ( 8) is set to be 1.5 and the linear bases are used in MLPG and MSLSM.To study the convergence behavior we define the following error norms in displacement and energy, respectively: where u is a vector collecting nodal displacement results where  is the infinitesimal strain tensor and  is the Cauchy stress tensor.The relative displacement error and energy error are calculated by where the superscripts num and exact refer to numerical solutions and exact (or reference) solutions, respectively.

A Constant Strain Patch Test.
A constant strain patch test [46] using three distributions of 7, 28, and 126 irregular nodes is shown in Figure 4. Young's modulus and Poission's ratio of the material are 1000 and 0.25, respectively.The thickness of the plate is taken as a unit following plane stress assumption.Since the exact solution is linear, a linear basis for the MSLS interpolation is able to represent this solution.The computational results in Table 1 show that the present MSLSM passes the patch tests exactly up to the double precision of the computer.

A Cantilever Beam.
A cantilever beam problem with dimensions of  = 8 m and  = 1 m, as shown in Figure 5, is tested first.The beam is subjected to parabolically distributed downward traction equivalent to a unit load  at the righthand end and is constrained at the left-hand end as shown in Figure 5.The elastic material properties used are  = 1 × 10 5 Pa and ] = 0.25, and the problem is solved under a plane strain assumption.We refer to the analytical solution of the problem given in [47].The convergence of the present   method is studied using four nodal arrangements with 50, 138, 402, and 965 nodes.The convergence rate is compared among FEM, MLPG, and the present MSLSM in Figures 6 and 7.It can be found that MSLSM shows a good accuracy and convergence rate.Figures 8 and 9 show the vertical displacement V and the normal stress   along the  = /2 indicating accuracy of results using irregular 138 nodes and the proposed formulation.
As has been highlighted in [36,37], the computational cost in obtaining the shape functions and its derivatives is much lower by the present LS interpolation than by the MLS approximation.This is shown in Table 2.It can also be observed from the table that the difference in computational   efficiency between the two interpolations gets bigger when the number of nodes increases.

4.
3. An Infinite Plate with a Circular Hole.The second example is an infinite plate with a circular hole of radius  = 1 m.The plate is subjected to a far field traction  = 1 Pa in the  direction.A finite portion of the plate is considered for analysis and, due to the symmetry of the problem, only a quarter of the portion requires modeling, as shown in Figure 10.The elastic material properties used are  = 3.0 × 10 7 Pa and ] = 0.3 and plane stress conditions are assumed.The stresses and displacements for this are given in an analytical solution in [47] as where  is the shear modulus and  is the Kolosov constant where  = (3 − ])/(1 − ]) for the plane strain assumption.Traction-prescribed boundary conditions consistent with the exact solution in ( 32) are applied at the top and right edges.Four distributions of 53, 188, 564, and 1012 nodes, which are shown in Figure 11, are employed for the convergence studies.Figures 12 and 13 show that the MSLS has a good convergence rate in the displacement and energy norm.
The convergence slope of the present method is similar to that of MLPG though the latter is seen to be more accurate.The displacement   obtained using MSLS and MLPG is shown in Figure 14.It is seen that the displacement values given by MSLS as well as MLPG are very close to the exact solution.

Conclusions
In this paper, we proposed a local weak form meshless Shepard and least-squares interpolation.The interpolation features the use of duo nodal supports for local approximation and global approximation, respectively.The present formulation satisfies the delta property at desired nodes without using singular weight functions.The shape functions constructed conform to the PU property.The local Petrov-Galerkin weak form is used so that there is no need for  background mesh.It has been tested that the present formulation is more efficient when compared to the currently widely used MLS approximation.For the 2D example of 965 nodes, the computing time of the present MSLS is close to half of the computing time of MLS.The accuracy of the present method is slightly lower than the MLPG but higher than the FEM when the same number of nodes and the same order of basis function, for example, linear or quadratic bases, are used (linear basis corresponding to triangular element in the FEM).The formulation here is derived for 2D analysis but is readily extendable to 3D and the essential ideas are the same.Further development of the present interpolation such as application to problems of changing geometry/moving boundaries, finite deformation, elastoplasticity, and threedimensional cracking is ongoing.

Figure 1 :Figure 2 :
Figure 1: The setting of the nodal subdomains in Ω.

Figure 4 :
Figure 4: Nodal arrangements for the constant strain patch tests.

Figure 6 :
Figure 6: Convergence of relative displacement error of the cantilever beam.

Figure 7 :
Figure 7: Convergence of relative energy error of the cantilever beam.

Figure 11 :
Figure 11: Nodal arrangements used for the infinite plate problem.

Figure 12 :Figure 13 :
Figure 12: Convergence of relative displacement error for the infinite plate problem.

Figure 14 :
Figure 14: Comparison of the horizontal displacement  along  = 0 by different methods.

Table 1 :
Results of the constant strain patch test.

Table 2 :
Comparison of computational time in obtaining the strain matrix (unit: second).