Advances in the Improved Element-Free Galerkin Methods: A Comprehensive Review

1  Introduction

As an important numerical computation method, the meshless (or meshfree) method operates solely on discrete points rather than predefined meshes. Consequently, it eliminates the need for remeshing when solving large deformation problems where mesh distortion typically occurs. This inherent advantage over the conventional finite element method (FEM) [1] has led to its widespread adoption across multiple scientific and engineering disciplines.

The development of the meshfree method originated in 1977 when Lucy pioneered the smoothed particle hydrodynamics (SPH) method [2], representing the first meshless computational method. Later in 1981, Lancaster and Salkauskas first proposed the moving least squares (MLS) approximation [3] for surface fitting applications. Building on this foundation, Nayroles et al. developed the diffuse element method (DEM) [4] in 1992, implementing MLS approximation for nodal approximation. Subsequently, in 1994, Belytschko et al. developed the element-free Galerkin (EFG) method [5], which also utilized MLS approximation with Lagrange multipliers for boundary condition enforcement. Since then, numerous other meshless methods have been proposed, including the reproducing kernel particle method (RKPM) [6], the finite point method [7], the partition of unity method [8], the radial basis function method [9], the boundary node method [10], the natural element method [11], the local Petrov-Galerkin method [12], the local boundary integral equation [13], the point interpolation method [14], the collocation meshless method [15], the hybrid boundary node method [16], the moving Kriging interpolation method [17], the meshless finite element method [18].

Among the various meshless methods, the EFG method proposed by Belytschko et al. [5] represents a fundamental and widely studied method in numerical computation. After the strong form of the governing differential equation is transformed into its weak integral form, and the shape function is introduced to derive the discrete equations, which are solved by using the computational tool of background integration mesh, ultimately yielding the numerical solution for the corresponding problem. The EFG method constructs shape functions using the MLS approximation [4], with error estimates analyzed by Li [19]. Noted for its foundation in least squares mathematical theory, the MLS approximation achieves high computational accuracy. A key advantage of the MLS method over the least squares method is its mobility and compact support. This allows for the use of a linear basis function to obtain stable and highly accurate numerical solutions. Currently, the EFG method has been applied to solve diverse problems, including transient heat transfer problem [20], hyperbolic problem [21], fracture mechanics [22,23], optimization design [24], dynamic boundary flow [25], magneto-electro-elastic plates [26], shrinkage analysis of round billets [27], free vibration of rotating nanobeams [28], wave propagation [29], Signorini-Tresca contact problem [30], topology optimization [31–33], free vibration [34], buckling [35] and nonlinear analyses [36,37]. Additionally, the coupled FEM-EFG method [38,39], the variational multiscale EFG method [40,41] and the enrichment of EFG method [42] have also been developed. However, the computation of the shape function in the EFG method demands significant computational resources, resulting in reduced efficiency. Additionally, it may encounter singular matrices during numerical implementation. These limitations necessitate further methodological improvements.

In 2003, Cheng and Chen [43] proposed the improved moving least squares (IMLS) approximation by replacing the basis function in MLS with orthogonal ones, resulting in higher computational efficiency. Later, in 2008, Zhang et al. studied the improved element-free Galerkin (IEFG) method based on the IMLS approximation and applied it to analyze 2D fracture problems. Subsequently, the IEFG method is successfully extended to solve various problems, including elasticity, potential problems [44,45], wave equation, heat conduction [46] and elastic dynamics problems. Cheng and Liew further applied the IEFG method to study the modified equal width (MEW) wave equation [47] and generalized Camassa and Holm equation [48]. Cheng et al. made significant contributions to advancing the IEFG method, applying it to problems such as the Schrödinger equation, viscoelasticity [49], fracture mechanics, topology optimization, elastoplasticity, and diffusional drug release [50] problems. Cheng et al. further integrated deep learning with the IEFG method for solving inverse potential problems [51]. Zhang et al. investigated the biological population problem using the IEFG method [52], while Li et al. conducted error analysis for the IEFG method [53] and later applied it to the elliptic problem [54], p-Laplacian problem [55], and Boussinesq equation [56]. Debbabi and BelhadjSalah selected the IEFG method for the thermo-elastic problem [57]. Mousavi et al. coupled it with the finite strip method to analyze plate vibration [58]. Guo et al. further extended its application to study vibration [59], nonlinear behavior [60], and flutter [61] in composite quadrilateral plates, in addition to the dynamics of cantilever plates [62]. The IEFG method offers the advantage of fast computational speed, effectively addressing the limitations of the EFG method, such as slow computation and the tendency to form singular matrices.

The shape functions in both EFG and IEFG methods lack interpolation properties, traditionally requiring the penalty method to impose the essential boundary conditions, which inevitably reduces computational efficiency. To address this limitation, Lancaster and Salkauskas proposed the interpolating MLS approximation [3], enabling direct imposition of Dirichlet boundary conditions without auxiliary methods. Ren et al. further simplified the derivation of the shape function and proposed the improved interpolating MLS approximation [63]. Wang and Sun analyzed its error estimates [64], while Li and Wang investigated its inherent instability [65]. By constructing shape function with this method, the improved boundary EFG method was proposed. Later, Ren et al. presented the interpolating EFG method for 2D elasticity [66], potential and heat conduction [67,68] problems. Liu and Cheng extended their application to study two-point boundary value, potential [69], heat conduction, elastoplasticity and large deformation [70] problems. Sun and Wang adopted the method for the regularized long wave equation [71]. Zhang and Peng applied it to the viscoelasticity problem [72]. Dehghan and Abbaszadeh solved the magnetohydrodynamics equation using the interpolating EFG method [73]. Shen et al. employed it for evolutionary variational inequalities [74] and the mixed complementarity problem [75]. Finally, Zhang and Li developed the variational multiscale interpolating EFG method [76,77] based on the improved interpolating MLS approximation [63]. Since the interpolating EFG method can directly impose boundary conditions like FEM and features a simpler formulation than the standard EFG method, it offers improved computational efficiency for solving scientific and engineering problems.

A limitation of the improved interpolating MLS method is that its weight function becomes singular at nodes, which inevitably introduces truncation errors. To overcome this issue, Wang et al. studied a novel interpolating MLS approximation [78] by replacing the singular weight function with the nonsingular alternative, subsequently developing the new interpolating EFG method [78]. Later, Sun et al. extended this method to solve elasticity [79] and elastoplasticity [80] problems. Dehghan et al. applied this method to diverse problems, including solutions for fractional diffusion-wave [81], inverse tempered fractional diffusion equations [82], and Oldroyd models [83], with additional contributions to decomposition variational multiscale interpolating EFG method [84]. Liu and Cheng further advanced the method for 2D large deformation analysis [85,86] and 3D elastoplasticity simulations [87]. Additionally, the interpolating boundary EFG method [88,89], variational multiscale interpolating EFG method [90], interpolating EFG scaled boundary method [91–93], coupled interpolating EFG method [94], regularized and adaptive orthogonal interpolating MLS approximation [95,96] were also developed by researchers. The interpolating meshless method based on non-singular weight functions offers a wide range of weight function options and overcomes the limitations associated with singular weight functions.

In 2022, Ma et al. combined the Hermite interpolation with MLS approximation to approximate the field variables, and then proposed the Hermite interpolation EFG method for solving the elasticity problem [97], functionally graded structures [98], piezoelectric materials [99], and functionally graded piezoelectric structures [100]. Its distinct advantage lies in providing superior accuracy for calculating function derivatives at any point within the solution domain. This stems from explicitly incorporating the normal derivatives of field functions during shape function construction, which is the distinctive feature absent in previous EFG improvements.

The EFG, IEFG, and three interpolation EFG methods all use scalar basis functions to establish shape functions. To reduce the number of undetermined coefficients in the shape function, Cheng and Li integrated the complex variable theory combined with the MLS approximation, developing the complex variable moving least squares (CVMLS) approximation. This innovation led to the complex variable element-free Galerkin (CVEFG) method [101] for solving elasticity problems. Typically, a linear basis function is selected for shape function construction. When a scalar basis is selected, the number of undetermined coefficients is three, which reduces to two for a vector basis function, thereby lowering the matrix order in the shape function computation.

The CVMLS approximation exhibits a minor limitation in that the mathematical and physical interpretation of its functional construction remains unclear. To address this issue, Bai et al. studied the improved complex variable moving least squares (ICVMLS) approximation [102] by replacing the basis function in the CVMLS approximation with a conjugate basis function. Subsequently, the improved complex variable element-free Galerkin (ICVEFG) method was successfully applied to solve a wide range of 2D problems including: elasticity [102], potential, advection-diffusion, Schrödinger, wave propagation, viscoelasticity [103], elastoplasticity, topology optimization [104], large deformation [105,106], as well as the bending problems of Kirchhoff plate and thin plate on elastic foundations [107]. Li and Tian further extended the ICVEFG method to large deformation [108–110] and fracture [111–113] problems. Li and Li conducted error estimation analyses for both the ICVMLS approximation and the ICVEFG method [114]. Additionally, Niu et al. integrated the conjugate basis function into the RKPM for solving the 2D elasticity problem [115]. However, the studies on the ICVEFG method in the aforementioned literature are confined to two-dimensional problems. To extend its applicability, Li and Li derived the 3D shape function calculation formula for CVMLS approximation [116], and Li subsequently developed the CVEFG method [117] for analyzing 3D partial differential equations. More recently, Yang et al. applied this method to solve 3D elastoplastic mechanics problems [118], thereby significantly expanding its scope of applicability.

The CVEFG and ICVEFG methods lack interpolation properties. Ren et al. addressed this limitation by developing the interpolating CVMLS approximation [119] through the incorporation of singular weight functions into the CVMLS framework. Building on this advancement, Deng and He successfully implemented the interpolating CVEFG method for temperature field analysis, and subsequently enhanced the approach by proposing an improved interpolating CVEFG method [120–124] by integrating the singular weight functions with ICVEFG approximations. This method addresses the inability of the ICVEFG method to directly apply boundary conditions and achieves improved computational efficiency.

Regardless of which of the aforementioned improved methods is employed, the influence domain of a point in 3D problems still constitutes a three-dimensional region. The number of points included within this domain is significantly larger than that in a 2D problem. Furthermore, shape functions must be computed for every point within the influence domain. As a result, the improvement in computational speed of the EFG method remains relatively constrained, particularly for 3D problems.

In 2017, Cheng et al. pioneered the integration of the dimension split method [125,126] with the ICVEFG method, leading to the development of the dimension splitting complex variable element-free Galerkin (DS-CVEFG) method [127]. This innovative approach demonstrated significant computational time savings compared to the IEFG method for solving 3D partial differential equations. The method operates by reducing the 3D problem to a series of 2D problems, which are then solved using the ICVEFG method with specially constructed 2D ICVMLS approximations. The resulting 2D discrete equations are further processed in the splitting direction through the simple finite difference method (FDM) discretization. Subsequent applications successfully addressed heat conduction [128], convection diffusion and elasticity problems with improved efficiency. Drawing on the successful research of DS-CVEFG method, researchers have progressively developed several related methods, including: the dimension splitting IEFG (DSIEFG) [129], the dimension splitting interpolating EFG methods employing singular weight functions [130], and their counterparts utilizing nonsingular weight functions [131–135], and the hybrid reproducing kernel particle method [136]. These studies collectively demonstrate that the dimension splitting meshless methods can subsequently enhance the calculation speed of traditional meshless methods for solving 3D problems.

A common feature of these dimension-splitting meshless methods is their use of the simple FDM along the splitting directions to discretize the equations derived from the meshless discretization of 2D problems. Considering that the interpolating meshless method and the FEM have broader applicability than the FDM and offer greater convenience in handling natural boundary conditions, Cheng et al. replaced the FDM-based splitting-direction discretization with the FEM while maintaining the IEFG method for 2D domain treatment, thereby developing the dimension coupling method [137,138]. Wang et al. innovatively applied the interpolating MLS method [78] for shape function construction in the splitting direction, leading to the development of dimension splitting interpolating MLS approximation [139–142]. This approach enabled the subsequent development of both the interpolating dimension splitting variational multiscale EFG [143,144] and generalized EFG [145,146] methods. These improved dimension splitting meshless methods can also improve the efficiency of traditional meshfree methods subsequently.

The dimension splitting meshless method can significantly enhance the computational efficiency of the IEFG method for solving 3D problems. This advantage is validated in the numerical examples section by solving the 3D Helmholtz equation. However, challenges arise when addressing nonlinear mechanical problems. Such complex problems currently require the traditional meshless methods for effective solutions. Therefore, the IEFG method is employed in the numerical examples to solve material and geometric nonlinear problems, demonstrating its superior computational speed over the standard EFG method.

This paper contributes a comprehensive review of various improved EFG methods, establishing learning foundations for new researchers while offering guidance for future advancements and research directions in EFG methods. The first section is dedicated to a literature review on enhanced EFG methods, analyzing their strengths and weaknesses. Since these advancements primarily focus on the shape function optimization, Section 2 provides detailed mathematical derivations of key shape function formulations. Section 3 introduces the IEFG method for solving 3D heat conduction, elastoplasticity and large deformation problems, as well as two types of dimension splitting meshless methods applied to the 3D Helmholtz equation. The final section presents conclusions and proposes four specific directions for future research.

2  Shape Functions

The main distinction between meshless methods and the FEM lies in their shape function construction. This section systematically presents the shape function formulation for the MLS, IMLS, interpolating MLS, and ICVMLS approximations.

2.1 The MLS Approximation

In 1981, Lancaster and Salkauskas proposed the MLS approximation [3], originally developed for surface fitting. After 1994, this technique became the foundation for constructing shape functions in the renowned EFG method.

The approximation uh(x) for the function u(x) is expressed as a polynomial, though the coefficients and degree of the polynomial are initially unknown. According to the research of Dolbow and Belytschko on the EFG method, even when using lower-order polynomial basis (such as linear basis functions), sufficiently accurate numerical solutions can be obtained by selecting appropriate weight functions [147].

When applying a polynomial to approximate an arbitrary function u(x), the local approximation function yields the following expression:

uh(x,x^)=∑i=1mqi(x^)ai(x)=qT(x^)a(x),x∈Ω,(1)

where x^ is a point within the local neighborhood of x; m is the number of basis functions used in approximation, and

a(x)=(a1(x),a2(x),⋯,am(x))T(2)

is unknown coefficient vector corresponding to the basis function. qT(x) denotes the basis function vector.

The coefficient vector a(x) can be obtained by minimizing the following functional

J=∑I=1nw(x−xI)[qT(xI)a(x)−uI]2=∑I=1nw(x−xI)[∑i=1mqi(xI)ai(x)−uI]2,(3)

where w(x−xI) (I=1,2,⋯,n) is weight function, xI refer to an arbitrary node in problem domain and the corresponding influence domain cover x.

Eq. (3) can be written as

J=(P(x)a(x)−u(x))TW(x)(P(x)a(x)−u(x))(4)

By minimizing the functional J, the expression for the unknown coefficient vector a(x) can be obtained. Let

∂J∂a=0,(5)

thus the matrix a(x) can be obtained as

a(x)=A−1(x)B(x)u(x),(6)

where

A(x)=[(q1,q1)(q1,q2)⋯(q1,qm)(q2,q1)(q2,q2)⋯(q2,qm)⋮⋮⋱⋮(qm,q1)(qm,q2)⋯(qm,qm)]=PT(x)W(x)P(x),(7)

B(x)=[w(x−x1)q(x1),w(x−x2)q(x2),⋯,w(x−xn)q(xn)]=PT(x)W(x),(8)

P(x)=[1q2(x1)⋯qm(x1)1q2(x2)⋯qm(x2)⋮⋮⋱⋮1q2(xn)⋯qm(xn)],(9)

u(x)=(u1,u2,⋯,un)T,(10)

W(x)=[w(x−x1)0⋯00w(x−x2)⋯0⋮⋮⋱⋮00⋯w(x−xn)].(11)

Typically, the cubic and quartic spline functions serve as the primary weight functions. The cubic spline function, defined as

w(d)={23−4d2+4d3d≤0.543−4d+4d2−43d30.5<d≤10d>1,(12)

while the quartic spline function is

w(d)={1−6d2+8d3−3d4d≤10d>1,(13)

where d is a function of x−xI, and

d=‖x−xI‖ρI,(14)

ρI=dmax⋅cI,(15)

the scales parameter dmax defines the size of the influence domain for a node, and cI denotes the distance between the point xI and its closest neighbor node.

Therefore, the approximate function can be written as

uh(x)=Φ(x)u=∑I=1nΦIuI,(16)

where

Φ(x)=(Φ1,Φ2,⋯,Φn)=q(x)A−1(x)B(x)(17)

is the shape function.

We arrange six nodes uniformly in [0, 1] and discretize the domain with five background integration cells, each having one Gauss point. Under this configuration, a cubic spline weight function is selected, and dmax is set to 2.0. The corresponding shape function and its derivative are shown in Figs. 1 and 2. It is shown that the shape function and its derivative are smooth at every node.

Figure 1: Schematic diagram of the MLS shape functions

Figure 2: Schematic diagram of the MLS derivative of shape functions

This shape function can achieve high computational accuracy due to its foundation in least squares mathematical theory. However, the inversion of matrix A(x) requires significant computational effort, leading to reduced efficiency. Additionally, it may encounter singular matrices during numerical implementation.

2.2 The IMLS Approximation

To simplify the high-computational cost are primarily focus of matrix A(x) inversion in the MLS approximation, Cheng and Chen [43] proposed the IMLS approximation. The formulation derivation proceeds as follows:

Applying Gram-Schmidt orthogonalization to the basis function qT (x) of MLS approximation yields:

pi=qi−∑k=1i−1(qi,pk)(pk,pk)pk,i=1,2,3,⋯,(18)

where

q=(qi)=(1,x1,x2,x3,x12,x22,x32,x1x2,x2x3,x3x1,⋯),(19)

and

(pi,pj)=0,i≠j(20)

Therefore, the matrix A in Eq. (7) can be simplified as the following new form

A~(x)=[(p1,p1)0⋯00(p2,p2)⋯0⋮⋮⋱⋮00⋯(pm,pm)].(21)

The resulting IMLS approximation shape function is therefore given by:

Φ∗(x)=(Φ1∗,Φ2∗,⋯,Φn∗)=pT (x)A~−1(x)B~(x),(22)

where

p(x)=(1,p2,p3,⋯,pm)T ,(23)

B~(x)=P~T(x)W(x),(24)

P~(x)=[1p2(x1)⋯pm(x1)1p2(x2)⋯pm(x2)⋮⋮⋱⋮1p2(xn)⋯pm(xn)].(25)

Evidently, the inverse of the modified matrix A~(x) in the IMLS approximation has a simpler form, which consequently accelerates the computation of shape function compared to the traditional MLS method.

2.3 The Interpolating MLS Approximation

Both the MLS and IMLS approximations lack interpolation properties, the penalty method is often used to impose the essential boundary condition in the corresponding meshless methods, which inevitably affects computational efficiency. To address this limitation, Ren et al. [63] proposed the interpolating MLS approximation based on the singular weight function.

When the zeroth-order basis function q = 1 is selected in the MLS method, thus the approximation function is

uh(x)=vT (x−xI)u(x)=b(x)u(x)A(x),(26)

where

vT (x−xI)=b(x)A(x)=(v(x−x1),v(x−x2),⋯,v(x−xn)),(27)

A(x)=∑I=1nw~(x−xI),(28)

b(x)=(w~(x−x1),w~(x−x2),⋯,w~(x−xn)),(29)

In the interpolating MLS approximation, the singular weight function can be selected as

w~(x−xI)={1|x−xI|β,|x−xI|≤ρI0,|x−xI|>ρI,(30)

or

w~(d)={1dβ(1−1d)2,d≤10,d>1,(31)

and β is a positive integer.

Because q1 is equal to 1 in the MLS approximation, after normalization processing of q1 at point x, we obtain:

q^1=q1||q1||x=[∑I=1nw~(x−xI)]−1/2,(32)

by orthogonalizing q2, q3, …, qm against q^1 at point x, we obtain

q^i=qi−(qi,q^1)(q^1,q^1)q^1=qi(x)−∑I=1nv(x−xI)qi(xI),i=2,3,4,⋯,(33)

The resulting interpolating MLS approximation shape function is therefore given by:

Φ^(x)=(Φ^1,Φ^2,⋯,Φ^n)=vT (x)+q⌢T (x)A⌢−1(x)B⌢(x),(34)

where

q⌢(x)=(q^2,q^3,⋯,q^m)T ,(35)

A⌢(x)=P⌢T(x)W~(x)P⌢(x),(36)

B⌢(x)=P⌢T(x)W~(x),(37)

P⌢(x)=[q^2(x1)q^3(x1)⋯q^m(x1)q^2(x2)q^3(x2)⋯q^m(x2)⋮⋮⋱⋮q^2(xn)q^3(xn)⋯q^m(xn)],(38)

W~(x)=[w~(x−x1)0⋯00w~(x−x2)⋯0⋮⋮⋱⋮00⋯w~(x−xn)].(39)

Therefore, the approximation function yields the following expression:

uh(x)=Φ^(x)u(x)=vT (x)u(x)+q⌢T (x)a⌢(x),(40)

where

a⌢(x)=A⌢−1(x)B⌢(x)u(x).(41)

The interpolating MLS approximation satisfies the Kronecker delta property, endowing its shape functions with interpolation characteristics that permit direct imposition of essential boundary conditions. Moreover, the unknown coefficient vector a⌢(x) in the shape functions contains one fewer element than in conventional MLS approximation, thereby reducing computational cost associated with shape function matrix operation.

2.4 The Interpolating MLS Approximation with Nonsingular Weight Function

A drawback of the interpolating MLS approximation proposed by Ren et al. [63] is weight function singularity in constructing shape function, leading to the inevitable truncation errors. To overcome this limitation, Wang et al. studied a new interpolating MLS approximation [78] by replacing the singular weight function with the nonsingular alternative.

Inspired from Eq. (34), we can derive a new set of basis function such that the approximation function satisfies the interpolation property. Let

q~i(x^)=qi(x^)−∑I=1nv~(x,xI)qi(xI),i=1,2,3,⋯,m,(42)

where v~(x−xI) must be satisfied as

v~(xI,xJ)=δIJ,(43)

and

∑I=1nv~(x,xI)=1.(44)

Therefore, the form of v~(x,xI) can be selected as

v~(x,xI)=ζ(x,xI)∑J=1nζ(x,xJ),(45)

where

ζ(x,xI)=∏I≠J‖x−xJ‖2‖xI−xJ‖2.(46)

From Eq. (42), we derive the same transformation for u(x^)

u~(x^)=u(x^)−∑I=1nv~(x,xI)u(xI),(47)

the corresponding local approximation function is

u~h(x,x^)=∑i=1mq⌣i(x^)a⌣i(x).(48)

From Eqs. (42) and (44), we have

q⌣1(x^)=1−∑I=1nv~(x,xI)=0.(49)

Thus Eq. (48) can be simplified as

u~h(x,x^)=∑i=2mq⌣i(x^)⋅a⌣i(x)=q⌣T(x^)a⌣(x),(50)

where

q⌣(x)=(q⌣2(x),q⌣3(x),⋯,q⌣m(x))T ,(51)

q⌣i(x)=qi(x)−∑I=1nv~(x,xI)qi(xI),(52)

a⌣(x)=(a⌣2(x),a⌣3(x),⋯,a⌣m(x))T.(53)

From Eqs. (47), (48) and (50), we have

uh(x,x^)=∑I=1nv~(x,xI)u(xI)+∑i=2mq⌣i(x^)⋅a⌣i(x).(54)

Thus the resulting shape function is therefore given by:

uh(x)=Φ~(x)u(x)=∑I=1nΦ~IuI,(55)

where

Φ~(x)=(Φ~1,Φ~2,⋯,Φ~n)=v~(x)+q⌣T (x)A⌣−1(x)B⌣(x),(56)

v~(x)=(v~(x,x1),v~(x,x2),…,v~(x,xn)).(57)

A⌣(x)=P⌣T(x)W(x)P⌣(x),(58)

B⌣(x)=P⌣T(x)W(x)(E−V(x)),(59)

   E is identity matrix, and

P⌣(x)=[q~2(x1)q~3(x1)⋯q~m(x1)q~2(x2)q~3(x2)⋯q~m(x2)⋮⋮⋱⋮q~2(xn)q~3(xn)⋯q~m(xn)],(60)

V(x)=[v~(x,x1)v~(x,x2)⋯v~(x,xn)v~(x,x1)v~(x,x2)⋯v~(x,xn)⋮⋮⋱⋮v~(x,x1)v~(x,x2)⋯v~(x,xn)].(61)

In contrast to the interpolating MLS approximation developed by Ren et al. [63], the shape function formulation proposed by Wang et al. [78] allows arbitrary weight function selection from standard MLS approximation. This method not only eliminates computational challenges caused by weight functions singularities but also avoids the associated truncation errors.

2.5 2D ICVMLS Approximation

The MLS, IMLS, and two type of interpolating MLS approximations are formulated for scalar function approximation. Bai et al. [102] developed the ICVMLS approximation to handle vector-valued function approximation.

For any function u(z), its approximation is

uh(z)=u1h(z)+iu2h(z)=q¯T(z)a(z)=∑i=1mq¯i(z)ai(z),z=x1+ix2∈Ω,(62)

q¯T(z) denotes the basis function vector, which is conjugate to qT(z),

qT(z)=(1,q2(z),⋯,qm(z)),(63)

and

aT(z)=(a1,a2,⋯,am).(64)

In general, the simple linear basis function is adopted for the ICVMLS approximation, which can be written as

q¯T(z)=(q¯1,q¯2)=(1,z¯)=(1,x1−ix2),(65)

Define the following functional

J=(Q¯(z)a(z)−u⌣(z))TW(z)(Q¯(z)a(z)−u⌣(z))¯,(66)

where

Q¯(z)=[11⋯1q¯2(z1)q¯2(z2)⋯q¯2(zn)]T ,(67)

u⌣(z)=(u1(z1)+i u2(z1), u1(z2)+i u2(z2),…, u1(zn)+i u2(zn))T,(68)

W(z)=[w(z−z1)0⋯00w(z−z2)⋯0⋮⋮⋱⋮00⋯w(z−zn)].(69)

Let

∂J∂a=0,(70)

we get

a(z)=A−1(z)B(z)u⌣(z),(71)

where

A(z)=QT(z)W(z)Q¯(z),(72)

B(z)=QT(z)W(z).(73)

Thus the approximation function can be obtained as

uh(z)=Φ(z)u⌣=∑I=1nΦI(z)u(zI),(74)

where the shape function is

Φ(z)=[Φ1,Φ2,⋯,Φn]=q¯T(z)A−1(z)B(z).(75)

Then Eq. (74) can be discretized as

u1h(z)=Re[Φ(z)u⌣],(76)

and

u2h(x)=Im[Φ(z)u⌣].(77)

If u(z) is a real valued function, uh(z)=u1h(z) in Eqs. (62) and (68) can be simplified as

u⌣(z)=(u1(z1), u1(z2),…, u1(zn))T.(78)

Since the ICVMLS method yields the approximation function, the expression in Eq. (77) tends to zero. This demonstrates that only the real component of shape function needs to be retained to approximate the real-valued functions.

The use of the ICVMLS approximation with linear basis function for 2D shape function result in only 2 undetermined coefficients, as opposed to the 3 required by the MLS approximation, thereby directly lowering the matrix order in the computation.

3  The Improved EFG Methods for Partial Differential Equations and Nonlinear Mechanics Problems

3.1 The IEFG Method for 3D Steady Heat Conduction Problem in Anisotropic Media

This section employs a steady heat conduction problem to validate both the effectiveness and computational efficiency of the IEFG method.

The governing equation is

k11∂2T∂x12+k22∂2T∂x22+k33∂2T∂x32+2k12∂2T∂x1∂x2+2k23∂2T∂x2∂x3+2k31∂2T∂x3∂x1=f(x),x=(x1,x2,x3)∈Ω,(79)

the Dirichlet and Neumann boundary conditions are

T(x)=T¯(x),x∈Γu,(80)

q(x)=k11∂T∂x1n1+k12∂T∂x2n1+k13∂T∂x3n1+k21∂T∂x1n2+k22∂T∂x2n2+k23∂T∂x3n2+k31∂T∂x1n3+k32∂T∂x2n3+k33∂T∂x3n3=q¯,x∈Γq(81)

where k11>0, k, k22>0 refer to thermal conductivity coefficients, and k12=k21, k23=k32, k31=k13, k11k22>k122, k22k33>k232, k33k11>k312. f(x) denotes the heat source density. Γ=Γu∪Γq, Γu∩Γq=∅, T¯ is known temperature and q¯ is given density of heat source. ni is component of the unit normal to boundary along i.

The functional corresponding to the governing heat conduction equation is constructed as follows:

Π=∫Ω12[k11(∂T∂x1)2+k22(∂T∂x2)2+k33(∂T∂x3)2]dΩ+∫Ω(k12∂T∂x1∂T∂x2+k23∂T∂x2∂T∂x3+k31∂T∂x3∂T∂x1)dΩ+∫ΩTfdΩ−∫ΓqTq¯dΓ+α2∫Γu(T−T¯)(T−T¯)dΓ(82)

where α is penalty factor.

The corresponding Galerkin weak form can be derived as

∫Ωδ⁢(L1⁢T)T⁢(L2⁢T)dΩ−∫ωδ⁢T⁢f⁡d⁢Ω−∫Γqδ⁢T⁢q―dΓ+α⁢∫Γuδ⁢T⁢TdΓ−α⁢∫Γuδ⁢T⁢T―dΓ=0,(83)

where

L1(⋅)=(k11∂∂x1,k22∂∂x2,k33∂∂x3,2k12∂∂x1,2k23∂∂x2,2k31∂∂x3)T(⋅),(84)

L2(⋅)=(k11∂∂x1,k22∂∂x2,k33∂∂x3,2k12∂∂x2,2k23∂∂x3,2k31∂∂x1)T(⋅).(85)

Assume that M discretized nodes xI (I=1,2,⋯,M) are distributed within the problem domain. We introduce the IMLS shape function to approximate T(x)

T(x)=Φ∗(x)T=∑I=1nΦI∗TI,(86)

T(x)=(T1,T2,⋯,Tn)T.(87)

Thus

L1T(x)=L1[Φ∗(x)T]=B1(x)T=∑I=1nB1I(x)TI,(88)

L2T(x)=L2[Φ∗(x)T]=B2(x)T=∑I=1nB2I(x)TI,(89)

where

B1(x)=(B11(x),B12(x),⋯,B1n(x)),(90)

B2(x)=(B21(x),B22(x),⋯,B2n(x)),(91)

B1I(x)=[k11ΦI,1∗(x)k22ΦI,2∗(x)k33ΦI,3∗(x)2k12ΦI,1∗(x)2k23ΦI,2∗(x)2k31ΦI,3∗(x)],(92)

B2I(x)=[k11ΦI,1∗(x)k22ΦI,2∗(x)k33ΦI,3∗(x)2k12ΦI,2∗(x)2k23ΦI,3∗(x)2k31ΦI,1∗(x)].(93)

By introducing the shape function, we obtain

∫Ωδ[B1(x)T]T[B2(x)T]dΩ+∫Ωδ[Φ∗(x)T]TfdΩ−∫Γqδ[Φ∗(x)T]Tq¯dΓ+α∫Γuδ[Φ∗(x)T]T[Φ∗(x)T]dΓ−α∫Γuδ[Φ∗(x)T]TT¯dΓ=0(94)

The final equation can be derived

KT=F,(95)

where

K=∫ΩB1T⁡(x)⁢B2⁡(x)dΩ+α⁢∫ΓuΦ∗T⁡(x)⁢Φ∗⁡(x)dΓ,(96)

F=−∫ΩΦ∗T⁡(x)⁢f⁡d⁢Ω+α⁢∫ΓuΦ∗T⁡(x)⁢T―dΓ+∫ΓqΦ∗T⁡(x)⁢q―dΓ.(97)

To verify the accuracy of the derived numerical solution formulation, we present the following numerical example, the governing equation is

∂2T∂x12+∂2T∂x22+∂2T∂x32+∂2T∂x1∂x2+∂2T∂x2∂x3+∂2T∂x3∂x1=0,(x1,x2,x3)∈[0,1]×[0,1]×[0,1],(98)

the boundary conditions are

T(0,x2,x3)=0.5x22−1.5x32+0.5x2x3+2,(99)

T(1,x2,x3)=0.5x22−1.5x32−0.5x2+x3+0.5x2x3+2.5,(100)