A Modification of the Moving Least-Squares Approximation in the Element-Free Galerkin Method

The element-free Galerkin (EFG) method is one of the widely used meshfree methods for solving partial differential equations. In the EFG method, shape functions are derived from a moving least-squares (MLS) approximation, which involves the inversion of a small matrix for every point of interest. To avoid the calculation of matrix inversion in the formulation of the shape functions, an improved MLS approximation is presented, where an orthogonal function system with a weight function is used. However, it can also lead to ill-conditioned or even singular system of equations. In this paper, aspects of the IMLS approximation are analyzed in detail. The reason why singularity problem occurs is studied. A novel technique based on matrix triangular process is proposed to solve this problem. It is shown that the EFG method with present technique is very effective in constructing shape functions. Numerical examples are illustrated to show the efficiency and accuracy of the proposed method. Although our study relies on monomial basis functions, it is more general than existing methods and can be extended to any basis functions.


Introduction
In recent years, the meshfree (meshless) method has been developed rapidly as a computational technique for solving partial differential equations.In the meshfree method, only a mesh of nodes and a boundary description are needed to develop the discrete equations.This makes it much flexible in solving problems with moving discontinuities (e.g., fracture of solids) and/or moving boundaries (e.g., shape optimization problems) that cannot be solved easily by conventional numerical methods, such as the finite element method (FEM), finite difference method (FDM), and boundary element method (BEM).
A group of meshfree methods have been proposed and developed, such as the diffuse element method (DEM) [1], the element-free Galerkin method (EFG) [2,3], the meshfree point interpolation method [4], the meshless method based on radial basis functions [5][6][7], the meshless local Petrov-Galerkin method [8], and the meshfree weak-strong (MWS) method [9].Many details about meshfree methods have been presented [10].Among these meshfree methods, the EFG method is a very promising method and is currently widely used in computational mechanics and other areas.
In the EFG method, moving least-squares (MLS) approximation is used to construct shape functions for each computational point.One disadvantage of MLS approximation is that a set of linear algebraic equations must be solved for every computational point of interest [11].If  basis functions were used to construct the shape functions, then an  ×  moment matrix, say A, must be inverted for every computational point when the discrete equations are assembled.With inappropriate nodes distribution and basis functions, A may be ill-conditioned or even be singular.Besides, in the postprocessing, A also must be inverted at each node when the displacements, strains, and stresses are computed.The computational cost associated with this is quite burdensome.Moreover, in order to retain the high accuracy of MLS approximation, the moment matrix A must be inverted accurately.If the inverse of A is obtained inaccurately, then the accuracy of EFG method is compromised.In order to ameliorate this shortcoming, an improved MLS (IMLS) approximation is described and studied in

The MLS Approximation with Orthogonal Functions
Consider a subdomain Ω x , the neighborhood of a point x and denoted as the domain of definition of the MLS approximation for the trial function at x, which is located in the problem domain Ω.To approximate the distribution of function  in Ω x , over a number of randomly located nodes {x  },  = 1, . . ., , the moving least-squares (MLS) approximation  ℎ (x) of , for all x ∈ Ω x , can be defined as where p  (x) = [ 1 (x),  2 (x), . . .,   (x)] are basis functions and a(x) is a vector containing the coefficient a  (x) ( = 1, 2, . . ., ), which are also functions of x.The basis functions can be chosen widely.In general, they can be chosen as the monomial bases, which are defined as follows: (i) linear basis (ii) quadratic basis The coefficient vector a(x) can be obtained at any point x by minimizing a weighted discrete  2 norm, which can be defined as where   (x) = (x−x  ) is the weight function associated with the node , with   (x) ≥ 0 for all x in the support domain of   (x), x  denotes the value of x at node , and  is the number of nodes in Ω x .In this paper, the cubic spline weight function is used where   = |x − x  | is the distance from node x  to point x and   is the size of the support domain for the weight function.
The matrices P and W in (4) are ]× (6) and û = [û where If A(x) is nonsingular, then from ( 7) we obtain Substituting a(x) into (1), the expression of the local approximation  ℎ (x) is thus where Φ(x) = { 1 (x) ⋅ ⋅ ⋅   (x)} = p  (x)A −1 (x)B(x) is the shape function and The derivative of these shape functions can also be obtained where A −1 , = −A −1 A , A −1 and the index following a comma is a spatial derivative.Since A(x) = PWP  , we should only solve the derivative of P and W to get the derivative of A(x).It should be noted that the above expression is a standard way to compute the derivative.Some novel and numerically accurate schemes for computations of derivatives of shape functions have also been studied; see [21].
As we mentioned above, (7) must be solved accurately for every point to retain the accuracy of the MLS approximation.In the meshfree methods, it may be time consuming and even cannot obtain desired accuracy when A(x) is ill-conditioned.To overcome these shortcomings, weighted orthogonal basis functions are proposed to derive the improved moving leastsquares approximation (IMLS).Using weighted orthogonal basis functions in the MLS approximation, the moment matrix A(x) in ( 7) is a diagonal matrix and the necessity for solving (7) can be eliminated.Moreover, the condition number of the matrix A(x) can be improved to some extent, but we cannot obtain a well-conditioned matrix A(x) when ill-conditioned matrix or singular problems occur; see the discussion in [16] and Section 3.
To diagonalize the moment matrix A(x) for arbitrary computational point x, the following orthogonality condition is imposed at any point x where a(x) is to be computed: where the function set { 1 (x),  2 (x), . . .,   (x)} can be termed a weighted orthogonal function set with weight functions {  } about points {x  }.
For the given arbitrary basis functions   (x) ( = 1, 2, . . ., ), the orthogonal basis functions   (x) can be obtained by using the Schmidt orthogonalization procedure as follows: where If the weighted orthogonal basis function set about {x  } is used, then (7) becomes where Ã(x) = QWQ  , B(x) = QW, and From ( 16), we can see that Ã(x) is a diagonal matrix.If (  ,   ) ̸ = 0 ( = 1, 2, . . ., ), then we can obtain the coefficients ã (x) directly: That is, where Equations ( 17) and ( 1) give the following expressions of the approximation function  ℎ (x) as where . ., ) are the shape functions.The derivative of these shape functions can also be obtained: where A , = −A Ã, A and the index following a comma is a spatial derivative.

A Modified Moving Least-Squares Approximation
In this section, more comprehensive analysis of the use of orthogonal basis functions in the MLS approximation is discussed, although some of them are described in [16].MLS approximation with orthogonal basis functions possesses some distinguished advantages, but there is also a possible singularity problem of the moment matrix Ã in its way.A typical example which shows the MLS approximation with orthogonal basis functions fails is presented.The main reason for the singularity problem is analyzed.Then a modified MLS approximation is proposed to overcome this problem.

Aspects of Use of Orthogonalization.
As discussed in [12,13,16,17], the aim of the orthogonalization process is to enhance the computational efficiency and remove inaccuracies in the MLS approximation by avoiding the inversion of an ill-conditioned matrix.Using the orthogonal basis functions, the moment matrix Ã is a diagonal matrix; thus the matrix inversion at each computational point is avoided, which can accelerate the computation especially for large problem.Moreover, the condition number of the diagonal moment matrix Ã can be improved greatly for a lot of problems.But this process fails when the original moment matrix A is a singular (ill-conditioned) matrix.That is to say the source of inaccuracy causing the ill-conditioning is not removed by this procedure.
It is shown in [12] that orthogonal basis functions   (x) can be obtained from Pascal basis   (x) by using the Schmidt orthogonalization.Thus, there is a linear relation between them; that is, each   (x) can be written as or in matrix form where G ∈ R × is a transform matrix mapping p(x) to q(x).
The Schmidt orthogonalization makes G a lower triangular matrix with all diagonal entries of unit value.Thus we obtain det(G) = 1, where det(♯) means the determinant of matrix ♯.Equation ( 24) makes the following equality hold: where P and Q are defined in ( 6) and ( 17), respectively.Then we have From ( 26), we can find the following.
(1) Ã is congruent with A, and these two matrices have the same inertia, which means these two matrices have same numbers of positive, negative, and zero eigenvalues.
(2) Solving ( 16) is equivalent to solving which means the matrix G can be viewed as a split preconditioner [22] for solving (7).Using a preconditioner, the condition number of the system matrix may be reduced.That is to say, the condition number of Ã = GAG  may be improved from matrix theory.This has been proved for a lot of problems in applications when A is ill-conditioned.(3) Since G is a lower triangular matrix with all diagonal entries of unit value, that is, det which means if A is ill-conditioned, then the determinant of A is close to zero.And we can get that det( Ã) is also close to zero.(4) If A and Ã are well defined, then the shape functions derived using either bases are identical; that is, Moreover, the derivatives of these shape functions are also identical; that is, From these aspects, we know that if A is not well defined (illconditioned or even singular), then we cannot obtain a wellconditioned diagonal moment matrix Ã.But the accuracy may be improved for many applications; this is because the transformation matrix G, which can be considered as a preconditioner, can reduce the condition number in some extent (see also Table 1).That is why many researchers believe that using the orthogonal basis functions in MLS approximation can improve the accuracy of solution [13][14][15].The following example also shows the above discussion.Consider a computational point x located at (1.5, 0.5), whose influence domain contains eight nodes; see Figure 1.These eight nodes are located in parallel lines in both  and  directions.The rectangle domain was chosen as the support domain of the computational point x.The radiuses of the support domain are   = 2.0 in  direction and   = 0.6 in  direction.These values are not of significance but mean that all eight nodes are located in the influence domain of x.
When quadratic basis and cubic weight functions are utilized to automatically construct the shape function, we find that the moment matrix is 01235 0.01852 0.00617 0.03189 0.00926 0.00617 0.01852 0.03189 0.00926 0.06019 0.01595 0.00926 0.00617 0.00926 0.00617 0.01595 0.00926 0.00617 0.03189 0.06019 0.01595 0.12139 0.03009 0.01595 0.00926 0.01595 0.00926 0.03009 0.01595 0.00926 0.00617 0.00926 0.00617 0.01595 0.00926 0.00617 Ã (x) = diag (0.01235, 0.00412, 0.00309, 0.00197, 0.00103, 0.00000) . (32) It is easy to see that Ã is also rank deficiency.This example shows that it is possible to form an ill-conditioned or singular linear system of equations, even if orthogonalization procedure is used.If A is singular or there is a zero on the diagonal of Ã, then no shape functions can be derived.If some nodes are perturbed, for example, node 8 is perturbed to node 8  (see Figure 1), then A becomes nonsingular but may be an ill-conditioned matrix.It should be noted that node 8  shown in Figure 1 is just one case; it can be perturbed in any direction indeed.In the perturbed case, using the orthogonalization process may be more effective than the traditional MLS approximation, since the condition number decreases and no matrix inversion is needed to compute.Table 1 presents the minimum values of the condition number of A and Ã, the minimum values of Ã , and the values of shape functions. 1 − 8 stand for the shape functions of node 1-node 8, respectively.Σ indicates the sum of these shape functions.

The Cause of Singularity Problem.
From above discussion, we know that the MLS approximation with orthogonal basis functions may also fail, although the number of nodes in the support domain of a computational point x is larger than the number of bases.The cause of singularity problem has been discussed in the original paper of the MLS approximation [11] and a recent paper [18], where the authors also presented the optimal radius of the support domain of a computational point for polynomial basis.In this subsection, we discuss the cause of singularity problem from matrix theory.
Close examination has shown that the problem of singularity arises from rank deficiency of matrix P. In fact, the matrix W = diag((x − x 1 ), . . ., (x − x  )) is a diagonal positive definite matrix for a given set of nodes in the support domain.If we further assume that the number of these nodes is larger than the number of monomial bases and P has full row rank, then we can obtain a nonsingular moment matrix A(x) from the matrix theory of view.The following proposition gives the nonsingular condition of the moment matrix A(x).The proof of this proposition is given in the appendix.

Proposition 1. The moment matrix A(x) is symmetric.
Assume that a set of nodes are located in the support domain of a computational node x; then the moment matrix A(x) is a nonsingular matrix provided that  ≥ , rank (P) = ; (33) that is, P has full row rank, where  is the number of nodes in the support domain and  is the number of basis functions.

Techniques to Avoid the Singularity.
To avoid a singular moment matrix, some techniques have been proposed, such as the perturbation method [16] and finding appropriate support domain for a computational point [18].Using a perturbation method, one could not determine the disturbance beforehand.If the disturbance is too small, the condition number will be very large; see Table 1 for a degradation example.If the disturbance is too large, the moment matrix may be well defined, but the nodes are not the original nodes; this will take large error for discrete linear systems.Finding appropriate support domain for a computational point is a good choice to avoid the singularity problem, but it is only efficient for regular nodes distribution and special basis functions, such as monomial bases.For irregular nodes distribution and general basis functions, this method is not very efficient.In this subsection, we discuss a new technique, which is similar to the matrix triangularization algorithm (MTA) [20].
As discussed above, the reason of the singularity problem of the moment matrix A(x) is the rank deficiency of P due to the improper enclosure of nodes and basis functions.For a given set of nodes, the basis functions are very important.The main idea of this technique is, therefore, to automatically find the monomials that are responsible for the rank deficiency for a given set of nodes.These monomials should be done specially to ensure a full rank moment matrix.But this technique should be an automatic procedure, and it has to be done without increasing too much computational cost and should not be too complex.
The details of the algorithm are as follows.
(1) Suppose that  monomials are chosen to obtain the basis functions and  nodes are selected in the support domain of an interpolation point x.Using ( 6), the matrix P can be obtained.It should be noted that the rows of P correspond to the monomials in the basis, and its columns associate with the nodes in the support domain.
(2) The matrix P is triangularized to determine the row rank .
(i) If  = , it means that P has full row rank.From Proposition 1, we know that the moment matrices A(x) and Ã(x) are generally well defined.In this case, the orthogonal basis functions q(x) can be obtained by using the Schmidt orthogonalization process.Thus, a diagonal moment matrix Ã(x) and its inversion A(x) are obtained.Therefore, shape functions can be easily got.In particular, the moment matrix A(x) may be ill-conditioned.In this case, using the Schmidt orthogonalization process is a good choice and shape functions can also be easily got (see discussion in Section 3.1).(ii) If  < , it means that there is a rank deficiency in P. Thus, the moment matrix is singular and we cannot obtain the shape functions.The reason of the rank deficiency is that some rows of P ( −  rows) are linearly related to other rows.Therefore, the  −  rows should be done with some special techniques.The procedure is as follows.
(a) In the row triangularization process, the permutations of rows are recorded.From the diagonal elements of the triangularized matrix, we can find out which rows (total  −  rows) are related to other rows.This means that these rows (which corresponded to the monomials) have no effect in forming shape functions.(b) We can do the Schmidt orthogonalization process for the remaining  basis functions.
Since basis functions are corresponding to the rows of P, we can set the corresponding  −  rows and  −  columns of A or Ã to be zero except the diagonal entry to be one and set the  −  rows of B or B to be zero.(c) Through the above operation, we can obtain a diagonal moment matrix Ã(x), which has full rank now.Thus the shape functions can be obtained easily.

Flowchart of the Technique.
The flowchart of the technique is briefly given as follows.
(1) Determine the support domain for point x to obtain an  ×  matrix P and an  ×  weight matrix W.  16)-( 18).(1) Construct the matrix P: using (6), the matrix P is constructed as follows:

Technique for an Eight-Node Approximation. In order to show how this technique works, an example of eightnode approximation, as shown in
It can be easy to find that rows of P correspond to the monomials in the bases, and its columns associate with the eight nodes in the support domain.(2) Row rank and row determination: using some permutation matrices, the matrix P is transformed to a row trapezoidal matrix to determine the row rank.The row trapezoidal matrix P is  It can be found that the sixth row of P is zero and the row rank is 5.It means that there is a rank deficiency in P. Therefore, there is a rank deficiency in the moment matrix A(x) from Proposition 1. From P we can also find that the sixth row is linearly related to other five rows.It should be noted that there is no row changed in the row triangularization process, which means that the sixth row of P is corresponding to the sixth row of P and the sixth row of P is also linearly related to other five rows.
(3) Since the sixth row is corresponding to the basis function  2 , the remaining basis functions {1, , ,  2 , } are used to construct orthogonal basis functions q(x) by Schmidt orthogonalization process.Now, the sixth row and the sixth column of Ã(x) are set to be zero except the diagonal entry to be one, and the sixth row of B() is set to be zero.In this case, we have the moment matrix Ã (x) = diag (0.01235, 0.00412, 0.00309, 0.00197, 0.00103, 1.00000) .
(4) Compute the MLS shape functions from ( 16)-( 18): after the above steps, the moment matrix Ã has full rank.The shape functions can be obtained easily.They are listed in Table 2.

Element-Free Galerkin Formulation
Consider the 2D problem of the deformation of a linear elastic medium from an undeformed domain Ω, enclosed by Γ: where  is the stress tensor corresponding to the displacement field u and b is a body force vector, with boundary conditions as follows: in which t and u are prescribed tractions and displacements, respectively, on the traction boundary Γ  and on the displacement boundary Γ  , and n is the unit outward normal matrix to the boundary Γ  .
Using the standard principle of minimum potential energy for (37)-(38), that is, to find u ∈ ( 1 ()) 2 , such that is stationary, where   () denotes the Sobolev space of order  and  is a function in the Sobolev space,  and  = D are strain and stress vectors, and D is the strain-stress matrix.MLS equation (10) or ( 21) is used to approximate the displacements in the Galerkin procedure.Then we can get Then substituting (40) into (39) leads to the following total potential energy in the matrix form, as and invoking Π = 0 results in the following linear system of u: in which the stiffness matrix K is built from 2 × 2 matrices K  and the right hand side vector f is built from the 2 × 1 vectors f  .These matrices and vectors are defined by where ] (for plane stress problem) . (44) In the sequel and unless mentioned otherwise, EFG method indicates the element-free Galerkin method with the original MLS approximation (10) and MEFG method stands for the element-free Galerkin method with the modified MLS approximation (21).These methods employ the moving least-squares (MLS) approximation or its modified form to approximate the trial functions.Another problem is that, in general, they do not have the property of nodal interpolants as in the FEM; that is,   (x  ) ̸ =   , where   (x  ) is the shape function corresponding to the node at x  , evaluated at a nodal point, x  , and   is the Kronecker data, unless the weight functions used in the MLS approximation are singular at nodal points [11].Therefore, essential boundary conditions should be imposed with special techniques [23,24].In this paper, we use the penalty method [23].

Numerical Experiments
In this section, two numerical examples are presented to show the efficiency of the EFG method with the modified

Cantilever Beam.
The first example is a cantilever beam of length  and height  subjected to traction at the free end (see Figure 2).The beam has a unit thickness and hence a plane stress problem is considered here.The closed-form solution is available for parabolic traction of force where the moment of inertia  of the beam is given  =  3 /12.
To evaluate the accuracy of the coupled method, the following error norms are used: where ‖♯‖ 2 is 2-norm of vector "♯", û and u are the approximation and exact solution of displacements, and σ and  are the approximation and exact value of stresses.
To show the effectiveness of EFG method with the modified MLS approximation studied in this paper, a typical regular node distribution (with 37 × 13 nodes) is shown in Figure 3.To evaluate the integrals, 15 × 8 background cells are used.For each background 4 × 4 Gauss points are employed.The rectangle domain is chosen as the support domain of the computational point x.We adopt two kinds of radius of the support domain.The first case is   = 2.2 in the  direction and   = 1.3 in the  direction.The second case is   =   = 3.0 in  and  directions, respectively.In the first case, singularity problem of the moment matrix A(x) occurs;  thus the EFG method fails, while the MEFG method can also obtain excellent agreement with the analytical solution; see Figures 4 and 5.In the second case, no singularity problem appears; that is, the moment matrix A(x) is always well defined.Thus both the EFG method and the MEFG method can obtain excellent agreement with the analytical solution; see also Figures 4 and 5.It should be also noted that, in the second case, the results of EFG method and MEFG method are almost the same, which confirms the theoretical analysis in Section 3.1.The convergence of the MEFG method and the accuracy of the MEFG method for irregular node distribution are also studied.Four different regular node distributions with 52 (13×4: 13 nodes in the  direction, 4 nodes in the  direction), 175 (25×7) nodes, 481 (37×13) nodes, and 784 (49×16) nodes and four different irregular node distributions are considered.Figure 3 also plots a typical irregular node distribution (with 481 nodes).For intergration for each problem, 5 × 4, 10 × 4, 15 × 8, 20 × 15 background cells are used.Here, the radius of the support domain is taken as   = 2.2 in the  direction and   = 2.3 in the  direction.The computational results are listed in Table 3.The convergence curves for displacements and stresses with different node distributions are plotted in Figure 6.In Figure 6, ℎ is the maximum size of node arrangement.From Table 3 and Figure 6, we can see that the MEFG method possesses convergence.Convergence for the regular node distribution problem is better than that for the irregular distribution problem.are also considered here.The computational results are listed in Table 4, where  1 is the relative error for the case   =   = 3.0 (which corresponds to the nonsingularity problem) and  2 is the relative error for the cases   = 3.0 and   = 1.5 (which corresponds to the singularity problem).The convergence curves are depicted in Figure 10.From Table 4 and Figure 10, we can also see that the MEFG method also possesses convergence.

Conclusion and Outlook
In this paper, aspects of MLS approximation with orthogonal basis functions, which are used to construct shape functions in the element-free Galerkin method, are analyzed in detail.
It is shown that this method has some advantages; for example, it avoids computing inversion of moment matrices; the condition number of moment matrix constructed by MLS approximation with orthogonal basis functions can be improved greatly for many problems.But it does not overcome the singularity problem of moment matrix, either.The reason why singularity problem occurs is studied.A novel technique based on matrix triangular process is proposed to solve this problem in detail.Our study relies on monomial basis functions, but it can be extended to any basis functions.Numerical examples are illustrated to show the effectiveness of the new approach.Although the proposed method has been investigated in the EFG method, it is also suitable for other meshfree methods based on MLS approximation.

Figure 1 :
Figure 1: An example with eight nodes and their coordinates.

( 2 )
Transform P to be a row trapezoidal matrix to get row rank  and record the exchanges of rows.(3) If  = , then go to the next step.If  < , then determine which  −  rows are related to other rows.The remaining  basis functions are used to construct orthogonal basis functions.Set these  −  rows and − columns of A or Ã to be zero except the diagonal entry to be one, and set the  −  rows of B or B to be zero.(4) Compute the shape functions from (

Figure 1 ,
is presented here in detail.These eight nodes are included in the support domain of an interpolation point x = [1.5, 0.5].The complete quadratic basis and the cubic weight function (5) are used.Thus  = 6 and  = 8.
1 , û2 , . . ., û ].Because the number of nodes  used in the MLS approximation is usually much larger than the number of unknown coefficients , the approximated function does not pass through the nodal values.Thus û ( = 1, 2, . . ., ) are the fictitious nodal values, not the nodal values in general.

Table 1 :
Degradation cases and perturbation of shape functions.

Table 2 :
Shape functions obtain by matrix triangular process.

Table 3 :
Convergence study for Example 4.1.