Advances in the Element-Free Galerkin Method: From Linear Solid Mechanics to Multi-Physics Applications and Hybrid Domain Coupling

1  Introduction

Recent advances in numerical methods have enabled the solution of a wide range of problems arising in science and engineering, leading to significant developments in computational mechanics and physics. Mesh-based methods such as the Finite Element Method (FEM) [1] and the Finite Volume Method (FVM) [2] remain the most widely employed techniques for the numerical approximation of governing equations, and their capabilities have been demonstrated in increasingly complex applications. These include linear and non-linear problems in solid mechanics [3], fluid dynamics [2,4], thermo-mechanics [3] and multiphysics [5,6], and even topology optimisation [7]. Such mesh-based numerical approaches rely on piecewise approximations of the field variables, which often pose challenges when capturing steep gradients [8–11] or ensuring the smoothness of derivatives [12–15]. Such a feature makes it difficult to handle problems involving large geometric deformations [16,17], moving boundaries [8,10,11], material discontinuities [11,18,19], or evolving interfaces that seldom align with the element or cell boundaries [11,15,20]. Although these issues can be overcome via adaptive re-meshing techniques, their implementation in complex geometries is computationally demanding and requires mapping of field variables between successive meshes. This process introduces additional costs and may degrade the accuracy and stability of the numerical solution [15,20,21]. To overcome these limitations, alternative strategies have been proposed and successfully applied. These include advanced mesh-handling techniques—such as arbitrary Lagrangian–Eulerian (ALE) formulations [22], sliding [23] and overlapping [24,25] meshes, immersed boundary methods [26], level set techniques [27,28], and also approaches that relax the dependence on mesh quality such as smoothed FEM variants [29] and meshfree methods [12]. Within the latter, the Element-Free Galerkin (EFG) method has emerged as a particularly powerful framework, offering flexibility in the discretisation of complex domains and the treatment of problems that challenge traditional mesh-based schemes. EFG methods were formally introduced by Belytschko et al. in 1994 [30] as an alternative framework for the numerical solution of the Galerkin weak form of the conservation equations, employing shape functions constructed from moving least squares (MLS) approximations. The method was used in the solution of benchmark elliptical-linear problems concerning elasticity, fracture mechanics and steady-state heat conduction. This pioneering work put into clear perspective several aspects concerning the accuracy and implementation of this then-emerging numerical approach, particularly the imposition of essential boundary conditions and the numerical integration requirements for improved accuracy. Since MLS approximations do not satisfy the Kronecker delta property, the use of Lagrange multipliers was proposed as the first strategy for enforcing essential boundary conditions. Furthermore, it was demonstrated that remarkable accuracy can be achieved through the use of high-order quadrature rules for numerical integration. Another major advantage highlighted was the absence of any post-processing requirement for the achievement of smooth physical fields obtained as derivatives of the primary variables, unlike in standard mesh-based techniques such as FEM. Subsequently, weighted orthogonal basis functions were introduced to construct MLS approximations at each quadrature point without the need to solve a system of linear equations [31]. This procedure was later termed improved MLS (IMLS), and the corresponding Galerkin weak formulation gave rise to the improved EFG (IEFG) method [32]. The first developments on EFG methods reflected an early recognition of their potential for practical adoption, mostly in elasticity problems. Within this framework, the intrinsic features of EFG enabled the design of simple and effective techniques to address classical challenges such as the capture of stress concentrations [30,31], material discontinuities [33], and the prevention of locking phenomena in nearly incompressible materials [30,34] and bending-dominated problems [31,35]. Although issues such as volumetric and shear locking may still arise in EFG formulations, these effects tend to be less severe compared to their FEM counterparts [30,31]. The growing dissemination, increasing recognition, and continuous refinement of EFG methods motivated works devoted to facilitating their numerical implementation, such as the development of dedicated libraries for the computation of MLS approximations [36] and also introductory contributions aimed at guiding researchers in programming the method [37]. Further efforts to improve the computational efficiency of EFG methods enabled the first implementations for three-dimensional problems [38], consolidating this mesh-less approach as a versatile framework for computational mechanics. As EFG methods have been applied to progressively larger and more complex problems, considerable and continuing efforts have been devoted to mitigating their computational challenges. The construction of shape functions for EFG methods involves the search for neighbouring nodes supporting each integration point in the problem domain, which is actually the most time-consuming stage in the computational implementation of EFG methods [39–42]. Moreover, both sparse matrix assembly and shape function construction become increasingly demanding as the number of supporting nodes per integration point grows [43]. This is mainly because achieving stable results and optimal convergence typically requires the use of high-order integration rules [44,45]. The approaches proposed to enhance computational performance in EFG methods include: the development of more efficient integration schemes such as post-processing techniques for nodal integrations based on Voronoi diagrams [46,47], virtual element decompositions [48,49], stabilised [50] and variationally consistent [51–53] formulations, reproducing kernel gradient smoothing frameworks [54–56] and also projection of EFG shape functions in FEM spaces [57]; improved implementation of sparse matrices assembly procedures and construction of basis functions [43]; computation of deferred correction vectors for the solution of non-linear problems as alternative to matrix reassembly [11]; application of reduced order modelling techniques [58–60]; more efficient neighbouring nodes searching procedures [40,43]; construction of approximations in complex variable spaces to reduce the dimensionality of basis functions in both MLS [61,62] and IMLS [63,64] approximations; and the implementation of techniques to reduce the number of nodes required to represent the problem domain, such as dimension splitting techniques for complex variable [65,66] and standard MLS [62,63,67], and overlapping discretisation approaches [19,41,42]. The ongoing development of procedures devote to improving efficiency, accuracy, and stability of EFG methods has enabled their extension to increasingly complex frameworks in computational mechanics, beyond linear elasticity and other standard elliptical problems. Some of these include nonlinear heat transfer with phase change involving transient conduction in fixed [68–71] and variable domains with moving boundaries [8,10,72], moving heat sources [41,42] and also orthotropic media [20], advection–diffusion [9,11,73,74] mechanisms, linear [18,75–77] and non-linear [16,78–80] thermo-mechanics, plasticity [16,81,82], elastoplasticity [83–85], large deformation for incompressible [86–89] and near-compressible [17,90–92] hyper-elasticity and plasticity [16,93], anisotropy [18,20,94–96], analysis of incompressible Stokes flow under two-level [97,98], generalised [67,99,100] and stabilised divergence-free [66,101] formulations, non-linear fluid-dynamics problems that include implementation of standard Characteristic Based Split techniques for transient analysis [102–104] within EFG frameworks, mixed u-p formulations [105,106], fluid-structure interaction [15], benchmark analyses at high Reynolds numbers [107,108], and free surface flow [109], coupled heat transfer and fluid flow [11,110–113], viscoelasticity [114] and viscoelastic flows [115], linear [116] and non-linear [112,117,118] porous media flow problems, other multiphysics scenarios such as magnetohydrodynamics for both moderate [60,119,120] and high Hartmann numbers [121–124] and phase field modelling in fracture mechanics [125–128], heat transfer [129–133], structural [133–137] and thermo-mechanical [137–140] topology optimisation, and even the solution of partial differential equations in modern physics [141–144] and biological [145–147] applications.

The studies reported so far put into an appropriate perspective how the various improvements and modifications introduced in EFG formulations enabled their widespread application beyond the linear elasticity framework within which they were originally developed. This proliferation has also been accelerated by the fact that most of the enhancements and adaptations implemented to extend EFG methods to transport phenomena and multiphysics problems have been largely based on its analogies with FEM, beyond the well-known differences concerning the construction of shape functions and the assembly of the algebraic system of equations. Consequently, most of the adaptations and improvements of EFG have been achieved through the incorporation of techniques already well established and extensively developed within the FEM framework. This has been achieved by retaining the distinctive advantages of EFG regarding (i) the possibility of more easily achieving higher-order approximations with continuous derivatives than FEM and other mesh-based techniques, and (ii) the greater flexibility in adding or removing nodes [12,19,42]. In contrast, achieving such continuity in FEM is challenging and typically requires either the use of highly complex finite elements with numerous nodal unknowns or post-processing to handle the discontinuous derivative fields generated by simpler elements [12,19,42,148]. These particular features have strongly encouraged the implementation of EFG methods for the solution of increasingly complex problems beyond standard linear elasticity and elliptic potential cases. The techniques inherited from FEM to address difficulties common to this well-established mesh-based framework include the implementation of dimension-splitting methods (DSM) to transform higher-dimensional heat transfer [62], fluid dynamics [67], advection-diffusion [149–151] or even wave propagation [63] problems into multiple lower-dimensional sub-problems sequentially coupled via finite differences; proper orthogonal decomposition (POD) to construct reduced-order models (ROM) with fewer degrees of freedom [58–60]; variational multiscale techniques to stabilise saddle-point problems arising in multiphysics from the coupling of fields such as velocity-pressure in linear [97,98] a non-linear [60,108,109,152] fluid dynamics, electric potential-current density-axial velocity for fully developed MHD flow in channels [119,123,124], temperature-pressure-velocity in coupled fluid flow heat transfer [110–112,153], and velocity-pressure coupling in porous media flow problems [116,118], and to suppress spurious oscillations under advection-dominated conditions in pure convection-diffusion [58,149,154–156], convection-diffusion-reaction [157–159], discontinuities capturing [160] and fluid dynamics [108,152]; reduced integration schemes to mitigate chequerboard-type instabilities in incompressible flow problems [15,106,107] and heat transfer coupling [11,113]; extended approximations with cover functions based on linear combinations of polynomials to satisfy the inf–sup condition in Stokes flow problems [67,99]; streamline-upwind Petrov–Galerkin (SUPG) stabilisation for advection-dominated regimes [9,11,15,42,113]; adaptive refinement strategies; and the implementation of immersed boundary techniques [15]. These developments have given rise to EFG-based adapted formulations such as the Variational Multiscale EFG (VMEFG) [109,111,119], the Generalised EFG (GEFG) [67,99], the Reduced Integration Penalty EFG (RIP-EFG) [11,15,106,107], the SUPG-stabilised EFG [9,11,15,42], and the Penalty-Based Immersed Boundary EFG (PBIB-EFG) [15].

The analogies between FEM and EFG formulations have also fostered the emergence of hybrid approaches that combine the advantages of both methods within a unified computational framework. Since most EFG formulations employ shape functions that do not satisfy the Kronecker delta property, the direct imposition of essential (Dirichlet-type) boundary conditions and the coupling with standard FEM approximations require special attention. Hybrid FEM–EFG techniques were therefore developed to exploit the flexibility and higher accuracy of EFG in regions demanding fine resolution, while retaining the computational efficiency of FEM elsewhere. Early coupling strategies, such as the non-overlapping procedure proposed by Belytschko et al. [161], used interface elements to bridge displacement fields and compensate for the lack of Kronecker delta property in MLS approximations. Although these methods did not ensure perfect continuity of shape function derivatives across the interface, the effect on overall accuracy was limited as the coupling affected only a small portion of the domain. Subsequent developments sought to remove the need for explicit interface elements, including hierarchical IEFG–FEM blending schemes [162], Lagrange multiplier-based couplings [163], and integration constraint formulations [164]. More recent advances have achieved seamless direct coupling by exploiting shape functions satisfying the Kronecker delta property, either fully or weakly. Zambrano-Carrillo et al. [19] employed Moving Kriging (MK) functions, Ullah et al. [165,166] used maximum-entropy (max-ent) functions with weak Kronecker delta behaviour at the boundaries, and Zhang et al. [135] introduced MLS approximations with linearly varying support sizes to recover the property at the interface. Beyond hybrid coupling, the lack of the Kronecker delta property in standard EFG formulations has given rise to alternative strategies conceived to allow the imposition of essential (Dirichlet-type) boundary conditions when direct prescription of nodal values is not possible. In this context, penalty-based methods [12–16], Lagrange multipliers [12,72,73,78,79], and Nitsche’s [11,55,144] method have been employed, with the latter providing a more consistent and accurate enforcement without requiring excessively large penalty parameters. Alternatively, the use of shape functions based on such as Moving Kriging (MK) [145,146,160], interpolating MLS [60,83,86], or Local Maximum Entropy (LME) [43,136,155,165,166] formulations enables the straightforward imposition of prescribed nodal values, while retaining the intrinsic advantages of EFG in terms of high-order continuity and flexible node placement.

The present review aims to provide a detailed yet conceptually accessible discussion of the main theoretical and computational aspects underpinning the development and practical implementation of EFG formulations. In particular, attention is given to the construction of shape functions—including MLS, IMLS, MK, and LME approximations—their numerical integration and computational efficiency, and the strategies to impose essential boundary conditions, including penalty, Lagrange multiplier, and Nitsche-based approaches. Moreover, the review addresses the extension of EFG methods from classical linear elasticity problems to more complex scenarios, such as non-linear heat transfer with phase change, coupled fluid–heat transport, and multiphysics applications, highlighting the tailored formulations designed for such purposes. Finally, recent Chimera-Type hybrid FEM–EFG approaches are discussed [19,42], emphasising their ability to achieve seamless coupling and accurate transfer of field variables between overlapping domains without requiring a prescribed topological relationship.

2  Shape Functions in EFG Methods

Element-Free Galerkin (EFG) formulations are developed upon the weak form of the governing equations, which necessitates the construction of shape functions for the spatial approximation of the field variables. Unlike FEM, the construction of shape functions in EFG methods is not subjected to a geometric parametrisation of the problem domain at the local or elemental level. There is no need for a prescribed connectivity, and the evaluation of physical variables is performed within a moving local support domain by continuously capturing and weighting the information of neighbouring nodes involved in the construction of the shape functions at each point of interest. This provides greater flexibility in node placement and facilitates the treatment of problems involving large deformations, evolving boundaries, or discontinuities. Nevertheless, it has also introduced several numerical challenges concerning the imposition of essential boundary conditions and the computational cost associated with neighbour searching and matrix assembly.

In this section, several advanced techniques for constructing shape functions in EFG methods are reviewed. The discussion begins with the MLS approach, which constitutes the original basis of the EFG method, and proceeds with its subsequent enhancements, including the improved MLS, interpolating MLS, and complex-variable-based MLS variants. Following these, the Reproducing Kernel Particle Method (RKPM) is examined. Although originally conceived as a standalone framework, it has been recently recognised as a procedure for shape functions construction within EFG frameworks. Further developments, such as MK and LME approximations, are also examined, with particular emphasis on their fulfilment of the Kronecker delta property, which is fully satisfied in MK formulations and only weakly enforced at the boundaries in max-ent schemes.

2.1 Moving Least Squares (MLS) Approximations

The MLS approximations constitute the foundational technique for the construction of shape functions in EFG formulations, and still remain one of the most widely employed schemes for that purpose. Its fundamental idea is to approximate each scalar field variable u(x→) at any point x→ of a given computational domain Ω as a locally weighted least-squares fit of nodal parameters within a moving local support domain.

The MLS approximation of a scalar field is expressed as [30,90]

uh(x→)=qT(x→)a(x→)∀x→∈Ω(1)

where q(x→)=[1,x,y,z,…]T is a complete basis of polynomial terms (typically up to first or second order), and a(x→) is a vector of coefficients that depend on the local coordinate x→. These coefficients are determined by minimising a weighted discrete least-squares functional [30,90]

J(a(x→))=∑I=1nW^(x→−x→(I))[qT(x→(I))a(x→)−u(I)]2,(2)

where n is the number of nodes I supporting the approximation at x→. W^(x→−x→(I)) is a non-negative weight function of the node I–evaluated at x→–, whereas uI is the corresponding nodal parameter.

Minimising J with respect to a(x→) yields the following system of equations:

A(x→)a(x→)=B(x→)u,(3)

with

A(x→)=∑I=1nW^(x→−x→(I))q(x→(I))qT(x→(I)),andB(x→)=[W^(x→−x→(1))q(x→(1)),…,W^(x→−x→(n))q(x→(n))].(4)

Assuming A(x→) is invertible, the coefficient vector is obtained as

a(x→)=A−1(x→)B(x→)u,(5)

where uT=[u(1),u(2),…,u(n)]. By substituting (5) in (1), the field variable can be expressed in terms of the nodal parameters as

uh(x→)=∑I=1nφ(I)(x→)u(I),(6)

where

φ(I)(x→)=qT(x→)A−1(x→)W^(x→−x→(I))q(x→(I))(7)

are the MLS shape functions. Although these approximations do not satisfy the Kronecker delta property, they possess high smoothness (typically C1 or higher) and exhibit excellent approximation capabilities with properly chosen weight functions W^(x→−x→(I)) and support sizes. It is worth noting that the shape and dimensions of the support domains may vary depending on the nodal distribution, and also the geometric and dimensional characteristics of the problem. Circular or spherical supports are often adopted for regular [19,41] or moderately irregular [11,42] nodal layouts, whereas rectangular [41,42,107,113], cylindrical [10,11,41,42], hexahedral [9,19,90], or even polygonal and polyhedral [156,157,167] supports are preferred when the nodal spacing exhibits marked directional variations. In such cases, the support size can be independently adjusted along each coordinate direction to capture the local anisotropy of the nodal arrangement. Accordingly, both the weight function W^(x→−x→(I)) and the corresponding normalised distance measure are defined according to the chosen support geometry. The support size is typically defined as a multiple of the local nodal spacing to ensure sufficient overlap between neighbouring supports, which is essential to guarantee numerical stability and accuracy.

It is important to distinguish between the support domain and the influence domain, two notions that are often used interchangeably in the literature but have distinct meanings depending on the adopted convention [12]. The support domainΩs(x→) is defined for each point of evaluation x→ as the region encompassing all the nodes whose weight functions are nonzero at that location, i.e.,

Ωs(x→)={x→(I)∈Rd | W^(x→−x→(I))≠0}.(8)

Conversely, the influence domain Ωi(I) is associated with each node I and represents the spatial region over which its weight function contributes to the approximation,

Ωi(I)={x→∈Rd | W^(x→−x→(I))≠0}.(9)

By definition, a node x→(I) influences a point x→ if and only if x→(I)∈Ωs(x→), which is equivalent to stating that x→∈Ωi(I). This dual relationship underpins the local nature of the Moving Least Squares (MLS) approximation, ensuring that only nearby nodes contribute to the construction of shape functions at a given point. This differentiation is not merely conceptual since it has a crucial role in defining the local nodal connectivity, and consequently in the assembly of the moment matrix involved in the MLS interpolation. The concepts of influence and support domain are depicted in Fig. 1.

Figure 1: Distinction between influence and support domains. Each node x→(1), x→(2), and x→(3) has an associated influence domain within which it can exert a contribution. Nevertheless, only nodes x→(2), and x→(3) are used in to define the support domain and to construct shape functions for approximation at the interest point x→. It is worth mentioning that x→ is not necessarily a node position, but any location within the problem domain.

The weight function W^(x→−x→(I)) is crucial in the definition of the local approximation, as it determines the extent and smoothness of the nodal contribution within its influence domain. Each node possesses an associated nodal basis function that is nonzero only inside its influence domain Ωi(I), where W^(x→−x→(I)) attains its maximum at the nodal position x→(I) and decreases monotonically and smoothly to zero at the influence domain boundary. Consequently, the weight function controls both the degree of locality and the smoothness of the resulting approximation. Conversely, the shape functionsφ(I)(x→) are constructed as combinations of these nodal basis functions through the MLS formulation, defining the actual approximation of the field variable at each evaluation point.

Commonly used weight functions include Gaussian distribution-based, quadratic, cubic, and quartic splines, and the compactly supported quartic weight function, which is very similar to the cubic splines and possesses second-order reproducing capacity [12]. The specific choice of weight function and its parameters has a strong influence on the approximation accuracy, the conditioning of the moment matrix, and the overall stability of the numerical formulation.

2.1.1 Improved MLS Approximations

The mainstream concept of IMLS approximations was first introduced by Lu et al. [31], who proposed the use of locally weighted orthogonal basis polynomials in order to simplify the evaluation of the moment matrix A(x→) in Eq. (4). The construction of an orthogonal vector of basis polynomials p(x→)=[p1(x→),p2(x→),…,pn(x→)]T within each local support domain gives rise to a diagonal moment matrix, which eliminates the need for its numerical inversion at every evaluation point. This improvement considerably enhances computational efficiency while preserving the smoothness and approximation accuracy characteristic of the standard MLS approaches. Building upon this idea, Kaljevié and Saigal [32] subsequently adopted the concept of Lu et al. and referred to it as IMLS approximations. The incorporation of this approach into the EFG framework was therefore termed the Improved EFG (IEFG) method, formalising the nomenclature for further works and demonstrating the practical benefits of this procedure in terms of computational efficiency. In this sense, further studies demonstrated that this procedure to construct shape functions is 30% faster compared to the standard MLS approximations [15,41,63].

The orthogonal polynomial basis vector is constructed from the components of a standard polynomial basis vector q(x→) as

p1=1,pi=qi−∑j=1i−1(qi,pj)(pj,pj),(10)

where the inner product between any pair of polynomial terms f(x→) and g(x→) is given by

(f,g)=∑I=1nW^(x→−x→(I))f(x→(I))g(x→(I)).(11)

The orthogonal polynomial basis vector constructed using (10) yields the diagonal moment matrix with components [11,15,168]

Aij(x→)=(pi,pj)=∑I=1nW^(x→−x→(I))pi(x→(I))pj(x→(I))={0,i≠j,(pi,pi),i=j,(i,j=1,2,3,…,m).(12)

Since the direct achievement of a diagonal moment matrix allows to circumvent the need for computing its inverse, the shape function for each node I can be merely computed as

φ(I)(x→)=W^(x→−x→(I))∑i=1mpi(x→)pj(x→(I))(pi,pi),(13)

where m is the number of polynomial terms. Please note that the definitions of n (number of nodes supporting the approximation at a position x→) and m (number of polynomial terms) will remain during the entire communication.

2.1.2 Interpolating MLS Approximations

In principle, the interpolating MLS approximations can be derived from the standard MLS by using singular weighting functions [83,86,169]. The development of interpolating MLS emerged in order to overcome the lack of Kronecker delta property in the standard MLS approximations, but the singularity of the weighting functions at the corresponding node positions of the yields to ill-conditioned moment matrices. A detailed discussion concerning this inherent instability of the interpolating MLS can be found in the work of Li and Wang [169], where a stabilised approach based on shifted and scaled polynomial basis is also proposed and analysed. Improved interpolating MLS (IIMLS) approximations, able to suppress the requirement of using singular weight functions, were proposed by Wang et al. [170], and this will be the procedure to be presented in this section.

The IIMLS starts from a standard vector of polynomial basis q(x→) with q1(x→)=1, from which a new set of polynomial terms is constructed at the support of a given point x→ as

pi(x→∗)=qi(x→∗)−∑I=1nv(x→,x→(I))qi(x→(I)),(14)

with:

v(x→,x→(I))=χ(x→,x→(I))∑J=1nχ(x→,x→(J))andχ(x→,x→(I))=∏J≠I‖x→−x→(J)‖2∏J≠I‖x→(I)−x→(J)‖2,(15)

where x→∗ is a point in the local support of x→. The same transformation given in (14) is now applied to the scalar field variable u(x→), i.e.,

u~(x→∗)=u(x→∗)−∑I=1nv(x→,x→(I))u(x→(I)).(16)

It is worth noting that v(x→(I),x→(J))=δIJ and ∑I=1nv(x→,x→(I))=1. The next step consists in performing the MLS approximation of u~(x→), which is

u~L(x→,x→∗)=∑i=2mpi(x→∗)a~i(x→)⏟pT(x→∗)a^(x→).(17)

Please note that the summation (17) start at i=2, since it can be demonstrated from (14) and the partition unity property of v(x→,x→(I)) that p1(x→∗)=0. The coefficients a~i(x→) are obtained by minimisation of the functional:

J^=∑I=1nW^(x→−x→(I))[pT(x→(I))a~(x→)−u~(x→(I))]2=∑I=1nW^(x→−x→(I))[pT(x→(I))a~(x→)+∑J=1nv(x→,x→(J))u(x→(J))−u(x→(I))]2.(18)

The minimisation of J~ with respect to a~(x→) yields

dJ^da^=0=∑I=1nW^(x→−x→(I))p(x→(I))pT(x→(I))a~(x→)+∑I=1nW^(x→−x→(I))p(x→(I))[∑J=1nv(x→,x→(J))u(x→(J))−u(x→(I))],(19)

which can be rewritten as

A(x→)a~(x→)=B(x→)[I−V(x→)]u∴a~(x→)=A−1(x→)B(x→)[I−V(x→)]u(20)

where:

A(x→)=∑I=1nW^(x→−x→(I))p(x→(I))pT(x→(I)),B(x→)=[W^(x→−x→(1))p(x→(1)),…,W^(x→−x→(n))p(x→(n))],(21)

whereas V(x→) is a n×n matrix made of n repetitions of the row

v(x→)=[v(x→,x→(1)),v(x→,x→(2)),…,v(x→,x→(n))],(22)

and I is the n×n identity matrix. Accordingly, the local approximation u~L(x→,x→∗) is

u~L(x→,x→∗)=pT(x→∗)A−1(x→)B(x→)(I−V(x→))u.(23)

Finally, the substitution of (23) in the transformation (16) yields

uh(x→)=[pT(x→)A−1(x→)B(x→)(I−V(x→))+v(x→)]⏟φ(x→)u.(24)

The achieved approximation fulfils the interpolating property and does not involve any singular weight function, so any standard weight function used in MLS can also be used for the IIMLS [149,170]. Although Wang et al. [170] claimed that the interpolating property of Eq. (24) is obvious, it is also fair to state that it is not strictly trivial. In fact, the operator (I−V(x→)) generates a vector of differences between the nodal field values and the local nodal value at x→=x→(I), i.e., (I−V(x→(I)))u=u−u(I)I. Since the MLS approximation exactly reproduces constant fields, any constant component of u (such as uII) is orthogonal to the polynomial space spanned by p(x→). Consequently, the term involving A−1(x→)B(x→)(I−V(x→)) vanishes at node x→(I), and only v(x→(I))u=u(I) remains, thus ensuring the interpolating property.

2.2 Reproducing Kernel Particle Method

Although the EFG method is traditionally built upon the MLS approximation, the Reproducing Kernel Particle Method (RKPM) has been recently recognised in the literature as a formally equivalent engine for constructing the underlying shape functions. As seen in the work of Cheng and Liew [171], RKPM is commonly treated as an independent meshfree method that utilizes a variational form to solve physical problems. On the other hand, Dehghan and Abbaszadeh [172] have formally categorised this approach as an RKPM-based EFG method. This shift acknowledges that when a global Galerkin weak form is used, the choice between MLS and RKPM becomes a matter of preference to construct the shape functions rather than a change in the fundamental numerical method.

The construction of RKPM shape functions starts from the definition of a continuous kernel approximation uh(x→) of a function u(x→) as follows [171,172]:

uh(x→)=∫ΩqT(x→−x→′)b(x→)W^(x→−x→′)u(x→′)dx→′(25)

where W^(x→−x→′) is the kernel function with a compact support, and qT(x→−x→′)b(x→) is the correction function expressed as a linear combination of the basis qT(x→−x→′). The vector of unknown coefficients b(x→) is determined by enforcing the reproducing conditions after applying a trapezoidal integration rule of the continuous (25) kernel approximation over the nodes x→(I), which yields

b(x→)=M−1(x→)H(26)

where the moment matrix M(x→) is

M(x→)=∑I=1nq(x→−x→(I))qT(x→−x→(I))W^(x→−x→(I))ΔV(I)(27)

and H=[1,0,0,…,0]T. Substituting the unknown coefficients in the trapezoidal integration of the continuous kernel approximation yields

uh(x→)=∑I=1nqT(x→−x→(I))M−1(x→)HW^(x→−x→(I))ΔV(I)⏟φ(I)(x→)u(I).(28)

In this formulation, ΔV(I) represents the volume often associated with the Voronoi cell of each node I. It is important to note that for uniform nodal distributions, ΔV(I) becomes constant and these shape functions become numerically identical to the standard MLS.

2.3 Moving Kriging Interpolation

The implementation of the MK interpolation within the EFG framework was first introduced in the pioneering work of Gu [173], allowing the construction of shape functions fulfilling the Kronecker delta property. The MK interpolation is based on a combination of a linear regression model and stochastic departures [145,160,173,174]:

uh(x→)=qT(x→)a+Z(x→),(29)

where q(x→) is the vector of polynomial basis functions, a is a vector of unknown coefficients, and Z(x→) is assumed to be a stochastic zero-mean process of variance σ2 and non-zero covariance given by

cov{Z(x→(I),x→(J))}=σ2R[R(x→(I),x→(J))],(30)

where R is the square n×n correlation matrix composed of the correlation functions R(x→(I),x→(J)) linking each pair (I,J) of the n nodes supporting the MK interpolation at x→. The most common and simplest correlation function is the Gaussian form:

R(x→(I),x→(J))=e−θ‖x→(I)−x→(J)‖2.(31)

The MK interpolation is obtained by minimising the following best linear unbiased predictor (BLUP) functional:

J[uh(x→)]=E[uh(x→)−u(x→)]2,(32)

subjected to the constraint

E[uh(x→)]=E[u(x→)],(33)

where E[(⋅)]=∫Ω(⋅)dΩ/∫ΩdΩ. The minimisation process, after substituting Eq. (29) into Eqs. (32) and (33), yields

uh(x→)=[qT(x→)A+rT(x→)B]⏟φ(x→)u,(34)

with

A=(QTR−1Q)−1QTR−1,andB=R−1(I−QA),(35)

where Q is an n×m matrix defined as

Q=[q(x→(1)),q(x→(2)),…,q(x→(n))]T,(36)

whereas the vector r(x→) has components r(I)(x→)=R(x→(I),x→), so that r(x→(J))=R. The matrix A has dimensions m×n, as can be verified by inspection of Eq. (35). The vector of shape functions φ(x→) has components

φ(I)(x→)=∑k=1mqk(x→)AkI+∑k=1nrk(x→)BkI,(37)

from which it is straightforward to prove the Kronecker delta property of φ(I)(x→), i.e., φ(I)(x→(J))=δIJ. Indeed, the value of φ(I)(x→(J)) is given by

φ(I)(x→(J))=∑k=1mqk(x→(J))AkI+∑k=1nrk(x→(J))BkI,(38)

which, by virtue of qT(x→(J))=Q and r(x→(J))=R, can be expressed in matrix form as

φ(x→(J))=QA+RB.(39)

Substituting B from Eq. (35) into Eq. (39) finally yields

φ(x→(J))=QA+RR−1(I−QA)=I,or equivalently,φ(I)(x→(J))=δIJ,(40)

thus confirming the Kronecker delta property. It should be noted that the construction of the MK shape functions is computationally more demanding than that of the standard MLS approximation. This is mainly due to the matrix inversions involved in the evaluation of both A and B in Eq. (35). In particular, the computation of A requires the inversion of the m×m matrix QTR−1Q, whereas B involves an additional n×n inversion through R−1. Consequently, the MK interpolation fulfils the Kronecker delta property, albeit at a higher computational cost. Additionally, Tu et al. [175] demonstrated that the matrices involved in the MK formulation become increasingly ill-conditioned as the nodal distribution is refined. To address this issue, they proposed a stabilised MK interpolation based on shifted and scaled polynomial basis functions:

q^T(x→)=qT(x→−x→(e)Chh),(41)

where x→(e) denotes the centroid of the n local nodes used to construct the MK interpolation at x→, i.e.,

x→(e)=1n∑I=1nx→(I).(42)

The parameters Ch and h correspond to a scaling constant and the minimum distance between any pair of the supporting nodes x→(I). The procedure to perform the MK interpolation remains the same, but with all operators defined in terms of the shifted and scaled polynomial basis. After performing the mapping q(x→)→q^(x→), Q→Q^, and R→R^, the components of Q^ and R^ are given by:

Q^=[q^(x→(1)),q^(x→(2)),…,q^(x→(n))]T,and[R^](IJ)=R^(x→(I),x→(J))=e−θ(Chh)2‖x→(I)−x→(J)‖2.(43)

According to Tu et al. [175], the optimal value of Ch for achieving the best conditioning of the matrices involved in the MK interpolation is 4.5. The correlation parameter θ has a strong influence on the accuracy of the MK approximation, whereby Zheng and Dai [174] proposed the expression θ=ω/lavg2 for this parameter. The measure lavg represents the average distance between the nodes supporting the approximation at x→, and good accuracy is commonly achieved with ω=0.1 [174–176].

2.4 Local Maximum Entropy Approximations

The construction of shape functions based on the principle of maximum entropy was first introduced by Sukumar [177] to address the underdetermined system of equations that arises in the construction of approximations over polygons with more than three sides in 2D problems, subjected to constant and linear reproducing constraints. Subsequently, Arroyo and Ortiz [178] proposed the Local Maximum Entropy (LME) formulation by introducing weight functions (or priors) into Sukumar’s original approach [177]. This modification enabled the construction of shape functions within a moving local support domain, Ωs(x→), defined by the neighbouring nodes with non-zero weight functions at the point of interest x→. Such approximations were then employed in the meshless solution of linear and non-linear elasticity problems. The weak Kronecker delta property of LME-based shape functions was also demonstrated, along with their correspondence with standard MLS approximations. Although Arroyo and Ortiz [178] used Gaussian weight functions in their pioneering work, the formulation was later generalised by Sukumar and Wright [179] to accommodate arbitrary weight functions.

The construction of LME approximations arises from the principle of determining a set of shape functions φ(I)(x→) for the nodes x→(I)∈Ωs(x→) such that

∑I=1nφ(I)(x→)=1and∑I=1nφ(I)(x→)x→(I)=x→,(44)

thereby fulfilling the constant and linear reproducing conditions [155].

For polygons with three sides (triangles) in 2D or polyhedra with four faces (tetrahedra) in 3D, Eq. (44) admits a unique solution. For cases involving a larger number of nodes defining the moving local support domain—as is typically the case in EFG methods—Eq. (44) becomes underdetermined, i.e., there are more unknowns φ(I)(x→) than available constraints. In such cases, the principle of maximum entropy is used to determine the unknown shape functions φ(I)(x→) by solving the following optimisation problem [165,180–182]:

maxφ[H(φ(x→))]=−∑I=1nφ(I)(x→)ln⁡(φ(I)(x→)W^(x→−x→(I))),(45)

subject to the constant and linear reproducing constraints in Eq. (44).

The rationale behind the LME formulation (45) is to ensure that the resulting shape functions φ(I)(x→) fulfil the reproducing conditions without introducing any additional, unwarranted assumptions about their values. In other words, the maximum-entropy solution represents the smoothest and most objective approximation consistent with the available information [178]. The solution of (44) subject to (45) yields [165,180–182]

φ(I)(x→)=W^(x→−x→(I))e−λ→⋅x→∗(I)∑J=1nW^(x→−x→(J))e−λ→⋅x→∗(J),(46)

where x→∗(I)=x→(I)−x→ are the shifted coordinates of the supporting nodes x→(I)∈Ωs(x→), and λ→ is the vector of Lagrange multipliers enforcing the linear reproducing constraints. This vector is λ→=[λx,λy] and λ→=[λx,λy,λz] in 2D and 3D problems, respectively. The vector of Lagrange multipliers is obtained from the dual of the optimisation problem (44) and (45), leading to the following non-linear system of equations [180,181]:

f(λ→)=−∑I=1nW^(x→−x→(I))e−λ→⋅x→∗(I)x→(I)∑J=1nW^(x→−x→(J))e−λ→⋅x→∗(J)=x→.(47)

This system is solved using the Newton–Raphson method, with the tangent operator (i.e., the Hessian matrix of the optimisation problem) given by [180,181]

H=∑I=1nφ(I)(x→)x→∗(I)⊗x→∗(I).(48)

A more detailed and practical explanation of the Newton–Raphson procedure for solving Eq. (47) can be found in the PhD thesis of Ullah [183].

3  Solving Numerical Difficulties in EFG Methods

Although Element-Free Galerkin (EFG) methods exhibit several noteworthy numerical advantages—such as [12,19,42] (i) the possibility of more easily achieving higher-order approximations with continuous derivatives, (ii) greater flexibility in adding or removing nodes, (iii) the absence of any post-processing requirement to obtain smooth physical fields derived from the primary variables, and (iv) improved accuracy compared with piecewise polynomial approximations used in standard FEM formulations—these methods also pose several implementation challenges.

EFG formulations can experience issues common to mesh-based techniques, such as volumetric locking in nearly incompressible media [34,106,184,185], shear locking in bending-dominated problems [31,35], and instabilities in advection-dominated transport phenomena scenarios [9,186–188]. EFG methods can also introduce additional difficulties specific to their formulation such as a higher computational cost mainly to the following aspects: the repeated inversion of local moment matrices [12,189,190], the high-order numerical integration required for non-polynomial shape functions [44–46,49], and the searching for neighbouring nodes at each integration point [39–41,43]. Other well-known issues include the need for special techniques to impose essential boundary conditions when using shape functions that do not satisfy the Kronecker delta property [11,12], the sensitivity of the results to the definition of the support domains—particularly in irregular nodal distributions [11,12,15]—and the treatment of material discontinuities [33,90,191,192]. This section discusses several strategies that have been developed to effectively address these difficulties.

3.1 Imposition of Essential Boundary Conditions and Treatment of Material Discontinuities

In EFG formulations based on shape functions that fulfil the Kronecker delta property (φ(I)(x→(J))=δIJ), the accomplishment of essential (Dirichlet) boundary conditions is straightforward, as they can be directly imposed in the form of prescribed nodal values. However, shape functions such as the MLS approximations do not possess this property. Consequently, the field variable at a node does not coincide with its nodal parameter, and essential boundary conditions cannot be imposed by merely prescribing nodal values as in mesh-based methods. To address this limitation, several specialised techniques have been developed for the imposition of Dirichlet boundary conditions. The most widely used approaches are the technique of Lagrange multiplier, the penalty method, and Nitsche’s formulation.

The method of Lagrange multipliers provides an exact enforcement of essential boundary conditions but enlarges the system of equations and may deteriorate its conditioning and sparsity. On the other hand, the penalty method preserves the system size and symmetry while maintaining a positive definite stiffness matrix. Nevertheless, its effectiveness depends on the appropriate choice of the penalty parameter. Excessively large values can lead to ill-conditioning, whereas smaller ones may result in a poor fulfilment of the essential boundary conditions. Finally, Nitsche’s method offers a more consistent and theoretically rigorous alternative by weakly enforcing the Dirichlet conditions through a symmetric modification of the weak form, providing improved accuracy and stability without introducing additional unknowns.

For the sake of generality, the discussion of these methods is presented with reference to a generic boundary value problem expressed in weak form as

a(δu,u)=ℓ(δu)∀δu∈𝒱0,(49)

where a(δu,u) denotes the bilinear form associated with the governing differential operator, and ℓ(δu) represents the corresponding linear functional.

3.1.1 Lagrange Multiplier Method

In this procedure, the weak form (49) is augmented, introducing a Lagrange multiplier field λ(x→) defined along the Dirichlet boundary Γu. The modified weak form reads [12,37,73,193]

a(δu,u)+∫Γuδλ(u−u¯)dΓ+∫ΓuλδudΓ=ℓ(δu),(50)

where λ enforces the constraint u=u¯ in a variationally consistent manner.

The Lagrange multiplier λ(x→) acts as an additional unknown field, which must be interpolated over the essential boundary. For a curvilinear coordinate s defined on Γu, the Lagrange multiplier field can be expressed as

λh(s)=∑J=1nsN(J)(s)λ(J),s∈Γu,(51)

where ns denotes the number of boundary nodes used for interpolation, and N(J)(s) are standard Lagrange interpolation functions.

Substituting the discrete forms uh of (6) and λh of (51) into the augmented weak form (50) leads to the following matrix system:

[KΛTΛ0][uλ]=[fg],(52)

where Λ represents the coupling between boundary and domain degrees of freedom.

This formulation ensures complete fulfilment of essential boundary conditions on Γu. However, the additional unknowns yield a larger system of equations. The matrix is no longer positive definite, whereby conditioning issues are introduced. These features have motivated the development of alternative, more efficient enforcement strategies, such as the penalty and Nitsche methods.

3.1.2 Penalty Method

The penalty method provides an alternative and simpler approach for the enforcement of essential boundary conditions in EFG formulations. In this technique, the constraints are imposed approximately by introducing a term that penalises any deviation of the numerical solution from the prescribed boundary values. Unlike the Lagrange multipliers technique, the penalty method does not introduce additional field variables that enlarge the algebraic system. This method instead modifies the stiffness matrix directly, maintaining its size and positive definiteness.

In the penalty approach, the weak form is augmented with a term that penalises the deviation between the approximate solution uh and the prescribed value u¯ along Γu as follows [12,13,16,106]:

a(δu,u)+∫Γuαδu(u−u¯)dΓ=ℓ(δu),(53)

where α>0 is the penalty parameter that controls the strength of the constraint enforcement. The higher the value of α, the closer u approaches u¯ on Γu. However, excessively large values may lead to ill-conditioning of the global system of equations.

The corresponding discrete system of equations can be expressed as

(K+Kα)u=f+fα,(54)

where K and f denote the global stiffness matrix and external force vector obtained from the domain integrals, whereas Kα and fα are the penalty contributions for the approximate imposition of the Dirichlet-boundary conditions.

The main advantages of the penalty method are its simplicity and computational efficiency, as the system matrix remains symmetric, positive definite, and of unchanged dimension. However, the method only enforces the essential conditions approximately. The accuracy and stability of the method is highly sensitive to an appropriate choice of the penalty parameter α since: excessively small values lead to poor constraint satisfaction, whereas very large ones can deteriorate the conditioning of the system. Nevertheless, an appropriate choice of the penalty parameter provides a suitable trade-off between numerical stability and the accuracy of boundary conditions enforcement.

3.1.3 Nitsche’s Method

Nitsche’s method can be interpreted as a consistent extension of the penalty formulation that allows the weak enforcement of Dirichlet boundary conditions, without either introducing additional unknowns or suffering from the conditioning issues associated with large penalty parameters. This approach has been successfully applied to both mesh-based and meshfree formulations due to its consistency, stability, and flexibility.

The method augments the weak form by adding terms that weakly enforce the Dirichlet condition, while preserving both the symmetry and consistency of the formulation. For a generic boundary value problem, the modified weak form reads [55,123,144,158]

a(δu,u)−∫Γuδu(κ∇u⋅n^)dΓ+∫Γu(κ∇δu⋅n^)(u−u¯)dΓ+∫Γuαδu(u−u¯)dΓ=ℓ(δu),(55)

where n^ denotes the unit outward normal vector to Γu, κ is the transport coefficient included in the diffusive term of a(δu,u), and α is a stabilisation parameter that performs a role similar to the penalty factor in Eq. (53).

The first two additional boundary terms in Eq. (55) ensure consistency of the weak formulation, whereas the third (symmetric) term provides stability. Nitsche’s formulation enforces the boundary conditions exactly in the limit of mesh refinement, without the need for excessively large stabilisation parameters.

The discrete form of Eq. (55) can be written as

(K+KN)u=f+fN,(56)

where KN and fN represent the symmetric Nitsche contributions associated with the boundary integrals. These contributions depend on both the shape functions and their normal derivatives evaluated along Γu.

The main advantages of Nitsche’s method are its consistency, improved accuracy, and the possibility of achieving optimal convergence rates without introducing additional variables. The system matrix remains both symmetric and well-conditioned, whereby this approach is often preferred over the penalty method for the weak imposition of Dirichlet boundary conditions in meshfree formulations. The stabilisation parameter α can be selected following analytical or numerical criteria to guarantee coercivity, typically as a function of the local stiffness and characteristic nodal spacing along the boundary.

In addition to the aforementioned variational methods, other strategies have been developed to address the lack of the Kronecker delta property in EFG approximations. The full transformation method [194,195] introduces a matrix that maps the meshfree degrees of freedom to the actual physical nodal values. This allows the essential boundary conditions to be imposed directly, as in conventional FEM, circumventing the need for additional interface integrals or multipliers. Alternatively, Kaljević and Saigal [32] proposed the use of singular weight functions to achieve interpolating properties with MLS approximations. Using kernels that tend to infinity at the nodal locations yields interpolating shape functions. Although this property can be enforced throughout the entire domain, its application can be specifically targeted at the boundary nodes to facilitate a straightforward enforcement of essential conditions without auxiliary variables. More recently, a variationally consistent framework has been proposed based on the reproducing kernel approach in the context of the Hellinger-Reissner variational principle [196]. This is a mixed formulation specifically developed in the framework of linear elasticity problems, where the displacement and stress fields are approximated independently to allow the essential boundary conditions to be naturally incorporated into the weak form through the stress-displacement coupling. This approach shares a similar structure with Nitsche’s method but offers the significant advantage of completely eliminating the need for the problem-dependent artificial parameters (α) required for stability in Eq. (55). Such parameter-free enforcement ensures stability and optimal convergence, representing a robust alternative for nearly incompressible elastic media and other elliptic systems where traditional stabilisation approaches may be challenging.

3.1.4 Application to Material Discontinuities

The aforementioned weak enforcement strategies can also be extended to impose interface conditions across material discontinuities. In such cases, the physical domain Ω is divided into two subdomains Ω+ and Ω−, separated by an internal interface Γi over which continuity of the field variable u and the corresponding diffusive flux −κ∇u⋅n^ must be guaranteed. The continuity conditions across Γi can be written as

⟦u⟧=u+−u−=0,⟦κ∇u⋅n^⟧=0,(57)

where ⟦⋅⟧ denotes the jump operator and n^ is the unit normal to Γi oriented from Ω− to Ω+.

The enforcement of the interface conditions in Eq. (57) follows the same principles used for Dirichlet boundaries:

•   Lagrange multipliers: The weak form is augmented with interface terms involving a multiplier field λ defined on Γi, enforcing ⟦u⟧=0 in a variationally consistent way [33,191]:

a(δu,u)+∫Γiδλ⟦u⟧dΓ+∫Γiλ⟦δu⟧dΓ=ℓ(δu).(58)

This approach ensures exact satisfaction of continuity but introduces additional unknowns associated with λ.

•   Penalty method: The interface continuity is enforced approximately by penalising the jump in u across Γi [12]:

a(δu,u)+∫Γiα⟦δu⟧⟦u⟧dΓ=ℓ(δu),(59)

where α is a penalty parameter controlling the strength of the constraint.

•   Nitsche’s method: A consistent and stable alternative is obtained by symmetrically adding interface terms analogous to those introduced for Dirichlet boundaries [197]:

a(δu,u)−∫Γi⟦δu⟧{{κ∇u}}⋅n^dΓ+∫Γi{{κ∇δu}}⋅n^⟦u⟧dΓ+∫Γiα⟦δu⟧⟦u⟧dΓ=ℓ(δu),(60)

where {{⋅}} denotes the average operator across Γi. The first two interface terms ensure consistency, whereas the last term guarantees stability.

These formulations enable a unified and flexible treatment of discontinuous material domains, allowing jumps in material parameters to be handled naturally within the weak formulation. The choice among Lagrange multipliers, penalty, or Nitsche’s approach depends on the desired trade-off between accuracy, computational cost, and conditioning of the resulting algebraic system.

It is important to remark that the bilinear and linear operators a(δu,u) and ℓ(δu) must be integrated considering the appropriate nodal support for each subdomain, as depicted in Fig. 2. The integrations over Ω− are performed using shape functions constructed only with the nodes belonging to Ω− and the interface nodes located on Γi, whereas those over Ω+ are evaluated using the nodes of Ω+ and the same shared interface nodes. This ensures that the support of each trial and test function remains confined to its corresponding subdomain. Integrals defined along the interface Γi involve quantities computed from both sides: the jump and average operators are evaluated using the nodal fields associated with Ω− (blue nodes) and Ω+ (red nodes), respectively. This distinction guarantees a consistent assembly of the weak form, preserving the correct coupling and continuity between the two material regions.

Figure 2: Representation of the computational domain divided into two subdomains, Ω− (blue) and Ω+ (red), separated by the internal interface Γi (black). The blue and red circles denote the nodes belonging exclusively to Ω− and Ω+, respectively, whereas black circles correspond to the nodes shared by both subdomains and located on Γi. Integrals over Ω− and Ω+ are evaluated using shape functions constructed from their respective nodal sets, including the interface nodes. Quantities involved in the interface integrals over Γi, such as jumps and averages, are computed using the blue-side values for the minus domain and the red-side values for the plus domain.

It is worth mentioning that the interface integrals involved in the methods of Lagrange Multipliers, Penalty, and Nitsche are primarily required for EFG formulations based on shape functions that do not fulfil the Kronecker delta property. Alternatively, material interfaces can be accurately modelled by using shape functions that possess this property at the interface Γi. While standard MLS approximations lack this feature, alternative schemes such as Moving Kriging [173], LME [178], or Interpolating MLS [170] allow for the direct nodal enforcement of continuity. The requirement ⟦u⟧=0 is naturally fulfilled through the shared nodal degrees of freedom at the interface Γi, effectively connecting the subdomains Ω+ and Ω− without the need for additional interface integrals or stabilisation parameters. In this framework, the gradient or flux discontinuity across the interface—arising from different material properties—is captured naturally by the assembly of independent shape function supports on either side of Γi. This approach significantly simplifies the implementation for multi-material problems, as the interface conditions are implicitly fulfilled by the interpolating nature of the approximation space.

3.2 Efficient Numerical Integration Schemes

In EFG methods, the accurate evaluation of domain and boundary integrals remains one of the main numerical challenges. The non-polynomial character of the shape functions and the overlap of support domains make the exact integration of the weak form virtually impossible, whereby the choice of a suitable quadrature scheme is crucial for ensuring accuracy, convergence, and computational efficiency. EFG methods often demand high-order quadrature schemes, since low-order integration rules can lead to non-converging and unstable solutions. Recent theoretical research by Wu and Wang [198] has established that this phenomenon arises from losing the Galerkin orthogonality condition due to integration errors. In such a fundamental work, it has formally been demonstrated that the overall convergence of EFG methods is governed by a synergy between the interpolation error (associated with the basis functions) and the integration error (associated with the quadrature). Consequently, the development of efficient domain integration algorithms specifically tailored for EFG methods is not only motivated by computational cost but by the mathematical requirement of matching the integration consistency with the approximation order to restore optimal convergence. The following section summarises the main approaches proposed in this context, highlighting their underlying principles and comparative performance.

The first attempt to improve the computational efficiency of domain integration in EFG methods was proposed by Beissel and Belytschko [199] in the context of elliptic problems in linear elasticity, and is known as the direct nodal integration (DNI) method. In this procedure, the integration domain associated with each node is defined through a Voronoi tessellation of the problem domain. Each Voronoi cell represents the region of integration for the corresponding node, and the integrals are approximated by single-point evaluations at the nodal positions. The integration weights are determined by the area (or volume) of the Voronoi cell for interior nodes and by the length (or surface area) of the boundary portions for nodes located on the traction boundaries. This strategy considerably simplifies the numerical integration process and removes the dependence on any auxiliary background mesh, thereby significantly reducing the computational cost. However, although the method proved highly efficient, it suffered from stability issues and exhibited sub-optimal convergence behaviour. This has motivated researchers to propose methods devoted to eliminate such instabilities, which include: adding squared residuals-based stabilisation terms to yield a residual-based DNI (RDNI) scheme [199], introducing the strain field from a first-order Taylor of the displacement field approximation to obtain a naturally stabilised DNI (NDNI) approach [200], or the least squares stabilisation term used in the modified DNI (MDNI) [201]. Despite the aforementioned improvements, these methods still inherit some fundamental limitations from the original DNI formulation. In particular, none of them is able to fully satisfy the linear patch test. This indicates that the consistency of the approximated strain field is not entirely recovered, although accurate displacement solutions can be obtained. The stress predictions remain sub-optimal, especially near the boundaries. Moreover, the computation of higher-order derivatives required in these formulations is often time-consuming and susceptible to accumulated numerical errors. A comprehensive discussion in this regard can be found in the work of Luo et al. [47]. These drawbacks have motivated the development of more robust integration strategies capable of restoring both stability and consistency, such as the so-called stabilised conforming nodal integration (SCNI) methods.

The SCNI approach introduced in the seminal work of Chen et al. [202] constitutes the first variationally consistent integration framework for EFG methods. In this approach, the shape functions gradient ∇φ(I) is replaced by a smoothed counterpart ∇~φ(I) defined within integration subdomain Ωg(I) associated with each node. It is worth remarking that the integration domain Ωg(I) of a node must not be confused with its influence domain Ωi(I), i.e., that region where φ(I)(x→)≠0. The smoothed operator is defined as an averaged (or “smoothed”) gradient obtained as

∫Ωg(I)∇~φ(I)(x→)dΩ=∫Γg(I)φ(I)(x→)⋅n^dΓ,(61)

where φ(I) denotes the shape function associated with node I, and n^ is the outward unit normal to the boundary Γg(I) of the nodal integration subdomain. The equality of Eq. (61) is known as the integration constraint, and ensures that the numerical solution of the weak form is consistent with the Galerkin formulation. This means that linear fields are exactly reproduced within each nodal subdomain, which guarantees first-order consistency of the discretisation and thus fulfilment of the linear patch tests. By satisfying this constraint, the SCNI method eliminates the instability and suboptimal convergence behaviour of direct nodal integration. However, the smoothed strain evaluated at the centre of each nodal integration cell can only reproduce a constant strain field within. The subdomain and corresponding boundaries of each node x(I) are commonly defined according to the conforming Voronoi diagram of the nodes distribution, as depicted in Fig. 3. Consequently, SCNI can only provide accuracy and convergence comparable to linear finite elements even if a quadratic approximation is used. Duan et al. [203] introduced the Quadratically Consistent Integration (QCI) framework to overcome such a limitation, which has been achieved by extending the variationally consistent concept of SCNI via the enforcement of both the differentiation of approximation consistency (DAC) and the discrete divergence consistency (DDC). The DAC condition requires that the derivatives of the shape functions φ(I) preserve the consistency of the polynomial basis p(x→), expressed as

∇p(x→)=∑I=1np(x→(I))∇φ(I)(x→),(62)

which ensures the exact reproduction of the chosen basis under differentiation. In turn, the DDC imposes that the numerical derivatives ∇φ(I) satisfy the divergence theorem in a discrete sense, providing compatibility between φ(I) and their gradients:

∫Ωg(I)∇φ(I)(x→)q(x→)dΩ=∫Γg(I)φ(I)(x→)q(x→)n^dΓ−∫Ωg(I)φ(I)(x→)∇q(x→)dΩ,(63)

where q(x→) is a smooth test function constructed from the derivatives of the polynomial basis p(x→).

Figure 3: Schematic representation of nodal integration subdomains defined through a conforming Voronoi tessellation of the node distribution. In the SCNI framework, each Voronoi cell serves as a smoothing domain Ωg(I) for the evaluation of the smoothed gradient ∇~φ(I). Similar Voronoi-based partitions have been subsequently adopted in other integration schemes such as QCI, RKGSI, and CEFG, owing to their convenient and mesh-independent definition of nodal subdomains.

In its original form, the method employed three quadrature points per integration cell (QC3) to compute smoothed nodal derivatives that fulfil both the quadratic DAC and DDC. Subsequently, Duan et al. [204,205] generalised this concept to a one-point quadrature scheme (QC1) by preserving quadratic consistency via the incorporation of higher-order derivatives evaluated at the cell centre. This formulation maintains second-order consistency while significantly reducing computational cost, and can exactly pass both the linear and quadratic patch tests. Based on these developments, a more general integration framework was later proposed by Chen et al. [51] under the name Variationally Consistent Integration (VCI). This approach extends the concept of variational consistency beyond the quadratic level achieved by QCI, providing a unified variational basis for constructing integration schemes of arbitrary order. Instead of relying on specific smoothing or correction procedures, VCI enforces that the discrete weak form exactly reproduces the continuous variational equations for polynomial fields up to a prescribed degree. Accordingly, it guarantees the exact satisfaction of higher-order patch tests while maintaining computational efficiency and stability. The SCNI and QCI schemes can be regarded as particular cases of the general VCI framework for first- and second-order consistency, respectively. Although the VCI framework provides a robust, theoretically general formulation, its practical implementation is considerably more involved than in SCNI or QCI. The method requires solving local variational problems within each integration subdomain to determine the consistent projection operators, which increases the formulation complexity and computational cost. A more general variationally consistent approach is provided by the Consistent Element-Free Galerkin (CEFG) method [52,53], which extends the ideas of SCNI and QCI to the framework of linear, quadratic, and cubic approximations based on the Hu–Washizu three-field variational principle. In this formulation, three independent fields are considered: the displacement u, an interpolated strain ε~, and an assumed Cauchy stress σ^. The main purpose of the method is enforcing the orthogonality condition as a discrete level

∫Ωg(I)σ^T:(ε−ε~)dΩ=0,(64)

which ensures that the Hu-Washizu weak form reduces to the standard Galerkin formulation. The CEFG method constructs the interpolated strain ε~ using corrected nodal derivatives, which are computed such that the orthogonality condition holds for each integration subdomain. In practice, this involves solving local projection problems over the integration cells. This is similar to SCNI and QCI, but generalized to higher-order approximations. The CEFG scheme naturally reduces to SCNI and the QCI formulations when using linear and quadratic polynomial bases, respectively. The underlying rationale of the CEFG is that enforcing the Hu–Washizu orthogonality condition allows the strain interpolation to be variationally consistent, independent of the order of approximation.

The integration schemes discussed thus far have been primarily developed and validated for second-order elliptic partial differential equations, where the weak form involves first-order derivatives. However, the requirement for variational consistency extends to higher-order elliptic systems [206,207]. The most notable example is the fourth-order elliptic equations governing thin-plate bending, where the Sub-domain Stabilized Conforming Integration (SSCI) proposed by Wang and Chen [208] represented a significant evolution of the SCNI concept. While standard SCNI ensures consistency for first-order derivatives, SSCI partitions the smoothing cells into triangular sub-domains to satisfy the integration constraints required for bending exactness [208,209]. This allows the EFG framework to maintain spatial stability and pass the pure bending test, demonstrating that stabilized integration principles can be successfully tailored for complex mechanical systems involving second-order derivatives in their weak formulation.

It should be noted that the aforementioned integration schemes have been primarily developed and validated in the framework of elliptic problems, such as linear elasticity and Poisson-type equations. The variational principles and consistency conditions underpinning these methods rely on the symmetric and positive-definite nature of the associated differential operators, which facilitates both the satisfaction of patch tests and the convergence of the numerical solution.

For problems beyond this class, including non-symmetric, advection-dominated, or multiphysics systems, the direct application of SCNI, QCI, or CEFG is generally limited. To overcome these restrictions, alternative projection-based integration schemes have been proposed. Among these, the Reproducing Kernel-based Gradient Smoothing Integration (RKSGI) [54] and the Projection Integration (PI) [57] methods stand out as flexible approaches capable of handling a broader spectrum of governing equations. In the RKGSI framework, a smoothed gradient of a nodal shape function φ(I)(x→) is constructed over non-overlapping integration cells as Ω≊∑l=1ncΩl, with the nodes as cell vertices. This smoothed gradient is obtained through a reproducing kernel approximation defined as

∇~φ(I)(x→)=∫Ωl(y→)qT(y→)a~(x→)∇φ(I)(y→)dΩ=a~(x→)∫Ωl(y→)q(y→)∇φ(I)(y→)dΩ,(65)

which allows the derivative of the field to be consistently projected over the nodal integration domain, preserving the compatibility with the discrete divergence. The coefficients vector a~(x→) is obtained from the reproducing conditions as a~(x→)=qT(x→)G−1, where

G=∫Ωl(y→)qT(y→)q(y→)dΩ.(66)

According to this, the smoothed gradient is given by

∇~φ(I)(x→)=qT(x→)G−1∫Ωl(y→)q(y→)∇φ(I)(y→)dΩ,(67)

which after integration by parts yields

∇~φ(I)(x→)=qT(x→)G−1[∫Γl(y→)q(y→)φ(I)(y→)n^dΓ−∫Ωl(y→)∇q(y→)φ(I)(y→)dΩ,].(68)

It is worth noting that the smoothed gradient constructed in Eq. (68) is specifically designed to satisfy the compatibility condition of Eq. (63) over each integration cell Ωl, rather than over the smoothing domains associated with the nodes as in conventional nodal integration schemes. The integration by parts transformation eliminates the need for explicitly computing the standard MLS derivatives, since the smoothed gradient is directly obtained from boundary and volumetric integrals involving only the shape function itself. This procedure yields a compatible gradient field that preserves the discrete divergence constraint while enabling the use of low-order quadrature rules, in a manner analogous to FEM. From a computational standpoint, the construction of the smoothed gradient proceeds locally within each integration cell Ωl, and is performed in a cell-wise manner as follows:

1.   Identify the vertices x→(I) associated with the integration cell Ωl.

2.   Compute the moment matrix G over Ωl, and evaluate the volumetric and boundary integrals in Eq. (68) for each node contributing to the cell.