An inherently consistent reproducing kernel gradient smoothing framework toward efficient Galerkin meshfree formulation with explicit quadrature

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 349 (2019) 628–672 www.elsevier.com/locate/cma

An inherently consistent reproducing kernel gradient smoothing framework toward efficient Galerkin meshfree formulation with explicit quadrature Dongdong Wang ∗, Junchao Wu Department of Civil Engineering, Xiamen University, Xiamen, Fujian 361005, China Received 21 November 2018; received in revised form 20 February 2019; accepted 20 February 2019 Available online 2 March 2019

Highlights • • • • •

A reproducing kernel gradient smoothing framework is proposed for meshfree methods. Integration consistency is an inherent property of the proposed methodology regardless of integration schemes. Explicit quadrature rules referring to 2D triangular and 3D tetrahedral integration cells are developed. Total number of sample points for the present quadrature rules is minimized from a global point of view. The proposed framework enables efficient Galerkin meshfree analysis with optimal convergence, particularly for higher order basis functions.

Abstract A reproducing kernel gradient smoothing framework with explicit quadrature rules is proposed for efficient Galerkin meshfree formulation. In this framework, the meshfree smoothed gradients are formulated via a reproducing kernel representation of the standard gradients of field variables. It is interesting to find that this reproducing kernel construction of smoothed gradients of meshfree shape functions enables an identical satisfaction of the integration constraint derived from Galerkin meshfree formulation. In other words, the integration consistency is an inherent property of the proposed reproducing kernel gradient smoothing methodology regardless of integration schemes. Consequently, without violation of the integration constraint, the conventional normal low order Gauss quadrature rules in finite element analysis now can be used to properly integrate the meshfree stiffness matrix simply through replacing the standard meshfree gradients with the reproducing kernel smoothed gradients at Gauss points. In order to efficiently compute the reproducing kernel smoothed gradients, a set of explicit quadrature rules referring to 2D triangular and 3D tetrahedral integration cells are systematically developed. The total number of sample points associated with these quadrature rules is minimized from a global point of view through introducing an equivalent number of sample points, which implies as many as sample points are shared by neighboring integration cells. The proposed methodology recovers the stabilized conforming nodal integration when linear basis function is employed, while the focus of the present work is multidimensional higher order basis functions such as quadratic and cubic ones. Superior convergence, accuracy as well as efficiency performances of the proposed reproducing kernel gradient smoothing framework are thoroughly demonstrated by a series of numerical examples. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Meshfree method; Numerical integration; Integration consistency; Reproducing kernel gradient smoothing; Higher order basis function; Efficiency

∗

Corresponding author. E-mail address: [email protected] (D. Wang).

https://doi.org/10.1016/j.cma.2019.02.029 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

629

1. Introduction Galerkin meshfree methods have attracted significant research attention and been widely used due to their excellent accuracy and sound stability characteristics [1,2]. An earlier Galerkin meshfree formulation is the diffuse element method developed by Nayroles et al. [3], which employs the moving least squares (MLS) [4] approximation and diffuse derivatives of shape functions. Later on, Belytschko et al. [5,6] proposed the elementfree Galerkin method via introducing full derivatives of MLS shape functions, background cell-based higher order Gauss integration and Lagrange multiplier enforcement of essential boundary conditions. On the other hand, through incorporating the correction function and reproducing conditions into the kernel estimation, Liu et al. [7] presented a reproducing kernel (RK) meshfree approximation and the related reproducing kernel particle method, whose discrete consistency was illustrated by Chen et al. [8]. It turns out that MLS and RK approximants are equivalent when monomial basis functions are used. Since then, Galerkin meshfree formulations have enjoyed great popularity and many versatile methods have been developed over the years, for instance, h-p cloud method [9], partition of unity method [10], natural element method [11], meshless local Petrov–Galerkin method [12], finite sphere method [13], radial point interpolation method [14], generalized finite-element method (GFEM) [15], reproducing kernel interpolation method [16], moving kriging interpolation element-free Galerkin method [17], reproducing kernel element method [18], max-entropy meshfree method [19,20], smoothed finite element method [21], optimal transportation meshfree method [22], semi-Lagrangian reproducing kernel particle method [23], generalized meshfree approximation method [24], reproducing kernel meshfree peridynamics [25], quasi-convex reproducing kernel meshfree method [26], direct displacement smoothing method [27], hierarchical partition of unity element-free Galerkin method [28], quasi-linear reproducing kernel particle method [29] and smoothed particle Galerkin method [30], among others. More detailed summaries on the advances of Galerkin meshfree methods can be found in [1,2,31–36]. For the weak form-based Galerkin meshfree methods, domain integration is inevitable. However, unlike the conventional polynomial type of finite element shape functions [37,38], the meshfree shape functions such as MLS/RK meshfree shape functions have a rational nature and their influence or support domains are usually overlapped each other [39]. Consequently, the normal low order Gauss quadrature rules designed for the polynomial shape functions are not well suitable for Galerkin meshfree methods and time-consuming higher order Gauss integration rules are often required to maintain the solution accuracy and convergence behavior [39–42]. Thus the development of efficient domain integration algorithms has been a very active topic for Galerkin meshfree methods [2,40–44]. Bessiel and Belytschko [45] presented a nodal integration method to improve the computational efficiency, where the nodally integrated weak form is stabilized via adding the equilibrium residual term multiplied by an artificial parameter. Employing additional stress points together with meshfree nodes as the integration locations offers another choice to retain the stability of Galerkin meshfree methods [46–49]. Moreover, the support domain integration method was also investigated in several works [12,13,50–53]. This method conveniently takes the support domains of meshfree shape functions as the integration cells, while in each cell higher order Gauss quadrature rules are still necessary to properly integrate the weak form. Along a different path, Chen et al. [40,41] presented the so-called integration constraint corresponding to the linear patch test for Galerkin meshfree methods, where it is shown that besides the linear reproducing or consistency condition of meshfree shape functions, the integration constraint has also to be satisfied for the domain integration of Galerkin meshfree methods in order to exactly pass the linear patch test. As a result, they proposed the stabilized conforming nodal integration (SCNI) by means of strain smoothing [54] for Galerkin meshfree methods. Both efficiency and stability are achieved in SCNI. Subsequently, SCNI has been further developed and generalized to many other problems, for examples, the locking-free plate and shell problems [55–59] large deformation failure simulations [41,60–62] and convection dominated problems [63], etc. A sub-domain stabilized conforming integration in conjunction with Hermite reproducing kernel meshfree approximation was introduced by Wang et al. [64–68] to deal with fourth order thin plate and shell problems. To accommodate the excessively large deformation analysis, Chen et al. [23,69] proposed a stabilized non-conforming nodal integration equipped with semi-Lagrangian meshfree formulation, where a regular domain is taken as the nodal representative cell for strain smoothing. Based upon SCNI, Duan et al. [70,71] investigated the integration constraints for higher order meshfree approximations and presented a quadratically consistent integration (QCI) method that is second order accurate for Galerkin meshfree formulations. In QCI, the smoothed derivatives of meshfree shape functions are not explicitly

630

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

constructed like those in SCNI, and a system of equations derived from the quadratic integration constraint need to be solved to attain the smoothed derivatives with respect to a specific coordinate component. For example, in 3D case, the meshfree smoothed gradients at a sampling point require solving three sets of 4 by 4 system of equations, which will introduce undesirable additional computational cost [71]. Meanwhile, an arbitrary order variationally consistent integration (VCI) has been proposed by Chen et al. [72,73] for Galerkin meshfree methods. The method of VCI enables the construction of various integration schemes up to the order of completeness in the meshfree approximation and it can also be used to correct other integration algorithms such as Gauss integration and stabilized non-conforming nodal integration [74]. Nonetheless, VCI falls into the category of Petrov–Galerkin formulation that involves non-symmetric stiffness matrices. Wang and Wu [75] proposed a quadratically exact nesting sub-domain gradient smoothing integration (NSGSI) algorithm for the stiffness matrix evaluation, which is established through optimally combining the contributions from the two-level triangular nesting sub-domains. It turns out that this approach is quite efficient, but its generalization to three dimensional formulation and higher order shape functions is not straightforward. Moreover, gradient stabilized meshfree nodal integration methods have also been developed, for example, the naturally stabilized nodal integration by Hillman and Chen [76] and the displacement smoothing stabilized nodal integration by Wu et al. [27], etc. In this study, we aim to develop a general reproducing kernel gradient smoothing framework that can provide explicit quadrature rules for efficient Galerkin meshfree formulation. In this framework, the meshfree smoothed gradients are represented by the standard gradients of field variables with a reproducing kernel formulation. This is different to the gradient reproducing formulation in [77–80], where the meshfree smoothed gradients are expressed by the field variables themselves. By construction, the smoothed gradients automatically satisfy the integration constraint resulting from Galerkin formalism. Thus an inherent integration consistency is embedded into the proposed meshfree smoothed gradients irrespective of integration methods. Consequently, similar to the finite element stiffness matrix integration, the conventional normal low order Gauss quadrature rules designed for polynomial shape functions [37,38] now can be employed to evaluate the proposed meshfree stiffness matrix without violation of the integration constraint. The only difference is to replace the standard gradients at the Gauss points with the proposed reproducing kernel smoothed gradients of meshfree shape functions. Subsequently, a set of explicit quadrature rules referring to 2D triangular and 3D tetrahedral integration cells are devised to compute the reproducing kernel smoothed gradients. In these quadrature rules, the number of integration sample points is minimized by introducing an equivalent number of sample points, as significantly improves the computational efficiency since as many as sample points are shared by neighboring integration cells. The proposed integration strategy is termed as reproducing kernel gradient smoothing integration (RKGSI). The formulation of RKGSI is applicable to arbitrary order multidimensional Galerkin meshfree methods. RKGSI will reduce to the method of SCNI when a linear basis function is used, and in this work particular emphasis is the Galerkin meshfree formulations with higher order basis functions such as quadratic and cubic basis functions. The efficacy of RKGSI is validated by numerical results. The remainder of this paper is organized as follows. Section 2 elaborates the Galerkin meshfree formulation with RK approximation and the integration constraint for numerical integration. Subsequently, a reproducing kernel gradient smoothing framework is presented in Section 3, where the inherent integration consistency is discussed. In Section 4, a set of explicit quadrature rules with particular reference to 2D triangular and 3D tetrahedral integration cells are developed to efficiently compute the reproducing kernel smoothed meshfree gradients, in which the number of sampling points is minimized from a global point of view. A series of numerical examples ranging from 1D to 3D are given in Section 5 to demonstrate the proposed methodology. Concluding remarks are finally drawn in Section 6. 2. Galerkin meshfree formulation and integration constraint 2.1. RK meshfree approximation For convenience of the subsequent development, the continuous form of reproducing kernel (RK) approximation defined in a problem domain Ω is first considered. According to the RK theory [7,8], the continuous RK approximant of a field variable u(x), denoted by u r (x), can be expressed as: ∫ u r (x) = c[ p] (x, s)φ(x, s)u(s)dΩ (1) Ω(s)

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

631

in which φ(x, s) is the non-negative kernel function located at s, φ(x, s) has a nodal influence domain Ω (s) and it will vanish when x ∈ / Ω (s). In this study, the classical cubic B-spline function is selected for the kernel function, in 1D case, the kernel function φ(x, s) takes the following form [1]: ⎧ 2 1 ⎪ ⎪ ⎪ − 4ℓ2 + 4ℓ3 ℓ≤ ⎪ ⎪ 2 ⎨3 |s − x| (2) φ(x, s) = 4 − 4ℓ + 4ℓ2 − 4 ℓ3 1 < ℓ ≤ 1 , ℓ = ⎪ r ⎪ 3 3 2 ⎪ ⎪ ⎪ ⎩ ℓ>1 0 where r denotes the support size measuring the influence domain of kernel function. A multidimensional kernel function φ(x, s) can be conveniently constructed via the tensor product of 1D kernel functions in different directions. c[ p] (x, s) is namely the correction function, which consists of a pth order monomial basis vector p[ p] and an unknown vector a: c[ p] (x, s) = p[ p]T (s)a(x)

(3)

where: { }T p[ p] (x) = 1, x, y, z, . . . , x 2 , y 2 , z 2 , . . . , x i y j z k , . . . , z p , 0 ≤ i + j + k ≤ p

(4)

The position dependent unknown vector a(x) in Eq. (3) can be solved by enforcing the pth order reproducing or consistency condition [7,8]: ∫ c[ p] (x, s)φ(x, s) p[ p] (s)dΩ = p[ p] (x) (5) Ω(s)

Substituting Eq. (3) into Eq. (5) leads to: A(x)a(x) = p[ p] (x)

(6)

where A(x) is the moment matrix defined as: ∫ p[ p] (s) p[ p]T (s)φ(x, s)dΩ A(x) =

(7)

Ω(s)

Thus from Eq. (6), we obtain the relationship of a(x) = A−1 (x) p[ p] (x) and accordingly the continuous RK approximant in Eq. (1) then becomes: ∫ −1 r [ p]T p[ p] (s)φ(x, s)u(s)dΩ (8) u (x) = p (x) A (x) Ω(s)

In discrete RK meshfree formulation, the problem domain Ω and the associated boundary Γ are discretized by a set of nodes or particles {x I } NP I =1 and each node x I is attached with a kernel function φ(x, x I ) and a meshfree shape function Ψ I (x). A discrete counterpart of the RK approximation in Eq. (8) can be established as [7,8]: u h (x) =

NP ∑

Ψ I (x)d I

(9)

I =1

with Ψ I (x) = p[ p]T (x) A−1 (x) p[ p] (x I )φ(x, x I ) A(x) =

NP ∑

p[ p] (x I ) p[ p]T (x I )φ(x, x I )

(10) (11)

I =1

It is straightforward to verify that the RK meshfree shape function Ψ I (x) of Eq. (10) meet the following pth order discrete reproducing or consistency condition [7,8]: NP ∑ I =1

Ψ I (x) p[ p] (x I ) = p[ p] (x)

(12)

632

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Meanwhile, Eq. (12) also implies the following reproducing condition for the meshfree shape function gradient Ψ I,i : NP ∑

[ p]

Ψ I,i (x) p[ p] (x I ) = p,i (x)

(13)

I =1

where the subscript comma means a partial differentiation with respect to given spatial dimensions. 2.2. Galerkin meshfree formulation Without loss of generality, we consider the following scalar potential problem as a model problem: ⎧ ⎨ki j u ,i j + b = 0 in Ω t = ki j n i u , j = t on Γ t ⎩ u=u on Γ g

(14)

where ki j is the material constant and b denotes the source term. t and u are the prescribed values on the natural boundary Γ t and the essential boundary Γ g , respectively. n i represents the component of the outward normal vector n of Γ t . The weak form of Eq. (14) is: ∫ ∫ ∫ δu ,i ki j u , j dΩ = δutdΓ + δubdΩ (15) Γt

Ω

Ω

Introducing the meshfree approximation of Eq. (9) into the weak form of Eq. (15) yields the following discrete Galerkin meshfree equation: Kd = f

(16)

where K is the stiffness matrix, d and f are the nodal coefficient vector and force vector, which can be assembled from their corresponding nodal contributions: NP

NP

NP

K = A [K I J ]; f = A [ f I ]; d = A [d I ] I,J =1

I =1

(17)

I =1

with A denoting the standard local–global assembly operator [37]. K I J and f I are given by: ∫ KI J = B TI D B J dΩ , B I = {Ψ I,i }, D = [ki j ] Ω ∫ ∫ Ψ I (x)b(x)dΩ Ψ I (x)t(x)dΩ + fI = Γt

(18) (19)

Ω

2.3. Integration constraint The integration constraint is the condition required for the Galerkin discrete formulation of Eq. (16) to exactly reproduce the solutions spanned by the basis vector p[ p] . Since p[ p] contains pth order monomial basis functions herein, any solution u(x) and its gradient u ,i (x) represented by p[ p] can be expressed as follows: u(x) = p[ p]T (x)α [ p]T

u ,i (x) = p,i

(20)

(x)α = p[ p−1]T (x)β i

(21)

where p is ( p − 1)th order monomial vector, α and β i are constant vectors. According to Eqs. (14), (20) and (21), the corresponding meshfree nodal coefficient d I , boundary flux t(x) and source term b(x) are: [ p−1]

d I = u(x I ) = p[ p]T (x I )α

(22) [ p]T

t(x) = ki j n i u , j (x) = ki j n i p, j (x)α = ki j n i p[ p−1]T (x)β j [ p]T

[ p−1]T

b(x) = −ki j u ,i j (x) = −ki j p,i j (x)α = −ki j p,i

(x)β j

(23) (24)

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

633

Now we would like that Eq. (16) is exactly satisfied with the quantities defined by Eqs. (22)–(24). Thus plugging Eqs. (22)–(24) into the left and right hand sides of Eq. (16) yields: ∫ NP ∑ ∑ NP K d = k Ψ (x) Ψ J, j (x) p[ p]T (x J )αdΩ I J J i j I,i J =1 Ω J =1 ∫ [ p]T (25) = ki j Ψ I,i (x) p, j (x)dΩ α ∫Ω = ki j Ψ I,i (x) p[ p−1]T (x)dΩ β j ∫ ∫ Ω [ p−1]T (x)dΩ ]β j (26) f I = ki j [ Ψ I (x)n i (x) p[ p−1]T (x)dΓ − Ψ I (x) p,i Γ

Ω

Consequently, an enforcement of the equality between Eqs. (25) and (26), namely, satisfying the discrete Galerkin meshfree equation of Eq. (16), leads to: ∫ Ψ I,i (x) p[ p−1] (x)dΩ = g i I (27) Ω

with ∫ Ψ I (x)n i p

gi I =

[ p−1]

∫

[ p−1]

(x)dΓ − Ω

Γ

Ψ I (x) p,i

(x)dΩ

(28)

Eq. (27) is the desired integration constraint which needs to be fulfilled for the numerical integration of Galerkin weak form. 3. Reproducing kernel gradient smoothing method 3.1. RKGS formulation In this section, a reproducing kernel gradient smoothing (RKGS) framework that identically meets the integration constraint of Eq. (27) is developed. In other words, the integration consistency is an inherent property of the present method. The proposed methodology begins with a reproducing kernel gradient smoothing representation of the field variable gradient u ,i : ∫ ˜(x, s)u ,i (s)dΩ ˜ u r,i (x) = c[ p−1] (x, s)φ (29) Ω(s)

where c[ p−1] (x, s) is the correction function for the reproducing kernel gradient approximation, ˜ u r,i is the reproducing ˜ kernel smoothed gradient corresponding to the standard gradient u ,i . φ (x, s) is the kernel function for the gradient reproducing operation. It is noted that the reproducing kernel gradients in [77–80] are expressed by the field variable u itself on the right hand side of Eq. (29), while the present reproducing kernel gradient smoothing is built upon the gradient of the field variable. Subsequently, in accordance with the reproducing kernel formulation [7,8], the correction function c[ p−1] (x, s) takes the following form: c[ p−1] (x, s) = p[ p−1]T (s)˜ a(x)

(30)

where the unknown vector ˜ a can be determined through imposing the ( p − 1)th order consistency condition: ∫ ˜(x, s) p[ p−1] (s)dΩ = p[ p−1] (x) c[ p−1] (x, s)φ (31) Ω(s)

Eq. (31) actually states that the reproducing kernel smoothed gradient u r,i should be able to exactly reproduce any ( p − 1)th order polynomials spanned by the monomial basis vector p[ p−1] , which is one order lower than the original monomial basis vector p[ p] used for the field variable approximation. Substituting Eq. (30) into Eq. (31) immediately gives the solution of the unknown vector ˜ a: ˜ a(x) = G −1 (x) p[ p−1] (x) where G(x) is the RKGS moment matrix: ∫ ˜(x, s)dΩ G(x) = p[ p−1] (s) p[ p−1]T (s)φ Ω(s)

(32)

(33)

634

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 1. Illustration of meshfree integration cells.

Now bringing Eq. (32) back into Eq. (29) yields the proposed reproducing kernel smoothed gradient expression: ∫ ˜(x, s)u ,i (s)dΩ ˜ u r,i (x) = p[ p−1]T (x)G −1 (x) p[ p−1] (s)φ (34) Ω(s)

Next, we proceed to evaluate Eqs. (33) and (34) in the context of meshfree numerical integration. Herein the triangular or tetrahedral integration cells with meshfree nodes as cell vertices are employed for the numerical integration of Galerkin meshfree methods. For clarity, the 2D meshfree integration cells are illustrated in Fig. 1, NC NC where the problem domain Ω is partitioned into a set of integration cells {ΩC }C=1 , i.e., ∪C=1 ΩC = Ω , N C represents the number of integration cells. Meanwhile, the following piecewise constant function is selected as the kernel function for gradient smoothing operation [40]: { ˜(x, s) = 1, x ∈ ΩC φ (35) 0, x ∈ / ΩC By the kernel function choice of Eq. (35), the gradient smoothing moment matrix G(x) of Eq. (33) would become a constant matrix within each integration cell ΩC , and for convenience of subsequent development it is then denoted by G C : ∫ GC = p[ p−1] (s) p[ p−1]T (s)dΩ , x ∈ ΩC (36) ΩC

It is noted that when the triangular or tetrahedral integration cells are employed, G C in Eq. (36) can be easily analytically evaluated. With the aid of Eqs. (35) and (36), when x ∈ ΩC , ˜ u r,i in Eq. (34) becomes: ∫ ˜ u r,i (x) = p[ p−1]T (x)G C−1 p[ p−1] (s)u ,i (s)dΩ Ω ∫ C ∫ (37) [ p−1] = p[ p−1]T (x)G C−1 [ p[ p−1] (s)n i u(s)dΓ − p,i (s)u(s)dΩ ] ΓC

ΩC

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

635

where the operation of integration by parts is used. Further introducing the meshfree approximation of Eq. (9) into Eq. (37) yields: ∫ ∫ NP NP ∑ ∑ [ p−1] p[ p−1] (s)n i ˜ u ,ih (x) = p[ p−1]T (x)G C−1 [ Ψ I (s)d I dΓ − p,i (s) Ψ I (s)d I dΩ ] ΓC

= =

NP ∑ I =1 NP ∑

p[ p−1]T (x)G C−1 [

ΩC

I =1

∫

p[ p−1] (s)n i Ψ I (s)dΓ −

ΓC

∫ ΩC

[ p−1]

p,i

I =1

(s)Ψ I (s)dΩ ]d I

(38)

˜I,i (x)d I Ψ

I =1

with ˜I,i (x) = p[ p−1]T (x)G −1 g iCI , x ∈ ΩC Ψ C ∫ ∫ [ p−1] g iCI = p[ p−1] (s)n i Ψ I (s)dΓ − p,i (s)Ψ I (s)dΩ ΓC

(39) (40)

ΩC

˜I,i (x) is the proposed reproducing kernel smoothed gradient. g C means that g i I in Eq. (28) is evaluated where Ψ iI ∑N C C in the integration cell ΩC , and it clear that by definition we have the identity of C=1 gi I = gi I . ˜ It turns out that the proposed reproducing kernel smoothed gradient Ψ I,i (x) defined by Eq. (39) identically meets the integration constraint described by Eq. (27). This fact can be straightforwardly proved by substituting Eq. (39) into Eq. (27): ∫ NC ∫ ∑ ˜I,i (x) p[ p−1] (x)dΩ = ˜I,i (x)dΩ Ψ p[ p−1] (x)Ψ Ω

= =

C=1 ΩC NC ∫ ∑ C=1 ΩC NC ∫ ∑ C=1  ΩC

=

NC ∑

p[ p−1] (x) p[ p−1]T (x)G C−1 g iCI dΩ p[ p−1] (x) p[ p−1]T dΩ G C−1 g iCI  

(41)

GC

g iCI

C=1

= gi I where the definition of g i I in Eq. (28) is utilized. Eq. (41) exactly is the integration constraint and thus the ˜I,i (x) has the inherent integration consistency. It is noted that proposed reproducing kernel smoothed gradient Ψ this inherent integration consistency does not rely on any specific integration algorithms. Moreover, the derivation of the integration constraint in Eq. (25) utilizes the gradient reproducing condition of Eq. (13), consequently the ˜I,i (x) should also meet this condition: reproducing kernel smoothed gradient Ψ NP ∑

˜I,i (x) p[ p] (x I ) = p[ p] (x) Ψ ,i

(42)

I =1

Eqs. (39) and (42) will serve a foundation herein to develop efficient explicit integration algorithms for meshfree methods. 3.2. Meshfree numerical integration with RKGS It is noted that RKGS has the attractive property of inherent integration consistency. Meanwhile, it is also noted ˜I,i is in fact a ( p − 1)th order polynomial and the that in Eq. (39), G C and g iCI are constants within ΩC , thus Ψ normal low order Gauss quadrature rules used in the finite element analysis can be employed to efficiently perform

636

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 2. Barycentric coordinates and integration cell mapping.

the meshfree numerical integration. Accordingly, the stiffness matrix of Eq. (18) can be computed using Gauss integration as follows: ∫ T ˜ B J dΩ B I D˜ KI J = Ω NC ∫ ∑ T ˜ B J dΩ = B I D˜ (43) Ω =

C C=1 NC ∑ NG ∑

T ˜ B I (x G ) D ˜ B J (x G )J (x G )ϖG

C=1 G=1

where x G ’s and ϖG ’s are the standard Gauss integration points and weights, for instance, the 3-point and 4-point Gauss quadrature rules can be used for 2D and 3D meshfree methods with quadratic basis functions [37,38], respectively. J is the determinant corresponding to the geometry mapping between the integration cell ΩC to ξ ˜I,i (x G ), can be the reference parametric triangular or tetrahedral cell ΩC . The component of ˜ B I (x G ), namely, Ψ conveniently calculated using Eq. (39). Besides the stiffness integration, another concern is the numerical integration for the force term [75]. While from the derivation of the integration constraint described by Eqs. (25)–(28), it can be seen that the force term is directly associated with g i I that also serves as a foundation of RKGS as shown in Eqs. (39)–(41). Consequently, if the force term employs the same integration strategy as that used for g i I in RKGS, the integration consistency is automatically fulfilled. Now the left problem is how to compute G C and g iCI given by Eqs. (36) and (40) in an efficient way. 4. Explicit quadrature rules for RKGS From Eqs. (36) and (40), it is clear that G C and g iCI involve the domain and boundary integrations. In this section, a set of explicit quadrature rules consisting of sampling points and corresponding weights are developed for efficient computation of G C and g iCI used in the proposed RKGS. As mentioned earlier, the 2D triangular and 3D tetrahedral integration cells enable an analytical calculation of the RKGS moment matrix G C , and thus these types of simplex integration cells are adopted for the subsequent development. Here we would like to emphasize that the meshfree stiffness matrix is computed using the standard low order Gauss integration rules as shown in Eq. (43), and the quadrature rules developed in this section are for the efficient evaluation of reproducing kernel ˜I,i of Eq. (39), which is required to complete the Gauss integration in Eq. (43). smoothed gradient Ψ 4.1. Formulation of RKGS with barycentric coordinates When simplex integration cells such as 2D triangular and 3D tetrahedral cells are adopted, it is very convenient to use barycentric coordinates [81] to formulate the quadrature rules. As shown in Fig. 2, an integration cell ΩC in ξ the physical space is denoted by ΩC in the parametric space. Accordingly, the relationship between the Cartesian

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

637

coordinate x = (x1 , . . . , xn sd ) and barycentric coordinate ξ = (ξ1 , . . . , ξn sd +1 ) for a generic integration cell ΩC is given by: x = x 1 ξ1 + x 2 ξ2 + · · · + x n sd +1 ξn sd +1

(44)

ξ1 + ξ2 + · · · + ξn sd +1 = 1

(45)

where n sd denotes the number of spatial dimensions. x M denotes the global coordinates of the Mth vertex of the integration cell ΩC . [ p] Meanwhile, the basis vector p[ p] (x) and its gradient p,i (x) can also be expressed by the monomial basis vector [ p] q with barycentric coordinates: p[ p] (x) = A[ p] q [ p] (ξ )

(46) [ p]

[ p]

[ p]

p,i (x) = A[ p] q , j (ξ )J ji−1 = A[ p] B j q [ p−1] (ξ )J ji−1

(47)

[ p]

where A[ p] and B j are some constant coefficient matrices. Ji j is the component of the Jacobi matrix J relating the ξ physical domain ΩC to the parametric domain ΩC . The details for the Jacobi matrices for triangular and tetrahedral [ p] integration cells are illustrated in Appendix A. q is the pth order monomial basis vector that has a similar structure as p[ p] : { }T j p , 0≤i + j +k ≤ p (48) q [ p] (ξ ) = 1, ξ1 , ξ2 , ξ3 , . . . , ξ12 , ξ22 , ξ32 , . . . , ξ1i ξ2 ξ3k , . . . , ξ3 With the notations in Eqs. (46) and (47), G C and g iCI in Eqs. (36) and (40) can be rewritten as: G C = A[ p−1] G C A[ p−1]T

(49)

g iCI = A[ p−1] g iCI

(50)

with ∫

q [ p−1] (ξ )q [ p−1]T (ξ )dΩ

GC = g iCI =

∫ ΩC

Ψ I (ξ )q [ p−1] (ξ )n i dΓ −

(51) ∫

[ p−1]

ΩC

ΓC

Ψ I (ξ )q , j

(ξ )J ji−1 dΩ

(52)

˜I,i of Eq. (39) can Consequently, within a generic integration cell ΩC the reproducing kernel smoothed gradient Ψ be constructed using the barycentric coordinates as: ˜I,i (x) = p[ p−1]T (x)G −1 g C Ψ C iI = q [ p−1]T (ξ )A[ p−1]T (A[ p−1] G C A[ p−1]T )−1 A[ p−1] g iCI −1 = q [ p−1]T (ξ )A[ p−1]T (A[ p−1]T )−1 G C (A[ p−1] )−1 A[ p−1] g iCI       (53) =I

=I

−1

= q [ p−1]T (ξ )G C g iCI ˜I,i (ξ ) = Ψ where I is the identity matrix. For the simplex integration cells used herein, G C in Eq. (51) can be analytically calculated [82]. In what follows, the computation of g iCI given by Eq. (52) by certain quadrature rules is discussed. ˜I,i identically meets the integration constraint provided It is noted that the reproducing kernel smoothed gradient Ψ that the gradient reproducing condition of Eq. (42) is fulfilled. Thus, substituting Eq. (53) into Eq. (42) yields: NP ∑

˜I,i (x) p[ p] (x I ) − p[ p] (x) = Ψ ,i

I =1

NP ∑

˜I,i (ξ ) p[ p] (x I ) − p[ p] (x) Ψ ,i

I =1

=

NP ∑

I =1 NP ∑

= [

I =1

= 0

−1

[ p]

p[ p] (x I )g iCIT G C q [ p−1] (ξ ) − A[ p] B j J ji−1 q [n−1] (ξ ) −1

[ p]

p[ p] (x I )g iCIT G C − A[ p] B j J ji−1 ]q [n−1] (ξ )

(54)

638

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

or NP ∑

[ p]

p[ p] (x I )g iCIT = A[ p] B, j J ji−1 G C

(55)

I =1

Eq. (55) actually is the condition to determine proper quadrature rules for the efficient computation of g iCI . With the aid of Eqs. (47) and (51), the right hand side of Eq. (55) can be rewritten as: ∫ [ p] [ p] B j q [ p−1] (ξ )q [ p−1]T (ξ )dΩ A[ p] B j J ji−1 G C = A[ p] J ji−1    ΩC q

1 = −A D ji n sd VC [ p] = A[ p] G j D ji [ p]

[ p]

∫, j

(ξ )

ΩC

(56)

[ p]

q , j (ξ )q [ p−1]T (ξ )dΩ

with [ p]

Gi

=−

1 n sd VC

∫

[ p]

ΩC

q ,i (ξ )q [ p−1]T (ξ )dΩ

(57)

where Di j is a parameter related to the integration cell geometry, which is given by Eqs. (A.6)–(A.8) in Appendix A. VC represents the length, area and volume of 1D, 2D and 3D integration cells, which are also defined by Eqs. (A.1)–(A.3) in Appendix A. Moreover, introducing certain numerical integration rules into the left hand side of Eq. (55) yields: ∫ ∑ NP NP ∑ CT [ p] = Ψ I (x) p[ p] (x I )n i q [ p−1]T (ξ )dΓ p (x I )g i I ΓC I =1

I =1





p[ p] (x)

NP ∑

∫ −

ΩC I =1

[ p−1]T

Ψ I (x) p[ p] (x I ) J ji−1 q , j



=

∑





L∈I

= A (

(58)



p[ p] (x) [ p] p (x L )n i q [ p−1]T (ξ L )ϖGb SC

[ p]

(ξ )dΩ

−

∑ L∈I

∑

δˆ j(n sd +1) q [ p] (ξ L )q [ p−1]T (ξ L )ϖ Lb +

L∈I [ p]

= A[ p] (R j + T

[ p] j )D ji

[ p−1]T

p[ p] (x L )J ji−1 q , j

(ξ L )ϖ Li VC

1 ∑ [ p] [ p−1]T q (ξ L )q , j (ξ L )ϖ Li )D ji n sd L∈I

with ⎧ [ p] ∑ R = δˆ j(n sd +1) q [ p] (ξ L )q [ p−1]T (ξ L )ϖ Lb ⎪ ⎪ ⎨ i L∈I 1 ∑ [ p] [ p−1]T [ p] ⎪ ⎪ T = q (ξ L )q , j (ξ L )ϖ Li ⎩ i n

(59)

sd L∈I

where I represents the group of integration sample points. ϖ Lb , ϖ Li are the normalized weights for the boundary and domain integrations, respectively. δˆ j(n sd +1) is a constant defined in Appendix A. SC denotes the area or length of ΓC . A combination of Eqs. (55), (56) and (58) gives the following equation which governs the quadrature rules used for the construction of reproducing kernel smoothed gradient: [ p]

Ri

[ p]

+Ti

[ p]

= Gi

Again it is noted that

[ p] Gi

(60) in Eq. (57) can be analytically obtained [82].

4.2. Explicit quadrature rules for RKGS evaluation In order to extract proper and efficient quadrature rules for the numerical integration used in RKGS, herein a twostep algorithm is introduced: (1) Determining the optimal number of integration sample points through minimizing

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

639

Fig. 3. Illustration of 1D sample point types.

Fig. 4. Illustration of 2D sample point types.

the number of sample points from a global point of view; (2) Computing the weights and locations of integration sample points through solving the system of equations deduced from Eq. (60). Another common guideline adopted here is that a boundary integration only involves the boundary sample points of integration cell. To improve the computational efficiency, a minimum number of integration sample points from a global point of view are very desirable. Thus the method of sample point grouping [83,84] is considered in this study. According to this method, a group of triangular or tetrahedral integration cells are packed together and the smallest equivalent number of sample points are searched in order to efficiently perform the numerical integration for RKGS. For convenience of the subsequent development, we define a representative sample point ξˆ = ⟨ξa , ξb , ξc , ξd ⟩ with 0 ≤ ξd ≤ ξc ≤ ξb ≤ ξa ≤ 1, which is a collection of a series of points whose coordinates can be obtained by only permuting the coordinate component sequence. ⟨•⟩ denotes the permutation operator that gives all possible distinct combinations of the coordinate components. Taking into the symmetry consideration, all the points represented by a single ξˆ have the same integration weight. For example, in 2D case, a representative sample point ξˆ = ⟨ξa , ξb , ξc ⟩ implies the following sample points: ⎧ for ξa = ξb = ξc ⎪ ⎪⟨ξa , ξb , ξc ⟩ : (ξa , ξb , ξc ) ⎨ ⟨ξa , ξb , ξc ⟩ : (ξa , ξb , ξc ), (ξb , ξa , ξc ), (ξb , ξc , ξa ) for ξa ̸= ξb = ξc (61) ⟨ξa , ξb , ξc ⟩ : (ξa , ξb , ξc ), (ξc , ξa , ξb ), (ξb , ξc , ξa ), ⎪ ⎪ for ξa ̸= ξb ̸= ξc ⎩ (ξa , ξc , ξb ), (ξc , ξb , ξa ), (ξb , ξa , ξc ) Moreover, for a given function f , ⟨ f (ξˆ )⟩ is defined as a summation of function values induced by all possible coordinates corresponding to ξˆ , for instance, ⟨ f (ξa , ξb , ξc )⟩, ξa ̸= ξb = ξc , means: ⟨ f (ξa , ξb , ξc )⟩ = f (ξa , ξb , ξc ) + f (ξb , ξa , ξc ) + f (ξb , ξc , ξa )

(62)

Meanwhile, as shown in Figs. 3–5 and Tables 1–3, based upon the symmetry requirement and the preference of sharing by neighboring integration cells all possible integration sample points are classified into several types, denoted by S A ’s, according to their locations. Without loss of generality, we further interpret the meaning of

640

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672 Table 1 Properties for different types of 1D sample points. Sample point type

Representative point ξˆ

rA

Unknown

S1

⟨1/2, 1/2⟩

1 = 1/1

ϖ Si

S2

⟨1, 0⟩

1 = 2/2

ϖ Sb , ϖ Si

S3

⟨ξa , ξb ⟩

2 = 2/1

ϖ Si , ξa

Table 2 Properties for different types of 2D sample points. Sample point type Representative point ξˆ

rA

Unknown

S1

⟨1/3, 1/3, 1/3⟩

1 = 1/1

ϖ Si

S2

⟨1, 0, 0⟩

0.5 = 3/6

ϖ Sb , ϖ Si

S3

⟨1/2, 1/2, 0⟩

1.5 = 3/2

ϖ Sb , ϖ Si

S4

⟨ξa , ξa , ξb ⟩ or ⟨ξa , ξb , ξb ⟩

3 = 3/1

ϖ Si , ξa

S5

⟨ξa , ξb , 0⟩

3 = 6/2

ϖ Sb , ϖ Si , ξa

S6

⟨ξa , ξb , ξc ⟩

6 = 6/1

ϖ Si , ξa , ξb

Table 3 Properties for different types of 3D sample points. Sample point typeRepresentative point ξˆ

rA

Unknown

S1

⟨1/4, 1/4, 1/4, 1/4⟩

1 = 1/1

ϖ Si

S2

⟨1, 0, 0, 0⟩

0.2 = 4/20 ϖ Sb , ϖ Si

S3

⟨1/3, 1/3, 1/3, 0⟩

2 = 4/2

ϖ Sb , ϖ Si

S4

⟨ξa , ξa , ξa , ξb ⟩ or ⟨ξa , ξb , ξb , ξb ⟩

4 = 4/1

ϖ Si , ξa

S5

⟨1/2, 1/2, 0, 0⟩

1.2 = 6/5 ϖ Sb , ϖ Si

S6

⟨ξa , ξa , ξb , ξb ⟩

6 = 6/1

S7

⟨ξa , ξb , 0, 0⟩

2.4 = 12/5 ϖ Sb , ϖ Si , ξa

S8

⟨ξa , ξa , ξb , 0⟩ or ⟨ξa , ξb , ξb , 0⟩

6 = 12/2

S9

⟨ξa , ξa , ξb , ξc ⟩ or ⟨ξa , ξb , ξb , ξc ⟩ or ⟨ξa , ξb , ξc , ξc ⟩12 = 12/1 ϖ Si , ξa , ξb

S10

⟨ξa , ξb , ξc , 0⟩

12 = 24/2 ϖ Sb , ϖ Si , ξa , ξb

S11

⟨ξa , ξb , ξc , ξd ⟩

24 = 24/1 ϖ Si , ξa , ξb , ξc

ϖ Si , ξa ϖ Sb , ϖ Si , ξa

S A ’s using the 2D triangular integration cells and sample points. As shown in Fig. 4, S1 denotes the centroid of integration cell and its location is known as (1/3, 1/3, 1/3) and the only unknown is the integration weight. In S1 , the representative sample point coincides with the only sample point, which is not shared by any other integration cells. S2 represents the vertices of integration cell. Clearly, in S2 there is one representative sample point ξˆ containing these sample points, whose locations are known and two integration weights, ϖ Sb for boundary integration and ϖ Si for domain integration, are unknowns. Each point in S2 for the triangular integration cell is shared by six neighboring cells for a generic patch of cells in an average sense [85]. S3 stands for the middle points of three sides of integration cell, which are shared by two adjacent cells. Thus S3 also contains one representative sample point on behalf of three sample points, for which there are two unknown integration weights. S4 implies the interior sample points which lie on one of the axes, say, only one coordinate component and the corresponding weight are unknowns. Clearly, by virtue of symmetry every ξˆ means three points in S4 which are not shared by other cells. Similarly, S5 is related to the boundary sample points with one barycentric coordinate component being zero, and thus only one coordinate component and two integration weights are left to be unknowns. The points in S5 are shared by two adjacent integration cells. The last type of sample points with two coordinate components and one weight as unknowns is indexed by S6 , which are inside the integration cell. According to the symmetry requirement, there are six points for one representative sample point ξˆ in S5 and S6 . Furthermore, it is noted that

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

641

Fig. 5. Illustration of 3D sample point types.

the numbers of representative sample points in S4 , S5 and S6 are unknowns as well. Along the same routine, the S A ’s for 1D and 3D sample points are illustrated in Figs. 3–5 and Tables 1–3, respectively. With the aid of Eqs. (61) and (62), owing to the symmetry property of sampling points in an integration cell ΩC , any integral in Eq. (60) can be concisely rewritten as: ∫ ∑ ∑ f (ξ )dΩ = ϖ L f (ξ L ) = ϖ S ⟨ f (ξˆ S )⟩ (63) ΩC

L∈I

S∈S A1 −A N

where ξˆ S and ϖ S stand for the representative sample point and its corresponding integration weight, S A1 −A N = S A1 ∪ S A2 ∪ · · · ∪ S A N , A N denotes the total number of point types, which equals 6 for the 2D example. Next, for clarity, we present the detailed 2D formulation with quadratic basis function in order to clearly illustrate the derivation of governing equations for the determination of integration points and weights. In 2D case, we have p = 2 and Eq. (60) leads to the following conditions: ∑

1 1 1 1 1 1 1 1 1 ϖ Lb δˆi3 q [3] (ξ L ) = −{0, δi1 , δi2 , δi1 , , δi2 , δi1 , , , δi2 }T 2 2 3 6 3 4 12 12 4 L∈I

(64)

1 1 1 1 1 ϖ Li q [2] (ξ L ) = {1, , , , , }T 3 3 6 12 6 L∈I

(65)

∑

where δi j is Kronecker delta, q [3] is the cubic basis vector referring to the barycentric coordinate system: q [3] (ξ ) = {1, ξ1 , ξ2 , ξ12 , ξ1 ξ2 , ξ22 , ξ13 , ξ12 ξ2 , ξ1 ξ22 , ξ23 }T

(66)

642

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

For instance, the second component in Eq. (64) implies: ∑ ∑ ϖ Lb ( δˆ1 ξ1 − δˆ3 ξ1 ) ϖ Lb δˆ13 ξ1 =  L∈I L∈I ∑ =0 =− ϖ Lb δˆ3 (ξ L )ξ1 (ξ L ) L∈I ∑ =− ϖ Sb ⟨δˆ3 (ξ S )ξ1 (ξ S )⟩ S∈S1−6

∑

=−

ϖ Sb ⟨δˆ3 (ξa , ξb , ξc )ξ1 (ξa , ξb , ξc )⟩

S∈S1−6

∑

=−

S∈S1

=0 ϖ Sb (δˆb ξa

∑

−

ϖ Sb δˆc ξa − 

∑

ϖ Sb (δˆb ξa + δˆc ξb + δˆa ξc )

(67)

S∈S2−4

+ δˆa ξb + δˆc ξb + δˆb ξc + δˆc ξa + δˆa ξc )

S∈S5,6

∑

=−

S∈S2

ϖ Sb ξa −  =1

∑ S∈S3

ϖ Sb ξa −  =1/2

∑ S∈S5

ϖ Sb (ξa + ξb )    =1

1 ∑ 1 =− ϖ Sb ⟨δˆ3 ⟩ = − 2 S∈S 2 2,3,5

or ∑

ϖ Sb ⟨δˆ3 ⟩ = 1

(68)

S∈S2,3,5

Moreover, the fourth component in Eq. (64) gives: ∑ ∑ ∑ ϖ Li ξ12 = ϖ Li ξ12 (ξ L ) = ϖ Si ⟨ξ12 (ξ S )⟩ L∈I

L∈I

∑

=

S∈S1−6 i 2 ϖ S ⟨ξ1 (ξa , ξb , ξc )⟩

S∈S1−6

=

∑ S∈S1

=

∑

ϖ Si ξa2 +  =(1/3)2

ϖ Si

S∈S1

1 + 9 S∈S

∑ S∈S2−4

∑ 2−4

∑

ϖ Si (ξa2 + ξb2 + ξc2 ) + 2

∑

ϖ Si (ξa2 + ξb2 + ξc2 )

S∈S5,6

ϖ Si [(ξa + ξb + ξc )2 − 2(ξa ξb + ξa ξc + ξb ξc )]    =1

+ ξb + ξc )2 − 2(ξa ξb + ξa ξc + ξb ξc )]   S∈S5,6 =1 ∑ ∑ 1 ∑ ϖ Si ξa ξb − 2 ϖ Si (ξa ξb + ξa ξc + ξb ξc ) = ϖ Si ⟨1⟩ − 2 3 S∈S S∈S3,4 S∈S1 ∑1−6 i − 2 ϖ S (ξa ξb + ξb ξa + ξa ξc + ξc ξa + ξb ξc + ξc ξb )

+ 2

ϖ Si [(ξa 

(69)

S∈S5,6

∑ 1 ∑ 1 = ϖ Si ⟨1⟩ − 2 ϖ Si ⟨ξa ξb ⟩ = 3 S∈S 6 S∈S1,3−6 1−6    =1

namely ∑ S∈S1,3−6

ϖ Si ⟨ξ1 ξ2 ⟩ =

1 12

(70)

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

643

Following the same way, it can be shown that along with Eqs. (68) and (70) there are totally four independent conditions arising from Eq. (60) for this 2D case. The other two independent equations read: ∑ 1 (71) ϖ Sb ⟨ξ1 ξ2 δˆ3 ⟩ = 6 S∈S3,5 ∑ ϖ Si ⟨1⟩ = 1 (72) S∈S1−6

Now the four conditions as defined by Eqs. (68)–(72) are used to determine the four unknowns corresponding to RKGS sample points and integration weights. To properly count the number of sample points globally, the following equivalent number of points r A is introduced for each representative sample point ξˆ in S A : number of points represented by ξˆ in S A (73) number of cells sharing a typical point in S A Accordingly, the r A ’s for various 2D and 3D cases are tabulated in Tables 1–3. The equivalent number of points r A will be used to minimize the total number of sample points. Subsequently, we would like to minimize the number of sample points from a global point of view as follows: ∑ argmin r Am A (74) rA =

A

s.t. c(m) ≥ 0

(75)

where m A stands for the minimum number of representative sample points in S A , m is a vector containing m ′A s, c(m) represents a set of constraints for m. It is noted that the minimization problem described by Eqs. (74) and (75) is a standard problem of linear and integer programming, and it can be effectively solved by many well established schemes, for example, the branch and bound algorithm [86]. Here we still use the previous 2D example to elaborate the constraints defined by Eq. (75). According to the definition of m A , we have: m1, m2, m3, m4, m5, m6 ≥ 0

(76)

m1, m2, m3 ≤ 1

(77)

Moreover, the possible total number of unknowns for all sets used in Eqs. (68)–(72) should not be less than the number of unknown sampling points and weights: ∑ ϑ A m A = m 1 + 2m 2 + 2m 3 + 2m 4 + 3m 5 + 3m 6 ≥ 4 (78) A

where ϑ A denotes the total number of unknowns corresponding to a representative point in S A , which can also be found from Tables 1–3. However, it is noted that these unknowns may not be always active. From Eqs. (68) and (71), it can be seen the numbers of active unknowns of S2 , S3 and S5 in these two equations are 1, 1 and 2, respectively. Thus the total numbers of active unknowns should also be greater than or equal to the number of associative equations, namely, 2, under this circumstance: m 2 + m 3 + 2m 5 ≥ 2

(79)

Similarly, the numbers of active unknowns of S1 and S3 −S6 in Eqs. (70) and (71) are 1, 2, 2, 3 and 3; the numbers of active unknowns of S3 and S5 in Eq. (71) are 1 and 2, which yield the following two additional constraints: m 1 + 2m 3 + 2m 4 + 3m 5 + 3m 6 ≥ 2

(80)

m 3 + 2m 5 ≥ 2

(81)

At this point, Eqs. (76)–(81) constitute the constraint conditions corresponding to Eq. (75). Following the same procedure, the independent governed equations of sample points and integration weights and the corresponding constraints for other 1D, 2D and 3D cases are provided in Appendix B.

644

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672 Table 4 1D quadrature rules for RKGS. Basis order

Figure

p=2

p=3

Point type

Coordinate

ϖ Sb

ϖ Si

S1

⟨ 12 , 12 ⟩

0

2 3

S2

⟨1, 0⟩

1

S2

⟨1, 0⟩

1

1 6 1 12

S3

⟨ 5+10 5 ,

0

5 12

√

√ 5− 5 10 ⟩

Thereafter, based upon the standard branch and bound algorithm [86], the optimal number of point sets can be obtained, for example, the result for 2D meshfree formulation with quadratic basis function is: { }T ∑ m= 0 1 1 0 0 0 , r Am A = 2 (82) A

As shown in Eq. (82), the optimal sample point types for 2D quadratic case are S2 and S3 , where each type contains one representative point with three sample points. From a global point of view the average number of sample points per integration cell is 2. The results of optimal sample points for other cases are described in Tables 1–3. Up to now, the optimal number of integration sample points has been achieved. Subsequently, the corresponding locations and weights can be attained through solving the associated governing equations, for example, Eqs. (68)–(72) for the 2D quadratic case. Substituting the unknowns defined by m in Eq. (82) into Eqs. (68)–(72) leads to: ⎧ b 2ϖS2 + ϖSb3 = 1 ⎪ ⎪ ⎪ ⎪ 1 1 ⎪ ⎨ ϖSb = 4 i 3 6 (83) 3ϖS2 + 3ϖSi 3 = 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎩1ϖi = 1 4 S3 12 Consequently, the following four integration weights are obtained: 2 1 1 (84) ϖSb2 = , ϖSb3 = , ϖSi 2 = 0, ϖSi 3 = 6 3 3 Finally, following exactly the same path, explicit quadrature rules with a minimum number of sample points can be established to efficiently carry out the RKGS computation for 1D, 2D and 3D meshfree methods with quadratic and cubic basis functions, respectively. The corresponding results for sample points and integration weights are detailed in Tables 4–6. 5. Numerical examples In this section, a series of numerical tests ranging from 1D to 3D are performed to verify the proposed methodology. Both quadratic and cubic basis functions are employed for the reproducing kernel meshfree approximation. For the patch tests and potential problems, ki j = δi j is used for in Eq. (14) and the following general form of analytical solutions is considered: ⎧ 1D ⎨{(i, j, k)|i = n, j = k = 0} ∑ ( j + 1)(k + 2) u= x i y j z k , C = {(i, j, k)|i + j = n, j ≤ i, k = 0} 2D (85) ⎩ 2 {(i, j, k)|i + j + k = n, j, k ≤ i} 3D (i, j,k)∈C where n = 1, 2, 3 are utilized in the linear, quadratic and cubic patch tests, respectively. For comparison purpose, L 2 and H1 error norms are calculated for all examples. In the numerical results, “GI-k” means using k-point Gauss quadrature rule in each integration cell, NSGSI stands for the nesting sub-domain gradient smoothing integration method [75] which has a property of quadratic

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

645

Table 5 2D quadrature rules for RKGS. Basis order p=2

Figure

Point type

Coordinate

ϖ Sb

ϖ Si

S2

⟨1, 0, 0⟩

1 6

0

S3

⟨ 21 , 12 , 0⟩

2 3

1 3

S1

⟨ 31 , 13 , 13 ⟩

0

9 20

S2

⟨1, 0, 0⟩

1 20

1 − 36

S3

⟨ 21 , 12 , 0⟩

16 45

4 135

S5

⟨ 7+14 21 ,

49 180

49 540

p=3

√

√ 7− 21 14 , 0⟩

exactness, and RKGSI represents the present reproducing kernel gradient smoothing integration method. Meanwhile, for convenience of comparison, similar to the conventional finite element analysis, here we also define the integration rules as normal ones if they can exactly integrate the stiffness matrices referring to the triangular cells when the shape functions are pth order polynomials. Accordingly, GI-2, GI-3 and GI-4 are the normal integration rules for 1D, 2D and 3D meshfree formulations with quadratic basis functions; GI-3, GI-6 and GI-11 are the normal rules for 1D, 2D and 3D meshfree formulations with cubic basis functions [38]. Of course, due to the rational nature of meshfree shape functions, the convergence rates by meshfree formulation using these so-called normal integration rules are usually far below than the optimal values. Higher order integration rules with much more Gauss points per cell will also be investigated in the numerical tests to compare and demonstrate the strength of the proposed RKGSI. On the other hand, we would like to emphasize that in the present RKGSI, as shown in Eq. (43), the stiffness matrix is evaluated using the normal Gauss integration rules corresponding to its basis order, but rather than using the direct shape function derivatives, the proposed reproducing kernel smoothed gradients are adopted at Gauss points, which can be efficiently computed within the RKGS framework as described in Section 4. 5.1. Patch tests (1) 1D patch tests In 1D patch tests regarding the potential problem, the problem domain is Ω = (0, 1) and the analytical solutions of linear, quadratic and cubic patch tests are defined by Eq. (85) with j = k = 0 and n = 1, 2, 3. The essential and natural types of boundary conditions following the analytical solutions are applied at the left and right end points. The essential boundary condition is enforced via the Nitsche’s method [87] throughout this study. As shown in Fig. 6, irregular discretizations with 11 and 21 nodes are adopted for the 1D patch tests with quadratic and cubic meshfree formulations. Normalized support sizes of 2.5 and 3.5 are employed for quadratic and cubic meshfree approximations, respectively. Here, for quadratic meshfree formulation, linear and quadratic patch tests are performed, while for cubic meshfree formulation, quadratic as well as cubic patch tests are studied. The 1D patch test results are presented in Tables 7 and 8, which indicate noticeable errors still exist even for higher order Gauss quadrature rules such as GI-5 for quadratic formulation and GI-9 for cubic formulation, and the errors are quite large when GI-2 and GI-3 are used for quadratic and cubic meshfree methods. NSGSI can pass the linear and

646

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672 Table 6 3D quadrature rules for RKGS. Basis order

Figure

Point type

Coordinate

ϖ Sb

ϖ Si

S2

⟨1, 0, 0, 0⟩

1 20

1 40

S3

⟨ 13 , 13 , 13 , 0⟩

2 15

9 40

S5

⟨ 12 , 12 , 0, 0⟩

9 20

0

S1

⟨ 14 , 14 , 14 , 14 ⟩

0

32 105

S2

⟨1, 0, 0, 0⟩

1 90

1 − 210

S3

⟨ 13 , 13 , 13 , 0⟩

81 320

243 4480

S5

⟨ 12 , 12 , 0, 0⟩

16 225

2 175

S8

⟨ 57 , 17 , 17 , 0⟩

2401 14400

343 9600

p=2

p=3

Table 7 Results of patch tests for 1D quadratic meshfree formulation. Method

GI-2 GI-5 NSGSI RKGSI

Linear patch test

Quadratic patch test

L 2 error

H1 error

L 2 error

H1 error

2.5349E−02 3.5735E−04 7.6444E−15 3.4517E−15

1.5295E−01 2.9346E−03 9.8592E−15 7.4629E−15

3.3460E−02 5.1410E−04 7.8613E−15 4.0030E−15

1.5593E−01 3.0177E−03 8.0852E−15 6.7969E−15

quadratic tests, but not the cubic test. This is because that NSGSI only meets the quadratic integration consistency. On the other hand, the proposed RKGSI can exactly pass the 1D patch tests for all cases, which is further verified by the gradient errors of patch tests depicted in Fig. 6. (2) 2D patch tests The 2D patch tests of potential problem are carried out in a square domain Ω = (0, 1) ⊗ (0, 1), where the top and bottom sides are prescribed by essential boundary conditions, and the left and right sides are subjected to natural boundary conditions. The boundary conditions are consistent with the analytical solutions defined by Eq. (85) with k = 0 and n = 1, 2, 3. The meshfree discretization with 81 nodes as shown in Fig. 7 is used for the 2D patch tests. In the computation, normalized support sizes of 2.5 and 3.5 are adopted for the kernel function of quadratic and cubic meshfree approximations. Tables 9 and 10 list the results of patch tests for 2D quadratic and cubic meshfree formulations, and Figs. 8 and 9 plot the gradient contours of patch tests. All the results well support

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

647

Table 8 Results of patch tests for 1D cubic meshfree formulation. Method

Quadratic patch test

Cubic patch test

L 2 error

H1 error

L 2 error

H1 error

GI-3 GI-6 GI-9 RKGSI

1.4381E−03 1.5557E−04 3.6050E−05 2.4427E−15

1.1767E−02 1.5079E−03 3.4561E−04 4.0235E−14

1.9635E−03 2.5832E−04 5.7308E−05 2.0209E−15

1.2680E−02 1.9159E−03 4.1597E−04 3.1661E−14

Table 9 Results of patch tests for 2D quadratic meshfree formulation. Method

GI-3 GI-13 NSGSI RKGSI

Linear patch test

Quadratic patch test

L 2 error

H1 error

L 2 error

H1 error

3.5116E−03 1.2634E−04 4.4466E−15 3.1692E−15

2.6492E−02 1.0405E−03 2.7518E−14 1.8522E−14

4.6009E−03 1.1399E−04 8.7591E−15 4.4470E−15

2.7591E−02 6.6605E−04 2.5700E−14 2.2258E−14

Table 10 Results of patch tests for 2D cubic meshfree formulation. Method

GI-6 GI-16 RKGSI

Quadratic patch test

Cubic patch test

L 2 error

H1 error

L 2 error

H1 error

2.4856E−03 6.3696E−05 2.9083E−14

1.7120E−02 4.7423E−04 2.3223E−13

4.9528E−03 1.2310E−04 2.3118E−14

2.6165E−02 6.9808E−04 1.5716E−13

Fig. 6. Results of 1D patch tests: (a) quadratic patch test for quadratic meshfree formulation; (b) cubic patch test for cubic meshfree formulation.

that the proposed RKGSI exactly passes various patch tests, in contrast to the obvious errors arising from meshfree formulations with higher order Gauss quadrature rules, say, GI-13 for quadratic formulation and GI-16 for cubic formulation, respectively. Again, it is evident that NSGSI at most passes the quadratic patch test as expected. (3) 3D patch tests For the 3D patch tests, a cube domain Ω = (0, 1) ⊗ (0, 1) ⊗ (0, 1) is considered, which is irregularly discretized by 355 meshfree nodes as shown in Fig. 10. The analytical solutions again are given by Eq. (85) with n = 1, 2, 3.

648

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 7. Irregular meshfree discretization for 2D patch tests.

Fig. 8. Quadratic patch test contour plots of u ,x for 2D quadratic meshfree formulation.

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 9. Cubic patch test contour plots of u ,y for 2D cubic meshfree formulation.

Fig. 10. Irregular meshfree discretization for 3D patch tests.

649

650

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 11. Quadratic patch test contour plots of u ,x for 3D quadratic meshfree formulation.

The essential boundary conditions according to the analytical solutions are enforced on the three planes of x = 0, y = 0 and y = 0, and the other surfaces of the cube are applied by the corresponding natural boundary conditions. Similar to the previous examples, the normalized support sizes of 2.5 and 3.5 are chosen for the quadratic and cubic meshfree approximations. The patch test results including L 2 and H1 errors are tabulated in Tables 11 and 12 and the related gradient contours are portrayed in Figs. 11 and 12. These 3D results once again demonstrate that the standard Gauss integration based meshfree methods cannot exactly pass the patch tests, and at the same time the exactness for patch tests is always ensured by the proposed RKGSI. 5.2. Convergence study for potential problems (1) 1D potential problem In order to investigate the convergence of the proposed method, a 1D example corresponding to Eq. (85) with j = k = 0 and n = 4 is first analyzed. The geometry and material properties are the same as those of 1D patch tests. The boundary conditions at the left and right ends are essential and natural types, respectively, which are deduced from the given analytical solution. The normalized support sizes of kernel function are 2.5 and 3.5 for the quadratic and cubic meshfree shape functions. Three progressively refined meshfree discretizations with 21, 41, 81 nodes are used to perform the convergence assessment. The convergence results of L 2 and H1 error norms are presented in Figs. 13–14. The results clearly show that the conventional Gauss quadrature based meshfree formulations cannot achieve the optimal convergence rates for both quadratic and cubic basis functions, even though higher order integration rules of GI-5 and GI-9 are used. On the contrary, optimal convergence rates of

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

651

Fig. 12. Cubic patch test contour plots of u ,x for 3D cubic meshfree formulation. Table 11 Results of patch tests for 3D quadratic meshfree formulation. Method

Linear patch test

Quadratic patch test

L 2 error

H1 error

L 2 error

H1 error

GI-4 GI-14 QCI RKGSI

7.4095E−03 8.7366E−04 1.0597E−14 7.9957E−15

4.5474E−02 5.4230E−03 9.6312E−14 9.5867E−14

1.0243E−02 1.1954E−03 8.6715E−15 7.6245E−15

5.5566E−02 6.5564E−03 8.2578E−14 7.1838E−14

(p + 1) for the L 2 error and p for the H1 error are produced by the proposed RKGSI. It is also noted that NSGSI performs similarly as RKGSI for quadratic formulation, while it is not applicable to higher order basis functions. (2) 2D potential problem The 2D potential problem also has similar statements as the previous 2D patch test example, but corresponds to an analytical solution of Eq. (85) with k = 0 and n = 4. Fig. 15 shows the three discretizations with 9 × 9, 17 × 17 and 33 × 33 nodes used for this 2D problem. Normalized support sizes of 2.5 and 3.5 are utilized for the kernel function to construct quadratic and cubic meshfree shape functions. The results regarding both convergence and efficiency are listed in Figs. 16–19, where GI-4, GI-14, NSGSI, QCI and RKGSI are adopted for the quadratic meshfree formulation, and GI-6, GI-16 and the present RKGSI are employed the cubic meshfree formulation, respectively.

652

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672 Table 12 Results of patch tests for 3D cubic meshfree formulation. Method

Quadratic patch test

Cubic patch test

L 2 error

H1 error

L 2 error

H1 error

GI-11 GI-20 RKGSI

2.0379E−03 5.2100E−04 3.8580E−14

1.8173E−02 3.7003E−03 4.2027E−13

4.4804E−03 1.0618E−03 5.3535E−14

3.2060E−02 5.8520E−03 3.9724E−13

Fig. 13. Convergence comparison for the 1D potential problem by quadratic meshfree formulation: (a) L 2 error; (b) H1 error.

Fig. 14. Convergence comparison for the 1D potential problem by cubic meshfree formulation: (a) L 2 error; (b) H1 error.

For the efficiency comparisons in Figs. 17 and 19, CPU cost vs. number of particles as well as L 2 error vs. CPU cost are examined. Among the different methods, the convergence results in Figs. 16 and 18 systematically exhibit the optimal convergence behavior of RKGSI. It is also noted that NSGSI and QCI yield similar convergence results as RKGSI for the quadratic meshfree formulation, but the proposed RKGSI is more preferable due to its higher efficiency as shown in Fig. 17. Moreover, in Fig. 19 the efficiency superiority of RKGSI is also well demonstrated for the cubic meshfree formulation. The higher efficiency of RKGSI is because it does not involve the derivative computation whose cost grows dramatically with the elevation of basis orders, Figs. 17 and 19 essentially show that

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

653

Fig. 15. Meshfree discretizations for the 2D potential problem: (a) 81 nodes; (b) 289 nodes; (c) 1089 nodes.

Fig. 16. Convergence comparison for the 2D potential problem by quadratic meshfree formulation: (a) L 2 error; (b) H1 error.

the absolute CPU cost of RKGSI is even less than GI-3 for quadratic formulation and GI-6 for cubic formulation, respectively. (3) 3D potential problem The 3D potential problem defined in a cube domain Ω = (0, 1) ⊗ (0, 1) ⊗ (0, 1), which has an analytical solution of Eq. (85) with n = 4. As for the boundary conditions computed from the analytical solution, the three surfaces of x = 0, y = 0 and y = 0 are subjected to the prescribed field variable values, while the rest three surfaces are under the prescribed boundary fluxes. As shown in Fig. 20, the meshfree discretizations with 5 × 5 × 5, 9 × 9 × 9 and 17 × 17 × 17 nodes are used for the meshfree computation, where normalized support sizes of 2.5 and 3.5 are employed for the quadratic and cubic approximation, respectively. The convergence and efficiency results for this 3D potential problem are summarized in Figs. 21–24, as consistently evince the superior performance of the proposed method of RKGSI regarding convergence rates, accuracy as well as efficiency, compared with the 3D normal Gauss integration rules of GI-4 and GI-11, and the higher order Gauss integration rules of GI-14 and GI-20 for the quadratic and cubic meshfree formulations. 5.3. Elasticity problems In this sub-section, two elasticity benchmark problems are presented to further examine the proposed method of RKGSI. It is noted that the proposed methodology here is directly applicable to elasticity problems, for example,

654

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 17. Efficiency comparison for the 2D potential problem by quadratic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

Fig. 18. Convergence comparison for the 2D potential problem by cubic meshfree formulation: (a) L 2 error; (b) H1 error.

the elasticity stiffness matrix can be obtained by simply replacing ˜ B I and D in Eq. (43) by their standard elasticity counterparts. (1) 2D elastic hollow cylinder problem Consider the 2D hollow elastic cylinder problem as shown in Fig. 25. This hollow cylinder has an inner radius a = 1 and an outer radius b = 2, which is subjected to both inner and outer external pressures of pi = 1000 and po = 500. For this problem, the material constants are Young’s modulus E = 2 × 107 and Poisson’s ratio ν = 0.3. Due to the symmetry, the hoop displacement component u θ = 0 and the analytical solution of this problem for radial displacement component u r is [88]: [ ] 1 a 2 b2 ( pi − po ) u r (r ) = (1 − ν)(a 2 pi − b2 po )r + (1 − ν) (86) r E(b2 − a 2 )

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

655

Fig. 19. Efficiency comparison for the 2D potential problem by cubic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

Fig. 20. Meshfree discretizations for the 3D potential problem: (a) 125 nodes; (b) 729 nodes; (c) 4913 nodes.

Fig. 21. Convergence comparison for the 3D potential problem by quadratic meshfree formulation: (a) L 2 error; (b) H1 error.

656

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 22. Efficiency comparison for the 3D potential problem by quadratic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

Fig. 23. Convergence comparison for the 3D potential problem by cubic meshfree formulation: (a) L 2 error; (b) H1 error.

where r and θ are the standard radial and angular coordinate components. E and ν are given by: ⎧ ⎨ E = E, ν=ν (plane stress) E (87) ν ⎩E = , ν = 1−ν (plane strain) 1 − ν2 For this elasticity problem, we take the advantage of the two-fold symmetry and only model a quarter part of the cylinder. The non-uniform meshfree discretizations with 153, 561 and 2145 nodes are illustrated in Fig. 26. For quadratic and cubic basis functions, the normalized support sizes of kernel function are set to be 2.7 and 3.8, respectively. The results for quadratic meshfree formulation are presented in Figs. 27 and 28, which show that for this problem with non-uniform discretizations, the normal rule of GI-4 even cannot ensure the convergence and the higher order rule like GI-13 still yields pretty low convergence rates. In contrast, optimal rates of convergence can be achieved by the methods of QCI, NSGSI and the present RKGSI, while RKGSI outperforms other methods when the efficiency is taken into account. Figs. 29 and 30 provide the results for cubic meshfree formulation, as also well testify that the present RKGSI performs superiorly compared with both GI-6 and GI-16. Besides, a detailed efficiency comparison for this problem using 2145 nodes is given Fig. 31, where the CPU costs for shape function

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

657

Fig. 24. Efficiency comparison for the 3D potential problem by cubic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

Fig. 25. Description of the hollow cylinder problem.

derivative computation and global stiffness assembly are emphasized. It can be clearly seen that the cost of shape function derivative computation dominates the overall cost. Since the proposed RKGSI does not involve the direct derivative calculation and also optimize the explicit quadrature rules with minimum sample points for reproducing kernel smoothed gradient evaluation, it significantly accelerates the meshfree computation in comparison with other methods as illustrated in Fig. 31. The contour plots for the stress components σx x and σx y in Figs. 32 and 33 further verify the superior performance of RKGSI. (2) 3D hollow sphere problem The last example is the hollow sphere elasticity problem depicted in Fig. 34. The inner and outer radii of the hollow sphere are a = 1 and b = 2. This hollow sphere is under an inner pressure of pi = 1000. For this problem, the material properties are Young’s modulus E = 2 × 107 , Poisson’s ratio ν = 0.3, and its analytical solution for

658

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 26. Meshfree discretizations for the hollow cylinder problem: (a) 153 nodes; (b) 561 nodes; (c) 2145 nodes.

Fig. 27. Convergence comparison for the hollow cylinder problem by quadratic meshfree formulation: (a) L 2 error; (b) H1 error.

Fig. 28. Efficiency comparison for the hollow cylinder problem by quadratic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

659

Fig. 29. Convergence comparison for the hollow cylinder problem by cubic meshfree formulation: (a) L 2 error; (b) H1 error.

Fig. 30. Efficiency comparison for the hollow cylinder problem by cubic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

the only nonzero displacement component u r is given by [88]: { } a3 b2 pi u r (r ) = 3 (1 − 2ν)r + (1 + ν) 2 b − a3 2r E

(88)

Due to the three-fold symmetry, only one-eighth of the hollow sphere is modeled and the corresponding nonuniform meshfree discretizations with 113, 428 and 2170 nodes are shown in Fig. 35. The quadratic and cubic meshfree formulations employ normalized support sizes of 2.1 and 3.1 for the kernel function, respectively. The convergence and efficiency comparison results for this 3D elasticity problem are presented in Figs. 36–39. The convergence results in Figs. 36 and 38 manifest that optimal convergence rates are consistently attained by the proposed RKGSI for both quadratic ( p = 2) and cubic ( p = 3) meshfree approximations, which are 3 and 4 for the L 2 error, and 2 and 3 for the H1 error, respectively. In the meantime, the efficiency results in Figs. 37 and 39 also display that compared with various Gauss integration rules, the computational cost is remarkably reduced by the proposed RKGSI for 3D meshfree formulations. This efficiency advantage of RKGSI is further disclosed by the detailed comparison of CPU costs for meshfree stiffness assembly and derivative computation with 2170 nodes as

660

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 31. Detailed efficiency comparison for the hollow cylinder problem: (a) Quadratic meshfree formulation; (b) Cubic meshfree formulation.

Fig. 32. Contour plots of σx x for the hollow cylinder problem by cubic meshfree formulation.

shown in Fig. 40. Moreover, the stress comparisons in Figs. 41 and 42 once again prove that the proposed RKGSI yields much more accurate results in comparison with both GI-11 and GI-20.

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 33. Contour plots of σx y for the hollow cylinder problem by cubic meshfree formulation.

Fig. 34. Description of the hollow sphere problem.

661

662

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 35. Meshfree discretizations for the hollow sphere problem: (a) 113 nodes; (b) 428 nodes; (c) 2170 nodes.

Fig. 36. Convergence comparison for the hollow sphere problem by quadratic meshfree formulation: (a) L 2 error; (b) H1 error.

Fig. 37. Efficiency comparison for the hollow sphere problem by quadratic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

6. Conclusions A general inherently consistent reproducing kernel gradient smoothing framework with explicit quadrature rules was proposed for efficient Galerkin meshfree analysis. In this framework, the meshfree smoothed gradients were

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

663

Fig. 38. Convergence comparison for the hollow sphere problem by cubic meshfree formulation: (a) L 2 error; (b) H1 error.

Fig. 39. Efficiency comparison for the hollow sphere problem by cubic meshfree formulation: (a) NP vs. CPU time; (b) CPU time vs. L 2 error.

built upon the standard gradients of field variables in the context of reproducing kernel formulation. It was shown that the present reproducing kernel smoothed gradients identically meet the integration constraint required to exactly reproduce any solutions spanned by the meshfree basis vector, no matter what integration schemes are used for Galerkin weak form. The proposed methodology reduces to the stabilized conforming nodal integration when linear basis function is used. The conventional normal low order Gauss integration algorithm now becomes well suitable for the stiffness matrix integration, except that the present reproducing kernel smoothed gradients are used at Gauss points, rather than the standard gradients of meshfree shape functions. The evaluation of reproducing kernel smoothed gradients was realized by a set of new explicit quadrature rules referring to 2D triangular and 3D tetrahedral integration cells. In these reproducing kernel gradient smoothing integration rules, an equivalent number of sample points was introduced to account the sample points from a global point of view and accordingly the total number of sample points was minimized in order to improve the computational efficiency. Particular attention was paid to the higher order basis functions, for example, quadratic and cubic basis functions. It is noted that the proposed explicit reproducing kernel gradient smoothing integration rules allow a very straightforward numerical implementation. A series of numerical examples ranging from 1D and 3D were studied to verify the proposed approach. Numerical results systematically demonstrated that Galerkin meshfree formulations using the present

664

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

Fig. 40. Detailed efficiency comparison for the hollow sphere problem: (a) Quadratic meshfree formulation; (b) Cubic meshfree formulation.

Fig. 41. Hollow sphere problem contour plots of stress σrr for the 3D cubic meshfree approximation.

reproducing kernel gradient smoothing integration exactly pass all patch tests up to the basis order of meshfree approximation. Meanwhile, very favorable performances in terms of convergence, accuracy as well as efficiency were consistently revealed for the proposed reproducing kernel gradient smoothing framework, in contrast to both the normal low order Gauss integration schemes and the costly higher order Gauss quadrature rules. Acknowledgments The support of this work by the National Natural Science Foundation of China (11772280, 11472233) and the Natural Science Foundation of Fujian Province of China (2014J06001) is gratefully acknowledged.

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

665

Fig. 42. Hollow sphere problem contour plots of stress σθθ for the 3D cubic meshfree approximation.

Appendix A In this Appendix, some geometric properties of simplex integration cells are summarized. Here VC represents the length, area and volume of simplex integration cell ΩC for 1D, 2D and 3D cases, respectively. Under different circumstances, VC can be expressed as follows: (1) 1D case [ VC = det

1 x1

1 x2

] (A.1)

(2) 2D case ⎡ 1 1 VC = det ⎣x1 2 y 1

1 x2 y2

⎤ 1 x3 ⎦ y3

⎡ 1 ⎢x1 1 VC = det ⎢ ⎣ y1 6 z1

1 x2 y2 z2

1 x3 y3 z3

(A.2)

(3) 3D case ⎤ 1 x4 ⎥ ⎥ y4 ⎦ z4

(A.3)

where x M = (x M , y M , z M ) denotes the global coordinates of the Mth vertex of the integration cell. Meanwhile, the Jacobi matrix J between Cartesian and barycentric coordinates is defined as: n

J = [Ji j ] = [

sd ∑ ∂ xi (ξ1 , . . . , ξn sd ) ], ξn sd +1 = 1 − ξi , i, j = 1, . . . , n sd ∂ξ j i=1

(A.4)

666

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

For simplex integration cells, the component of the inverse Jacobi matrix J −1 , denoted by Ji−1 j can be analytically computed as: 1 Di j , i, j = 1, . . . , n sd (A.5) Ji−1 j =− n sd VC where Di j is a parameter related to the integration cell geometry and it is given by [89]: (1) 1D case { } { } 1 D11 = −1 D21

(A.6)

(2) 2D case } { } { } { } { } { } { y3 − y2 D21 y1 − y3 D31 y2 − y1 D11 = , = , = x2 − x3 D22 x3 − x1 D32 x1 − x2 D12 (3) 3D case ⎡ ⎧ ⎫ ⎧ ⎫ 2 3 ⎨ Di1 ⎬ 1 ∑ ⎨ yc(i, j) z c(i, j+1) − yc(i, j+1) z c(i, j) ⎬ ⎢1 Di2 = z c(i, j) xc(i, j+1) − z c(i, j+1) xc(i, j) , c = ⎢ ⎣1 ⎩ ⎭ 2 ⎩ ⎭ Di3 xc(i, j) yc(i, j+1) − xc(i, j+1) yc(i, j) j=1 1

3 4 2 2

⎤ 4 3⎥ ⎥ 4⎦ 3

(A.7)

(A.8)

By definition, it is straightforward to show that Dki satisfies the following relationship: n sd +1

∑

Dki = 0

(A.9)

k=1

Moreover, for the outward normal vector n(ξ ) of integration cell surface, its component n i (ξ ) can be expressed via Dki as: 1 n i (ξ ) = δˆk (ξ )Dki , i = 1, . . . , n sd , k = 1, . . . , n sd + 1 (A.10) S(ξ ) or 1 n i (ξ ) = δˆ j(n sd +1) (ξ )D ji i, j = 1, . . . , n sd (A.11) S(ξ ) where δˆk and δˆkl (ξ ) are defined by: { ˆδk (ξ ) = 1 ξk = 0 , δˆkl (ξ ) = δˆk (ξ ) − δˆl (ξ ) (A.12) 0 ξk ̸= 0 S(ξ ) stands for the boundary length or area of simplex integration cell and it takes the following form: n sd ∑ √ S(ξ ) = δˆk Dki Dki

(A.13)

k=1

where S = 1 in 1D case. For example, the red circle point with ξ = (0, 1/2, 1/2) in Fig. 2 is the mid-point of one of the integration cell boundary and thus we have δˆ1 (ξ ) = 1, δˆ2 (ξ ) = 0 and δˆ3 (ξ ) = 0 in accordance with Eq. (A.12). Subsequently the corresponding outward normal vector and length of the specific boundary containing the given point can be computed according to Eqs. (A.11) and (A.13) as follows: n sd +1

S(ξ ) =

∑

√ δˆk Dki Dki

k=1 √ √ √ = δˆ1 D1i D1i + δˆ2 D2i D2i + δˆ3 D3i D3i    =0 =0 √=1 2 2 = D11 + D12 √ = (y3 − y2 )2 + (x2 − x3 )2

(A.14)

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

⎧ ⎪ ⎪ ⎪ ⎪ ⎨n 1 (ξ ) =

1 D11 1 δˆ j(n sd +1) D j1 = (δˆ1 D11 + δˆ2 D21 + δˆ3 D31 ) = S(ξ ) S(ξ ) S(ξ )

⎪ ⎪ ⎪ ⎪ ⎩n 2 (ξ ) =

1 1 D12 δˆ j(n sd +1) D j2 = (δˆ1 D12 + δˆ2 D22 + δˆ3 D32 ) = S(ξ ) S(ξ ) S(ξ )

667

(A.15)

Appendix B In this Appendix, all independent governing equations to determine the sample points and weights for RKGS integration are listed for 1D, 2D and 3D meshfree formulations using quadratic and cubic basis functions: (1) 1D quadratic case ⎧ ∑ ⎪ ⎪ ϖ Sb ⟨δˆ2 ⟩ = 1 ⎪ ⎪ ⎪ ⎪ S∈I S2 ⎪ ⎪ ⎪ ⎪ ⎨ ∑ ϖ Si ⟨1⟩ = 1 ⎪ S∈I S1−3 ⎪ ⎪ ⎪ ⎪ ∑ ⎪ 1 ⎪ ⎪ ϖ Si ⟨ξ1 ξ2 ⟩ = ⎪ ⎪ ⎩ 6 S∈I

(B.1)

S1,3

(2) 1D cubic case ⎧ ∑ ⎪ ⎪ ⎪ ϖ Sb ⟨δˆ2 ⟩ = 1, ⎪ ⎪ ⎨ S∈I S2

∑ 1 ⎪ ⎪ ⎪ ϖ Si ⟨ξ1 ξ2 ⟩ = , ⎪ ⎪ 6 ⎩ S∈I S1,3

(3) 2D quadratic case ⎧ ∑ ⎪ ⎪ ⎪ ϖ Sb ⟨δˆ3 ⟩ = 1, ⎪ ⎪ ⎨ S∈I S2,3,5

∑ ⎪ ⎪ ⎪ ϖ Si ⟨1⟩ = 1, ⎪ ⎪ ⎩ S∈I S1−6

∑ S∈I S1−3

∑

ϖ Si ⟨ξ12 ξ22 ⟩

S∈I S1,3

∑

ϖ Sb ⟨ξ1 ξ2 δˆ3 ⟩ =

∑

ϖ Si ⟨ξ1 ξ2 ⟩ =

S∈I S1,3−6

∑

1 6 (B.3)

1 12

ϖ Sb ⟨ξ1 ξ2 δˆ3 ⟩ =

S∈I S3,5

∑

(B.2)

1 = 30

S∈I S3,5

(4) 2D cubic case ⎧ ∑ ⎪ ϖ Sb ⟨δˆ3 ⟩ = 1, ⎪ ⎪ ⎪ ⎪ ⎪ S∈I S2,3,5 ⎪ ⎪ ⎪ ⎪ ∑ ϖ b ⟨ξ 2 ξ 2 δˆ ⟩ = 1 , ⎪ ⎪ S 1 2 3 ⎪ 30 ⎪ ⎨ S∈I S 3,5 ∑ 1 ⎪ ⎪ ϖ Si ⟨ξ1 ξ2 ⟩ = , ⎪ ⎪ 12 ⎪ ⎪ S∈I S1,3−6 ⎪ ⎪ ⎪ ∑ ⎪ 1 ⎪ ⎪ ϖ Si ⟨ξ12 ξ22 ⟩ = ⎪ ⎪ 180 ⎩ S∈I S1,3−6

ϖ Si ⟨1⟩ = 1

1 6

ϖ Si ⟨1⟩ = 1

S∈I S1−6

∑ S∈I S1,4,6

ϖ Si ⟨ξ1 ξ2 ξ3 ⟩ =

1 60

(B.4)

668

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

(5) 3D quadratic case ⎧ ∑ ⎪ ϖ Sb ⟨δˆ4 ⟩ = 1, ⎪ ⎪ ⎪ ⎪ ⎪ S∈I S2,3,5,7,8,10 ⎪ ⎪ ⎪ 1 ⎨ ∑ , ϖ Sb ⟨ξ1 ξ2 ξ3 δˆ4 ⟩ = 60 S∈I ⎪ S ⎪ 3,8,10 ⎪ ⎪ ∑ 1 ⎪ ⎪ ⎪ ϖ Si ⟨ξ1 ξ2 ⟩ = ⎪ ⎪ 20 ⎩ S∈I

∑

ϖ Sb ⟨ξ1 ξ2 δˆ4 ⟩ =

S∈I S3,5,7,8,10

∑

1 12

ϖ Si ⟨1⟩ = 1

(B.5)

S∈I S1−11

S1,3−11

(6) 3D cubic case ∑ ⎧ ϖ Sb ⟨δˆ4 ⟩ = 1, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ S∈I S2,3,5,7,8,10 ⎪ ⎪ ∑ ⎪ 1 ⎪ ⎪ ϖ Sb ⟨ξ1 ξ2 ξ3 δˆ4 ⟩ = , ⎪ ⎪ 60 ⎪ S∈I ⎪ S3,8,10 ⎪ ⎪ ⎪ ∑ ⎪ 1 ⎨ ϖ Sb ⟨ξ12 ξ22 ξ3 δˆ4 ⟩ = , 210 S∈I ⎪ S3,8,10 ⎪ ⎪ ⎪ ∑ 1 ⎪ ⎪ ⎪ , ϖ Si ⟨ξ1 ξ2 ⟩ = ⎪ ⎪ 20 ⎪ ⎪ S∈I S1,3−11 ⎪ ⎪ ⎪ ∑ ⎪ 1 ⎪ ⎪ ϖ Si ⟨ξ12 ξ22 ⟩ = , ⎪ ⎪ ⎩ 210 S∈I S1,3−11

∑

ϖ Sb ⟨ξ1 ξ2 δˆ4 ⟩ =

S∈I S3,5,7,8,10

∑

ϖ Sb ⟨ξ12 ξ22 δˆ4 ⟩ =

S∈I S3,8,10

∑

1 12

1 180

ϖ Si ⟨1⟩ = 1

(B.6)

S∈I S1−11

∑

ϖ Si ⟨ξ1 ξ2 ξ3 ⟩ =

S∈I S1,3,4,6,8−11

∑ S∈I S1,4,6,9−11

ϖ Si ⟨ξ1 ξ2 ξ3 ξ4 ⟩ =

1 120 1 840

Moreover, the constraints corresponding to Eq. (75) for the sample point optimization problem are provided for meshfree formulations using quadratic and cubic basis functions with details as follows: (1) 1D quadratic case ⎧ m A ≥ 0, A = 1, 2, 3 ⎪ ⎪ ⎪ ⎪ ⎨m 1 , m 2 ≤ 1 m 1 + 2m 2 + 2m 3 ≥ 3 ⎪ ⎪ ⎪m 1 + 2m 3 ≥ 1 ⎪ ⎩ m2 ≥ 1 (2) 1D cubic case ⎧ ⎪ ⎪m A ≥ 0, A = 1, 2, 3 ⎪ ⎪ ⎨m 1 , m 2 ≤ 1 m 1 + 2m 2 + 2m 3 ≥ 4 ⎪ ⎪ m ⎪ 1 + 2m 3 ≥ 2 ⎪ ⎩ m2 ≥ 1 (3) 2D quadratic case ⎧ m A ≥ 0, A = 1, 2, . . . , 6 ⎪ ⎪ ⎪ ⎪ m ⎪ 1, m2, m3 ≤ 1 ⎪ ⎨ m 1 + 2m 2 + 2m 3 + 2m 4 + 3m 5 + 3m 6 ≥ 4 m 2 + m 3 + 2m 5 ≥ 2 ⎪ ⎪ ⎪ ⎪ m 1 + 2m 3 + 2m 4 + 3m 5 + 3m 6 ≥ 1 ⎪ ⎪ ⎩ m 3 + 2m 5 ≥ 1

(B.7)

(B.8)

(B.9)

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

(4) 2D cubic case ⎧ ⎪ ⎪m A ≥ 0, A = 1, 2, . . . , 6 ⎪ ⎪ ⎪ ⎪m 1 , m 2 , m 3 ≤ 1 ⎪ ⎪ ⎨m 1 + 2m 2 + 2m 3 + 2m 4 + 3m 5 + 3m 6 ≥ 7 m 2 + m 3 + 2m 5 ≥ 3 ⎪ ⎪ m 1 + 2m 3 + 2m 4 + 3m 5 + 3m 6 ≥ 4 ⎪ ⎪ ⎪ ⎪ m 1 + 2m 4 + 3m 6 ≥ 1 ⎪ ⎪ ⎩ m 3 + 2m 5 ≥ 2 (5) 3D quadratic case ⎧ m A ≥ 0, A = 1, 2, . . . , 11 ⎪ ⎪ ⎪ ⎪ m ⎪ 1, m2, m3, m5 ≤ 1 ⎪ ⎪ ⎪ ⎨(m 1 + 2m 2 + 2m 3 + 2m 4 + 2m 5 + 2m 6 + 3m 7 + 3m 8 + 3m 9 + 4m 10 + 4m 11 ) ≥ 5 (m 1 + 2m 3 + 2m 4 + 2m 5 + 2m 6 + 3m 7 + 3m 8 + 3m 9 + 4m 10 + 4m 11 ) ≥ 3 ⎪ ⎪ m 2 + m 3 + m 5 + 2m 7 + 2m 8 + 3m 10 ≥ 3 ⎪ ⎪ ⎪ ⎪ m 3 + m 5 + 2m 7 + 2m 8 + 3m 10 ≥ 2 ⎪ ⎪ ⎩ m 3 + 2m 8 + 3m 10 ≥ 1 (6) 3D cubic case ⎧ m A ≥ 0, A = 1, 2, . . . , 11 ⎪ ⎪ ⎪ ⎪ m1, m2, m3, m5 ≤ 1 ⎪ ⎪ ⎪ ⎪ (m 1 + 2m 2 + 2m 3 + 2m 4 + 2m 5 + 2m 6 + 3m 7 + 3m 8 + 3m 9 + 4m 10 + 4m 11 ) ≥ 10 ⎪ ⎪ ⎪ ⎪ ⎨(m 1 + 2m 3 + 2m 4 + 2m 5 + 2m 6 + 3m 7 + 3m 8 + 3m 9 + 4m 10 + 4m 11 ) ≥ 7 m 1 + m 3 + 2m 4 + 2m 6 + 3m 8 + 3m 9 + 4m 10 + 4m 11 ≥ 5 ⎪ ⎪ ⎪m 2 + m 3 + m 5 + 2m 7 + 2m 8 + 3m 10 ≥ 5 ⎪ ⎪ ⎪ ⎪m 3 + m 5 + 2m 7 + 2m 8 + 3m 10 ≥ 4 ⎪ ⎪ ⎪ ⎪m 3 + 2m 8 + 3m 10 ≥ 3 ⎪ ⎩ m 1 + 2m 4 + 2m 6 + 3m 9 + 4m 11 ≥ 1

669

(B.10)

(B.11)

(B.12)

References [1] T. Belytschko, Y. Kronggauz, D. Organ, M. Fleming, P. Krysl, Meshless methods: an overview and recent developments, Comput. Methods Appl. Mech. Engrg. 139 (1996) 3–47. [2] J.S. Chen, M. Hillman, S.W. Chi, Meshfree methods: progress made after 20 years, J. Eng. Mech.-ASCE 143 (2017) 04017001. [3] B. Nayroles, G. Touzot, P. P. Villon, Generalizing the finite element method: diffuse approximation and diffuse elements, Comput. Mech. 10 (1992) 307–318. [4] P. Lancaster, K. Salkauskas, Surfaces generated by moving least squares methods, Math. Comp. 37 (1981) 141–158. [5] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Internat. J. Numer. Methods Engrg. 37 (1994) 229–256. [6] T. Belytschko, Y.Y. Lu, L. Gu, Crack propagation by element-free Galerkin methods, Eng. Fract. Mech. 51 (1995) 295–315. [7] W.K. Liu, S. Jun, Y.F. Zhang, Reproducing kernel particle methods, Internat. J. Numer. Methods Fluids 20 (1995) 1081–1106. [8] J.S. Chen, C. Pan, C.T. Wu, W.K. Liu, Reproducing kernel particle methods for large deformation analysis of nonlinear structures, Comput. Methods Appl. Mech. Engrg. 139 (1996) 195–227. [9] C.A. Duarte, J.T. Oden, An hp adaptive method using clouds, Comput. Methods Appl. Mech. Engrg. 139 (1996) 237–262. [10] I. Babuška, J.M. Melenk, The partition of unity method, Internat. J. Numer. Methods Engrg. 40 (1997) 727–758. [11] N. Sukumar, B. Moran, T. Belytschko, The natural element method in solid mechanics, Internat. J. Numer. Methods Engrg. 43 (1998) 839–887. [12] S.N. Atluri, T. Zhu, A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics, Comput. Mech. 22 (1998) 117–127. [13] S. De, K.J. Bathe, The method of finite spheres, Comput. Mech. 25 (2000) 329–345. [14] J.G. Wang, G.R. Liu, A point interpolation meshless method based on radial basis functions, Internat. J. Numer. Methods Engrg. 54 (2002) 1623–1648. [15] T. Strouboulis, I. Babuška, K. Copps, The design and analysis of the generalized finite element method, Comput. Methods Appl. Mech. Engrg. 181 (2000) 43–69. [16] J.S. Chen, W. Han, Y. You, X. Meng, A reproducing kernel method with nodal interpolation property, Internat. J. Numer. Methods Engrg. 56 (2003) 935–960. [17] L. Gu, Moving kriging interpolation and element-free Galerkin method, Internat. J. Numer. Methods Engrg. 56 (2003) 1–11.

670

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

[18] W.K. Liu, W. Han, H. Lu, S. Li, J. Cao, Reproducing kernel element method. Part I: theoretical formulation, Comput. Methods Appl. Mech. Engrg. 193 (2004) 933–951. [19] N. Sukumar, Construction of polygonal interpolants: a maximum entropy approach, Internat. J. Numer. Methods Engrg. 61 (2004) 2159–2181. [20] M. Arroyo, M. Ortiz, Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods, Internat. J. Numer. Methods Engrg. 65 (2006) 2167–2202. [21] G.R. Liu, K.Y. Dai, T.T. Nguyen, A smoothed finite element method for mechanics problems, Comput. Mech. 39 (2007) 859–877. [22] B. Li, F. Habbal, M. Ortiz, Optimal transportation meshfree approximation schemes for fluid and plastic flows, Internat. J. Numer. Methods Engrg. 83 (2010) 1541–1579. [23] P.C. Guan, S.W. Chi, J.S. Chen, T.R. Slawson, M.J. Roth, Semi-Lagrangian reproducing kernel particle method for fragment-impact problems, Int. J. Impact Eng. 38 (2011) 1033–1047. [24] C.T. Wu, C.K. Park, J.S. Chen, A generalized approximation for the meshfree analysis of solids, Internat. J. Numer. Methods Engrg. 85 (2011) 693–722. [25] M.A. Bessa, J.T. Foster, T. Belytschko, W.K. Liu, A meshfree unification: reproducing kernel peridynamics, Comput. Mech. 53 (2014) 1251–1264. [26] D. Wang, P. Chen, Quasi-convex reproducing kernel meshfree method, Comput. Mech. 54 (2014) 689–709. [27] C.T. Wu, M. Koishi, W. Hu, A displacement smoothing induced strain gradient stabilization for the meshfree Galerkin nodal integration method, Comput. Mech. 56 (2015) 19–37. [28] P. Metsis, N. Lantzounis, M. Papadrakakis, A new hierarchical partition of unity formulation of EFG meshless methods, Comput. Methods Appl. Mech. Engrg. 283 (2015) 782–805. [29] E. Yreux, J.S. Chen, A quasi-linear reproducing kernel particle method, Internat. J. Numer. Methods Engrg. 109 (2017) 1045–1064. [30] Y. Wu, C.T. Wu, Simulation of impact penetration and perforation of metal targets using the smoothed particle Galerkin method, J. Eng. Mech.-ASCE 144 (2018) 04018057. [31] S.N. Atluri, S.P. Shen, The Meshless Local Petrov–Galerkin (MLPG) Method, Tech Science, 2002. [32] I. Babuška, U. Banerjee, J.E. Osborn, Survey of meshless and generalized finite element methods: a unified approach, Acta Numer. 12 (2003) 1–125. [33] S. Li, W.K. Liu, Meshfree Particle Methods, Springer-Verlag, 2004. [34] X. Zhang, Y. Liu, Meshless Methods, Tsinghua University Press & Springer-Verlag, 2004. [35] V.P. Nguyen, T. Rabczuk, S. Bordas, M. Duflot, Meshless methods: a review and computer implementation aspects, Math. Comput. Simulation 79 (2008) 763–813. [36] G.R. Liu, Meshfree Methods: Moving beyond The Finite Element Method (2nd edition), CRC Press, 2009. [37] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, Mineola, New York, 2000. [38] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The Finite Element Method: Its Basis and Fundamentals (7th Edition), Elsevier, Singapore, 2015. [39] J. Dolbow, T. Belytschko, Numerical integration of the Galerkin weak form in meshfree methods, Comput. Mech. 23 (1999) 219–230. [40] J.S. Chen, C.T. Wu, S. Yoon, Y. You, A stabilized conforming nodal integration for Galerkin meshfree methods, Internat. J. Numer. Methods Engrg. 50 (2001) 435–466. [41] J.S. Chen, S. Yoon, C.T. Wu, Nonlinear version of stabilized conforming nodal integration for Galerkin meshfree methods, Internat. J. Numer. Methods Engrg. 53 (2002) 2587–2615. [42] I. Babuška, U. Banerjee, J.E. Osborn, Q.L. Li, Quadrature for meshless methods, Internat. J. Numer. Methods Engrg. 76 (2008) 1434–1470. [43] A. Karatarakis, P. Metsis, M. Papadrakakis, GPU-Acceleration of stiffness matrix calculation and efficient initialization of EFG meshless methods, Comput. Methods Appl. Mech. Engrg. 258 (2013) 63–80. [44] C. Peco, D. Milian, A. Rosolen, M. Arroyo, Efficient implementation of Galerkin meshfree methods for large-scale problems with an emphasis on maximum entropy approximants, Comput. Struct. 150 (2015) 52–62. [45] S. Beissl, T. Belytschko, Nodal integration of the element-free Galerkin method, Comput. Methods Appl. Mech. Engrg. 139 (1996) 49–64. [46] C.T. Dyka, P.W. Randles, R.P. Ingel, Stress points for tension instability in SPH, Internat. J. Numer. Methods Engrg. 40 (1997) 2325–2341. [47] T. Rabczuk, T. Belytschko, S.P. Xiao, Stable particle methods based on Lagrangian kernels, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1035–1063. [48] C.T. Wu, S.W. Chi, M. Koishi, Y. Wu, Strain gradient stabilization with dual stress points for the meshfree nodal integration method in inelastic analyses, Internat. J. Numer. Methods Engrg. 107 (2016) 3–30. [49] C.T. Wu, Y. Wu, Z. Liu, D. Wang, A stable and convergent Lagrangian particle method with multiple nodal stress points for large strain and material failure analyses in manufacturing processes, Finite Elem. Anal. Des. 146 (2018) 96–106. [50] A. Carpinteri, G. Ferro, G. Ventura, The partition of unity quadrature in meshless methods, Internat. J. Numer. Methods Engrg. 54 (2002) 987–1006. [51] M. Griebel, M.A. Schweitzer, A particle-partition of unity method. Part II: efficient cover construction and reliable integration, SIAM J. Sci. Comput. 23 (2002) 1655–1682. [52] K.C. Kwon, S.H. Park, S.K. Youn, The support integration scheme in the least-squares mesh-free method, Finite Elem. Anal. Des. 43 (2006) 127–144.

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

671

[53] Y. Liu, T. Belytschko, A new support integration scheme for the weakform in mesh-free methods, Internat. J. Numer. Methods Engrg. 82 (2010) 699–715. [54] J.S. Chen, C.T. Wu, T. Belytschko, Regularization of material instabilities by meshfree approximations with intrinsic length scales, Internat. J. Numer. Methods Engrg. 47 (2000) 1301–1322. [55] D. Wang, J.S. Chen, Locking-free stabilized conforming nodal integration for meshfree Mindlin-Reissner plate formulation, Comput. Methods Appl. Mech. Eng. 193 (2004) 1065–1083. [56] J.S. Chen, D. Wang, A constrained reproducing kernel particle formulation for shear deformable shell in Cartesian coordinates, Internat. J. Numer. Methods Engrg. 68 (2006) 151–172. [57] D. Wang, S.B. Dong, J.S. Chen, Extended meshfree analysis of transverse and inplane loading of a laminated anisotropic plate of general planform geometry, Int. J. Solids Struct. 43 (2006) 144–171. [58] D. Wang, Y. Wu, An efficient Galerkin meshfree analysis of shear deformable cylindrical panels, Interact. Multiscale Mech. 1 (2008) 339–355. [59] D. Wang, Y. Sun, A Galerkin meshfree method with stabilized conforming nodal integration for geometrically nonlinear analysis of shear deformable plates, 8 (2011) 685–703. [60] D. Wang, Z. Li, L. Li, Y. Wu, Three dimensional efficient meshfree simulation of large deformation failure evolution in soil medium, Sci. China Technol. Sci. 54 (2011) 573–580. [61] Y. Wu, D. Wang, C.T. Wu, Three dimensional fragmentation simulation of concrete structures with a nodally regularized meshfree method, Theor. Appl. Fract. Mech. 72 (2014) 89–99. [62] Y. Wu, D. Wang, C.T. Wu, H. Zhang, A direct displacement smoothing meshfree particle formulation for impact failure modeling, Int. J. Impact Eng. 87 (2016) 169–185. [63] M. Hillman, J.S. Chen, Nodally integrated implicit gradient reproducing kernel particle method for convection dominated problems, Comput. Methods Appl. Mech. Engrg. 299 (2016) 381–400. [64] D. Wang, J.S. Chen, A Hermite reproducing kernel approximation for thin-plate analysis with sub-domain stabilized conforming integration, Internat. J. Numer. Methods Engrg. 74 (2008) 368–390. [65] D. Wang, Z. Lin, Free vibration analysis of thin plates using Hermite reproducing kernel Galerkin meshfree method with sub-domain stabilized conforming integration, Comput. Mech. 46 (2010) 703–719. [66] D. Wang, Z. Lin, Dispersion and transient analyses of hermite reproducing kernel Galerkin meshfree method with sub-domain stabilized conforming integration for thin beam and plate structures, Comput. Mech. 48 (2011) 47–63. [67] D. Wang, H. Peng, A hermite reproducing kernel galerkin meshfree approach for buckling analysis of thin plates, Comput. Mech. 51 (2013) 1013–1029. [68] D. Wang, C. Song, H. Peng, A circumferentially enhanced hermite reproducing kernel meshfree method for buckling analysis of Kirchhoff–love cylindrical shells, Int. J. Struct. Stab. Dyn. 15 (2015) 1450090. [69] J.S. Chen, Y. Wu, Stability in Lagrangian and semi-Lagrangian reproducing kernel discretizations using nodal integration in nonlinear solid mechanics, in: Advances in Meshfree Techniques, Springer, 2007, pp. 55–76. [70] Q. Duan, X. Li, H. Zhang, T. Belytschko, Second-order accurate derivatives and integration schemes for meshfree methods, Internat. J. Numer. Methods Engrg. 92 (2012) 399–424. [71] Q. Duan, X. Gao, B. Wang, X. Li, H. Zhang, A four-point integration scheme with quadratic exactness for three-dimensional element-free Galerkin method based on variationally consistency formulation, Comput. Methods Appl. Mech. Engrg. 280 (2014) 84–116. [72] J.S. Chen, M. Hillman, M. Rüter, An arbitrary order variationally consistent integration for Galerkin meshfree methods, Internat. J. Numer. Methods Engrg. 95 (2013) 387–418. [73] M. Hillman, J.S. Chen, S.W. Chi, Stabilized and variationally consistent nodal integration for meshfree modeling of impact problems, Comput. Part. Mech. 1 (2014) 245–256. [74] M. Ruter, M. Hillman, J.S. Chen, Corrected stabilized non-conforming nodal integration in meshfree methods, Lect. Notes Comput. Sci. Eng. 89 (2013) 75–93. [75] D. Wang, J. Wu, An efficient nesting sub-domain gradient smoothing integration algorithm with quadratic exactness for Galerkin meshfree methods, Comput. Methods Appl. Mech. Engrg. 298 (2016) 485–519. [76] M. Hillman, J.S. Chen, An. accelerated, convergent, An accelerated convergent and stable nodal integration in Galerkin meshfree methods for linear and nonlinear mechanics, Internat. J. Numer. Methods Engrg. 107 (2016) 603–630. [77] S. Li, W.K. Liu, Synchronized reproducing kernel interpolant via multiple wavelet expansion, Comput. Mech. 21 (1998) 28–47. [78] S. Li, W.K. Liu, Reproducing kernel hierarchical partition of unity part I: formulation and theory, Internat. J. Numer. Methods Engrg. 45 (1999) 251–288. [79] S. Li, W.K. Liu, Reproducing kernel hierarchical partition of unity part II: applications, Internat. J. Numer. Methods Engrg. 45 (1999) 289–317. [80] S.W. Chi, J.S. Chen, H.Y. Hu, J.P. Yang, A gradient reproducing kernel collocation method for boundary value problems, Internat. J. Numer. Methods Engrg. 93 (2013) 1381–1402. [81] E.L. Wachspress, A Rational Finite Element Basis: Mathematics in Science and Engineering, Elsevier, 1975. [82] M.A. Eisenberg, L.E. Malvern, On finite element integration in natural co-ordinates, Internat. J. Numer. Methods Engrg. 7 (1973) 574–575. [83] D.A. Dunavant, High degree efficient symmetrical Gaussian quadrature-rules for the triangle, Internat. J. Numer. Methods Engrg. 21 (1985) 1129–1148. [84] M. Gellert, R. Harbord, Moderate degree cubature formulas for 3-D tetrahedral finite-element approximations, Commun. Appl. Numer. Methods 7 (1991) 487–495.

672

D. Wang and J. Wu / Computer Methods in Applied Mechanics and Engineering 349 (2019) 628–672

[85] J.H. Conway, S. Torquato, Packing, tiling, and covering with tetrahedra, in: Proceedings of National Academy of Sciences of USA, Vol. 103, 2006, pp. 10612–10617. [86] A. Schrijver, Theory of Linear and Integer Programming, John Wiley & Sons, 1998. [87] S. Fernández-Méndez, A. Huerta, Imposing essential boundary conditions in mesh-free methods, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1257–1275. [88] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity (3rd Edn), McGraw-Hill, 1970. [89] P. Wei, W. Xiao, Area calculation of three dimensional polygon, Chin. Math. Bull. 2 (1984) 18–21.