Generalized reproducing kernel peridynamics: unification of local and non-local meshfree methods, non-local derivative operations, and an arbitrary-order state-based peridynamic formulation

State-based peridynamics is a non-local reformulation of solid mechanics that replaces the force density of the divergence of stress with an integral of the action of force states on bonds local to a given position, which precludes differentiation with the aim to model strong discontinuities effortlessly. A popular implementation is a meshfree formulation where the integral is discretized by quadrature points, which results in a series of unknowns at the points under the strong-form collocation framework. In this work, the meshfree discretization of state-based peridynamics under the correspondence principle is examined and compared to traditional meshfree methods based on the classical local formulation of solid mechanics. It is first shown that the way in which the peridynamic formulation approximates differentiation can be unified with the implicit gradient approximation, and this is termed the reproducing kernel peridynamic approximation. This allows the construction of non-local deformation gradients with arbitrary-order accuracy, as well as non-local approximations of higher-order derivatives. A high-order accurate non-local divergence of stress is then proposed to replace the force density in the original state-based peridynamics, in order to obtain global arbitrary-order accuracy in the numerical solution. These two operators used in conjunction with one another is termed the reproducing kernel peridynamic method. The strong-form collocation version of the method is tested against benchmark solutions to examine and verify the highorder accuracy and convergence properties of the method. The method is shown to exhibit superconvergent behavior in the nodal collocation setting.


Introduction
Meshfree methods for continuum mechanics can generally be cast under two formulations: Galerkin meshfree methods and collocation meshfree methods [1]. Among these, the governing equations solved under these frameworks can also be different. On one branch, there are meshfree discretizations of the classical local model, namely the differential equation for linear momentum. Using the Galerkin formulation with a meshfree discretization of this equation results in methods such as the diffuse element method [2], the element free Galerkin method [3], the reproducing kernel particle method [4,5], and the many methods that ensued thereafter [6][7][8][9][10][11][12][13]. On the other hand, due to the global smoothness that can easily be attained with meshfree approximations, collocation of the strong form is straightforward, and several methods using this technique have been proposed starting with radial basis functions [14,15], and later with the reproducing kernel and moving least squares approximations [16][17][18], among others.
Recently, a non-local reformulation of continuum mechanics has been proposed [19,20] called peridynamics, in order to circumvent difficulties in treating discontinuities in the local models, first under a bond-based framework [19], and later under the so-called state-based framework [20]. This theory replaces the differentiation in classical continuum mechanics, i.e., the force density of the divergence of stress, with an integral of the actions of bonds to a given position, and admits discontinuities by precluding derivatives in the governing equations. In addition, the action is considered over a finite distance, embedding a length scale in the governing equations, resulting in a formulation which is non-local in nature. While the peridynamic equations can be solved under the Galerkin formulation [21,22], this approach requires evaluation of a double (six-dimensional) integral [23], which results in considerable computational expense. Therefore, for practical applications, the strong form version is often employed. Due to the simplicity of its implementation and relatively low-computational cost, the peridynamic strong form equations are generally solved by a node-based meshfree approach [24], which is based on nodal collocation and nodal discretization of the integral terms in the peridynamic equations. In this meshfree approach, the unknowns are therefore associated with the nodal points, which are both collocation and quadrature points, similar to smoothed particle hydrodynamics [25,26]. In this paper, discussions are focused on the collocation meshfree implementation of the state-based version under the correspondence principle [20], which constructs a non-local deformation gradient to facilitate the use of classical constitutive models, and is most closely related to meshfree discretizations of local models as first shown in [27] for uniform discretizations and infinite domains. In particular, the main focus of this study is the accuracy and convergence properties of the state-based peridynamic method, and enhancement thereof.
Theoretical analysis of the accuracy and convergence of local meshfree methods is well established, for instance see [28][29][30][31] for Galerkin analysis, and [18,[32][33][34][35] for collocation analysis. Sufficiently smooth problems solved with monomial basis vectors exhibit algebraic convergence in both strong formulations [35], and weak formulations [29], while global approximations such as radial basis show exponential convergence [33]. In particular, Galerkin meshfree methods with nth monomial completeness exhibit a rate of n ? 1 in the displacement solution. On the other hand, meshfree collocation approaches exhibit a rate of n ? 1 for a least-squares formulation using more collocation points than source points (approximation functions) [18], while using an equal number of source and collocation points exhibits an odd-even phenomenon where the rate of n is obtained for even orders, and n -1 for odd orders [17,[36][37][38], which has also been observed in isogeometric collocation [39,40]. A recently developed recursive gradient formulation [37] has been developed that exhibits superconvergence, that is, rates of n and n ? 1 for even and odd orders of approximations, respectively. Notably this allows linear basis to converge in collocation analysis, in contrast to direct gradients [18]. Finally, it should be noted that the accuracy of numerical integration in the weak-form based versions can heavily influence theoretical rates [41], although several approaches are available to rectify this situation (the interested reader is referred to [1] and references therein for details).
For peridynamics, the concept of convergence can be understood in several ways; see [42,43] for details: (1) Nconvergence, when the non-locality of the continuous peridynamic problem is kept fixed while the discretization is refined (convergence to the non-local continuum solution); (2) d-convergence, when the nonlocality is reduced for a fixed discretization (convergence to the discrete local solution); or (3) N-d convergence, when both discretization and nonlocality approach the vanishing limit simultaneously (convergence to the continuum local solution). In this work, the third type of convergence is studied. It has been noted that this type of convergence is linked to the concept of asymptotic compatibility, with discretizations being asymptotically compatible if they converge to the correct local model and solution associated with the nonlocal model [43,44].
The accuracy and convergence properties of state-based peridynamics has been studied in several works. Material models play an important role in convergence study of peridynamics since the results using different models converge to different solutions [45]. In [45], three different material models were tested to attempt to reproduce the solution of static linear elasticity. Among these, only the deformation gradient-based model [20] gave a convergent solution to all problems tested, with a first-order convergence rate in the L 2 error norm.
The effect of numerical integration in the non-local integrals has also been a focus in studies, as it may also have a strong effect on convergence in peridynamics [46][47][48][49][50]. In the first meshfree implementation of peridynamics [24], the peridynamic equation of motion was discretized by nodal integration with the full physical nodal volume as the integration weight, resulting in the so-called full volume (FV) integration. The FV integration shows erratic convergence behavior, both converging and diverging with refinement [46,51]. Several studies suggest this issue is due to rough approximation of the integration weights near the edges of the integration domains in the peridynamic equations [46][47][48][49]. So-called partial volume (PV) integration schemes have been proposed, including approximate PV [47,48] and analytical PV [49], to more accurately compute the partial volumes intersecting with neighborhoods that serve as integration weights for particles. Partial volume schemes have been shown to improve the accuracy and yield more consistent convergence rates compared with the FV integration [46]. Influence functions Computational Particle Mechanics that smoothly decay to zero at the boundary of the horizon have also been investigated under the FV integration framework, and it has been shown that this enhances the accuracy and can also give more consistent convergence behavior [46]. The idea behind smoothly decaying influence functions is to reduce the influence of the particles near the neighborhood boundary, mitigating the error due to numerical integration. These techniques exhibit firstorder convergence in displacements [46]. The limitation of first-order convergence was attributed in [50] to the piecewise constant nature of the approximation employed, where the reproducing kernel approximation was introduced in the peridynamic displacement field to increase the convergence rate. However, high-order Gauss integration was employed to achieve the integration accuracy necessary to avoid the aforementioned oscillatory convergence behavior, resulting in an increased computational cost. A high-order non-local deformation gradient was proposed in [52] for stability reasons, but was limited to uniform discretizations, and meanwhile the convergence properties of this method were not tested.
Taking another point of view, the non-local deformation gradient and force density in state-based peridynamics under the correspondence principle can be viewed as mathematical operators with certain approximation properties. In [27,53] it was shown that for uniform discretizations away from the boundary, the accuracy in the non-local deformation gradient is second-order. This is confirmed by other studies, where it has been further shown that near the boundary the accuracy will be first-order [52,54].
The precise relationship between the meshfree peridynamic method and the classical meshfree methods discussed has not been made clear. So far, one effort [27] has attempted to examine the relationship between the peridynamic approximation of derivatives via the non-local deformation gradient and the traditional meshfree approximation of derivatives. There it was shown that the statebased peridynamic formulation based on correspondence is equivalent to employing a second-order accurate implicit gradient reproducing kernel approximation [55], for both the deformation gradient operation on displacement, and force density operation on the stress, but this equivalence was established only for uniform discretizations, away from a boundary. However, this relationship also implies that these operations are both second-order accurate, at least in uniform discretizations, and away from the influence of a boundary.
In this paper, the precise relationship between meshfree methods for local models and non-local peridynamic meshfree discretizations under the correspondence principle is introduced, for general non-uniform discretizations, and finite domains. A generalized approximation which unifies these approaches is introduced termed the reproducing kernel peridynamic approximation, under both continuous (integral form) and discrete frameworks. It is shown that this approximation can yield four distinct cases: implicit gradients, the traditional non-local deformation gradient, as well as an arbitrary-order accurate non-local deformation gradient, and arbitrary-order accurate non-local higher-order derivatives. A formulation is then proposed called the reproducing kernel peridynamic (RKPD) method, consisting of the high-order accurate non-local deformation gradients, in conjunction with a high-order accurate force density, which results in an arbitrary-order accurate state-based peridynamic method. The formulation is tested under the node-based collocation framework, although a weak formulation is also possible. In contrast to the original formulation, the method is shown to exhibit convergent solutions with and without ghost boundary nodes, under both uniform and nonuniform discretizations, with superconvergent solutions for odd orders of accuracy.
The remainder of this paper is organized as follows. In Sect. 2, the governing equations for classical local methods and state-based peridynamics are briefly reviewed. The integral forms for the reproducing kernel and implicit gradient approximation are given in Sect. 3, and compared with the integral forms of the state-based peridynamic equations. In addition, the equivalence of implicit gradients and the peridynamic differential operator [56] is established. These formulations are then compared and contrasted, and the orders of accuracy are assessed. In Sect. 4, the continuous reproducing kernel peridynamic approximation is given, which unifies the two formulations, and provides arbitrary-order accurate non-local deformation gradients, and arbitrary-order accurate higher-order nonlocal derivative approximations. The discrete implicit gradient and peridynamic approximations are then discussed and compared in Sect. 5, with the order of accuracy in the discrete case assessed. Section 6 introduces the discrete reproducing kernel peridynamic approximation. The collocation implementation of the proposed formulation, the reproducing kernel peridynamic method, is then summarized in Sect. 7, and numerical examples are given in Sect. 8. Conclusions, and discussions on implications and possible future work are given in Sect. 9.

Governing equations
In this section, the governing equations for classical continuum mechanics and state-based peridynamics under the correspondence principle are briefly reviewed.

Classical continuum mechanics
The equation of motion for finite-strain continuum mechanics problems stated in the reference configuration X at material position X at time t is € uðX; tÞqðXÞ ¼ r Á r T ðX; tÞ þ bðX; tÞ ð 1Þ where € u Dv=Dt is the material time derivative of the velocity v, q is the material density in the undeformed configuration, r is the 1st Piola-Kirchhoff (PK) stress tensor (r T is the nominal stress), r denotes the del operator with respect to the undeformed configuration, and b is the body force in the undeformed configuration. In this work a Lagrangian description is adopted, as state-based peridynamics under correspondence relates the 1st PK stress to the force density in the governing equations [20].
Given a strain energy density function WðFÞ, kinetic variables such as the first 1st PK stress r can be obtained as r ¼ oWðFÞ=oF. In state-based peridynamics, an analogous relationship exists between the kinematic and kinetic entities, as described in the next section.

State-based peridynamics
In order to deal with discontinuities, peridynamics [19,20] has been introduced which precludes the differentiation involved in the governing equations for classical continuum mechanics (1). In state-based peridynamics, the force density of the divergence of nominal stress in (1) is replaced by an integral of force states T [20]: where H X is the so-called neighborhood of the particle X, which is often defined by a sphere encompassing the point X with radius d called the horizon, as shown in Fig. 1. The notation employed for the mathematical entity of states is that angle brackets denote the operation on that variable, while square brackets denote the dependence on the variable. A fundamental kinetic entity is the force state ThÁi, rather than for instance, the 1st PK stress in the Lagrangian formulation of classical solid mechanics. The quantity X 0 À X is said to be a ''bond'' of the points X 0 and X. Thus, it can be seen when comparing (1) to (2) that the force density of the divergence of stress is replaced by an integral of the action of states T on bonds local to X.
A fundamental kinematic entity in peridynamics is the deformation state YhÁi which maps (possibly nonlinearly and discontinuously) a bond in the undeformed configuration X 0 À X, to a bond in the current configuration x 0 À x, as shown in Fig. 1: Analogous to the dependence of stress measures on strain measures in classical mechanics, the force state T depends on the deformation state. On the other hand, in order to facilitate the use of constitutive models in the local theory, a non-local deformation gradient F can be obtained through a principle called reduction [20], which is briefly reviewed in Sect. 3.5. The 1st PK stress can be obtained via F , and can then be related to the force state T by means of energy principles.

Continuous reproducing kernel and peridynamic approximations
As will be demonstrated, the implicit gradient counterpart [36,55,57,58] to the reproducing kernel (RK) approximation [4,5] is most closely related to the non-local deformation gradient employed in peridynamics. In this section, the continuous (integral) form of the reproducing Computational Particle Mechanics kernel and peridynamic approximations are analyzed and compared. A brief review of both is given. The continuous versions must be discretized in practice; the discrete versions will be analyzed and compared in Sect. 5. Remarkably, if the quadrature and discretizations are consistent with one another, the main results of these analyses are the same, with minor exceptions.

Continuous reproducing kernel approximation
The continuous RK approximation of a function uðxÞ on a domain X & R d is constructed by the product of a kernel function U a with compact support, and a correction function composed of a linear combination of basis functions in the following form [4,5]: where HðxÞ is a column vector of complete nth order monomials (although other bases could be employed), and bðxÞ is a column vector of associated coefficients to be determined. The dependence of these vectors in the RK approximation on the free parameter n is to be understood herein for notational simplicity.
It should be noted that, in much of the literature, the shifted basis term x À x 0 is employed, while here the basis using x 0 À x is employed in order to unify the reproducing kernel approximation and peridynamic derivative approximation later in the text. The choice is arbitrary, and only results in sign differences in gradient reproducing conditions.
The kernel function U a has compact support with measure a, and the smoothness of the approximation is inherited from the kernel. For example, using C 2 kernels yields C 2 continuity of the approximation. In this work, the cubic B-spline kernel is employed for kernel functions and influence functions, which play the analogous role in peridynamics as discussed in Sect. 3.5. The one-dimensional cubic B-spline kernel shown in Fig. 2 is constructed as: In multiple dimensions, one may construct a kernel by tensor product yielding a box or cuboid support: or by defining U a ðx 0 À xÞ ¼ U a ðzÞ with z ¼ jx 0 À xj=a, yielding a spherical support. When monomials are employed for HðxÞ, the coefficients bðxÞ are determined by enforcing nth order accuracy of the approximation in (4). This can be achieved by directly enforcing the so-called reproducing conditions (discussed later), or by using a Taylor expansion. In this work, the latter approach is employed in order to fully illustrate the meaning and construction of an implicit gradient. The Taylor series expansion of uðx 0 Þ around x truncated to order n is: where a ¼ ða 1 ; . . .; a d Þ is a multi-index in R d of non-negative integers equipped with the notation a In matrix form (7) can be expressed as: where DðxÞ is a row vector of fo b uðxÞg n jbj¼0 and J is a diagonal matrix with entries f1=b!g n b j j¼0 . Substituting (8) into (4) yields The nth order accuracy of the approximation requires that R ½n fu x ð Þg ¼ u x ð Þ in the above. Examining (9), this can be phrased as the following vanishing moment conditions: Z where the fact that f1=b!gj b j j¼0 ¼ 1 has been employed. Solving for b x ð Þ from (10), the continuous RK approximation is obtained as where W x;x 0 À x ð Þis the nth order continuous reproducing kernel function with the dependency on a and n implied, and is the so-called moment matrix. In arriving at (11) the symmetry of M was employed, although it is not necessary to construct the approximation. In addition to nth order accuracy, the approximation can be shown to satisfy the following equivalent reproducing conditions: which is often employed as the condition (for better conditioning of the moment matrix): As previously mentioned, it can be seen that the reproducing kernel approximation (11) can also be obtained by directly imposing the reproducing conditions (14) with the construction in (4).

Continuous implicit gradient
An implicit gradient directly approximates the derivative o a u x ð Þ of a function u x ð Þ for some fixed a in the same form of the reproducing kernel [55,57]: where the notation of multi-indices with parenthesis ðaÞ is introduced to indicate evaluation with a fixed value of a and to distinguish between terms of the form x a . Substitution of the Taylor series expansion in (8) into (15) yields The condition to reproduce the gradient o a u up to nth order accuracy can be expressed as all moments vanishing except the moment corresponding to a, that is, requiring D ðaÞ ½n u x ð Þ f g ¼ o a uðxÞ for some given a, which can be expressed as: where H ðaÞ is a column vector of fa!d ab g n jbj¼0 (emanating from J À1 ): To make the multi-index notation clear, consider H x ð Þ and H ðaÞ in the construction of D ðaÞ ½n fuðxÞg in two dimensions with n = 2, and a ¼ ð1; 0Þ, that is, to approximate first order derivatives with respect to x 1 with second-order accuracy: Solving for b ðaÞ from (17), the continuous implicit gradient approximation for o a uðxÞ is obtained as where again the symmetry of M has been employed. From the above, it can be seen that the RK approximation (11) can be considered a special case of (20) with jaj ¼ 0, i.e., W ð0;0;0Þ ¼ W which has been observed in the early history of meshfree methods [57,59]. This also demonstrates that since a simple change of H 0 ð Þ to H ðaÞ can approximate derivatives, the matrix M x ð Þ contains information about derivatives as well as the function itself, and provides an efficient way to obtain derivative approximations rather than direct differentiation of (11), for which the cost is not trivial [35]. This fact has been leveraged for solving partial differential equations more efficiently than using direct differentiation [36,60].

Computational Particle Mechanics
In addition, the implicit gradient approximation has been utilized for strain regularization to avoid ambiguous boundary conditions [55], avoid differentiation in stabilization for convection dominated problems [58], among other applications [60], and historically, the implicit gradient has in fact been widely used (it its discrete form) to solve partial differential equations: the generalized finite difference method [61], synchronized derivatives [57], as well as the pioneering work of the diffuse element method [2] all utilize approximations which are essentially coincident with the implicit gradient (for details, see [1]). Finally as will be discussed in Sect. 3.3, the so-called peridynamic differential operator [56] is also the implicit gradient approximation with the selection of the same bases and weighting functions.
Analogous to the RK approximation, the implicit gradient can be shown to satisfy the following gradient reproducing conditions [55]: Or equivalently, Z Similar to the reproducing kernel approximation, it can be seen that the implicit gradient approximation (20) can be also obtained by directly imposing (22) on the approximation (15).

Equivalence of the peridynamic differential operator and the implicit gradient
In this section, the implicit gradient approximation reviewed in Sect. 3.2, and the peridynamic differential operator introduced in [56] are compared. Let us start by recasting the implicit gradient from (15) as: whereÛ a represents the corrected kernel function: and H a is the kernel support, i.e., the portion of the domain where U a ðx 0 À xÞ is non-zero. Note thatĤ ¼ H in (15), that is, monomial bases of order n are employed, while herê H is used to indicate that a generic basis vector can be employed.
Directly imposing the condition of gradient reproduction (22) on (23) can be written as Z which is, for 0 b j j n, Z H a Also, from (26) we get the following system: Z or which, when we select the basis vectorĤ to be H, becomes the moment matrix employed in Sects. 3.1 and 3.2, and leads to the construction in (20). Consider now the peridynamic differential operator introduced in [56]. A Taylor expansion of a function f ðxÞ is first considered: . . .
where n ¼ x 0 À x and the remainder Rðn; xÞ is considered negligible.
The main idea of the peridynamic differential operator is to define orthogonal functions g p 1 p 2 ...p d n n ð Þ, where p i , is akin to a i (i.e., p and a are the same), and represents the order of differentiation with respect to x i , with i ¼ 1; . . .; d, such that the following gradient reproducing conditions are imposed: where and H x is the peridynamic neighbourhood of x. By comparing (23) of the implicit gradient and (30) of the peridynamic differential operator we see that they both Computational Particle Mechanics aim at reproducing derivatives of a function of x through a convolution by finding an appropriate kernel function. In fact: • H x in (30) is the same as H a in (23), as long as the RK kernel support and peridynamic neighbourhood coincide • f ðxÞ in (30) is uðxÞ in (23): they both represent a generic function of x • g p 1 p 2 ...p d n ðx 0 À xÞ in (30) plays the role ofÛ a ðx; x 0 À xÞ in (23).
The functions g p 1 p 2 ...p d n ðnÞ are found by imposing satisfaction of the following orthogonality property: that is, using the notation introduced in Sect. 3 whereñ ¼ ðn 1 ; n 2 ; . . .; n d Þ: This is the same condition that is imposed on the corrected kernelÛ a ðx; x 0 À xÞ of the RK implicit gradient (26). Therefore, ifÛ a x; Þ then the reproducing kernel implicit gradient and the peridynamic differential operator coincide. In [56] g p n x 0 À x ð Þis defined as: which can be rewritten as: and, for a given p, a p ð Þ is a column vector of unknown coefficients fa ðpÞ q g n q j j¼0 . For example, in two dimensions (d ¼ 2): By comparing (34) and (24) we notice that the definitions ofÛ a x; x 0 À x ð Þand g p 1 p 2 ...p d n x 0 À x ð Þ are analogous: both are the product of a weighted basis vector (34), respectively) and unknown coefficients to be determined (i.e., b ðaÞ x ð Þ and a p ð Þ ). It is therefore clear that, if the same weighted basis vectors are selected, the two are the same, meaning that the RK implicit gradient and the peridynamic differential operator are the same operator. For example, to select the same weighted basis vectors: • Due to the arbitrariness of the choice of basis, one can , theñ Þ, the same weighted basis vector is employed. Now, the unknown coefficients of the peridynamic differential operator are found by substituting the definition of which, for some given p and for 0 ñ j j n, leads to or [56]. Again, we can see that for a given gradient to be reproduced, if the same weighted basis vector is chosen (i.e., (38) and (39), and (27) and (28) are the same, respectively, and this leads (39) is for reproducing the gradient o p f for a given p. Since A is independent of p, the equations associated with reproducing gradients o p f for 0 p j j n can be combined [56]. For example, in the two-dimensional case we can write In conclusion we have shown, even though the peridynamic differential operator represents arbitrary derivatives of functions through a convolution over the peridynamic neighbourhood, while the RK implicit gradient performs it over the RK kernel support, the two operators and the idea behind them, i.e., correcting a convolution operator using a weighted basis vector to obtain gradients, are the same.

Deformation gradient under continuous implicit gradients
We now return to the analysis and comparison of the implicit gradient and the way in which state-based peridynamics constructs the non-local deformation gradient. In the Lagrangian RK approximation [62], shape functions are constructed with reference to the material coordinate X, and the approximation to displacement u h is constructed as: For simplicity, the dependence of these constructions on t will be implied for other expressions. The deformation gradient F is the gradient of the motion of the body x ¼ uðx,tÞ and is constructed as where I is the identity tensor. If implicit gradients (20) are employed under the Lagrangian formulation, the deformation is approximated in the material coordinate as . . .; 0 corresponds to the case of H ðaÞ with jaj ¼ 1, for approximating first order derivatives with respect to X j , i.e.: From the derivation by the Taylor expansion in Sect. 3.2, it can be inferred (or directly shown) that the deformation gradient constructed by implicit gradients possess nth order accuracy (or nth order consistency) without additional analysis needed.

Deformation gradient under continuous peridynamics
To relate the theory of peridynamics to classical continuum mechanics and provide the ability to employ conventional constitutive models, a principle called reduction can be employed [20] to relate the kinematic entity of the state Y to a non-local version of a deformation gradient F , which yields where K is the reference shape tensor that describes the undeformed configuration around the point X, and S is the deformed shape tensor which describes the deformed configuration around the point X. In the above, the function w d ðX 0 À XÞ is called the influence function, which has compact support with measure d. Thus the influence function in the non-local deformation gradient plays a closely analogous role to the kernel function in the construction of the deformation gradient by the implicit gradient (45), since both control the locality of the approximation of deformation. With F in hand, the associated stress is calculated as in classical continuum mechanics, and is then related to the force state T, which will be discussed in Sect. 4.3.

Analysis of the continuous peridynamic deformation gradient
We next examine the properties of the nonlocal deformation gradient F . First, considering that the displacement u ¼ x À X, using the definition of the deformation state Y X 0 À X h i¼ x 0 À x and the reference position state X X 0 À X h i¼ X 0 À X, the deformed shape tensor S can be expressed as: Computational Particle Mechanics Using (50), the non-local deformation gradient (47) can then be expressed as In component form, the above is, A Taylor expansion on u i ðX 0 Þ gives, after some algebra, Since the fourth term is an even function about X for symmetric functions w d , we have the following truncation error for the continuous non-local deformation gradient: When the neighborhood H X and influence function w d ðX 0 À XÞ are centred around X and symmetric about each axis (as in the case of spherically-shaped influence functions that are purely a function of X 0 À X j j), away from the influence of the boundary of the domain (maintaining a perfectly spherical neighborhood) the third term on the right-hand-side vanishes when integrated, since it is then an odd function centred around X and thus and the continuous form of the non-local deformation gradient in peridynamics is second-order accurate.
Near the boundary of the domain, or in the case that the shape of H X or the influence function w d ðX 0 À XÞ is not symmetric about each axis, the third term does not vanish, and we have: and the continuous non-local deformation gradient is firstorder accurate. Thus in finite domains where the neighborhood is not symmetric near the boundary, and in a general case of arbitrary influence functions and neighborhood definitions, the continuous non-local deformation gradient is globally first-order accurate.

Comparison between continuous implicit gradients and peridynamics
In [27], it was shown that in the interior of a domain (away from the boundary, or in an infinite domain) with uniformly discretized state-based peridynamics using the correspondence principle, the discretized non-local deformation gradient (47) is equivalent to employing a local deformation gradient by the discretized form of implicit gradients (45). However, it will be demonstrated that in the general case, this is not true for both the discretized form and continuous form.
In order to facilitate a more general comparison between the non-local deformation gradient by peridynamics and the RK approximation, we first express the shape tensors (49) in matrix form: where PðXÞ ½ X 1 X 2 X 3 T can be considered a type of ''basis vector'' consisting of monomials of order one. When the shape tensor is expressed this way, it is immediately apparent that it is coincident with the Lagrangian RK moment matrix with a ¼ d, U a X À X 0 ð Þ¼w d X À X 0 ð Þ, and linear basis, but omitting the unity term in the RK basis vector HðXÞ. More discussion on this point will follow later in the text.
The deformed shape tensor can also be expressed in matrix form as Noting that PðxðX 0 Þ À xðXÞÞ ¼ PðX 0 À XÞ þ uðX 0 ÞÀ ð uðXÞÞ we have Introducing a vector P r j ¼ ½d 1j ; d 2j ; d 3j and using the symmetry of K, the expression (60) can be recast in indicial notation as: For comparison, the local deformation gradient F ij calculated by implicit gradients (45) Þcan be expressed as: Thus, the implicit gradient can be viewed as a type of nonlocal operation with length-scale a, which is not surprising since it approximates differentiation by integration, just like the non-local deformation gradient in peridynamics (see [27] for additional discussions). Two key differences can be observed however. One is that the ''basis'' in peridynamics PðXÞ omits the unity term in HðXÞ. If PðXÞ were to be employed in the implicit gradient approximation (20), partition of nullity [the completeness condition for 0th order accuracy in (13)] would not be able to be satisfied. This fact however seems to be ''compensated for'' in the peridynamic gradient by the convolution with u i ðX 0 Þ À u i ðXÞ rather than u i ðX 0 Þ alone. That is, if u ¼ constant then the non-local deformation gradient (61) yields the correct result of F ij ¼ I ij (0th order accuracy). Thus it can be seen that the form (61) is inherently firstorder accurate, and in special cases, as has been demonstrated, is second-order accurate.
Another interesting point is that in examining (61), the non-local peridynamic calculation of a gradient uses values of u near X, except the actual value at X, while the implicit gradient still uses the value of u at X. Thus one could interpret the peridynamic operation ''more non-local'' versus the implicit gradient approximation. Indeed, when X 0 ¼ X in (61), the peridynamic ''kernel'' in the convolution also vanishes since P only contains first-order monomials.
Finally, it should be emphasised that the deformation gradient by implicit gradients, and the non-local deformation gradient by peridynamics, are clearly not the same. In [27] however, an equivalence was established in the special case of a uniform discretization, and away from the influence of the boundary.

Continuous reproducing kernel peridynamic approximation
In this section, the continuous reproducing kernel peridynamic approximation is presented, which unifies the way in which state-based peridynamics under correspondence approximates gradients, and the implicit gradient approximation. The unification also provides two other distinct cases which will be discussed.

Continuous reproducing kernel peridynamic approximation
The convolution operations for the approximation of the gradient of a function in (61) and (62) can be unified as follows. First, consider a kernel estimate of the type (15) with a basis of monomials from order m to n to estimate gradients of a scalar field u X ð Þ, with the convolution of uðX 0 Þ À uðXÞ rather than uðX 0 Þ as in (15), and a general weighting function x l with measure l: where Q ½m;n ðXÞ is a column vector of the set of monomials fX b g n jbj¼m , and the dependency of the operator on m and l is implied for notational simplicity. To facilitate nth order accuracy in this approximation, taking the Taylor expansion on uðX 0 Þ in (7) yields: where DðXÞ is a row vector of fo b uðXÞg n jbj¼1 and J is a diagonal matrix with entries f1=b!g n b j j¼1 . Examining (64), it Computational Particle Mechanics is apparent that in order to reproduce gradients D ðaÞ ½n fuðXÞg ¼ o a uðXÞ up to nth order accuracy, we have the following vanishing moment conditions: where Q ðaÞ ½m;n is a column vector of fa!d ab g n jbj¼m : If m = 1, then the system in (65) has a unique solution. Alternatively, if m = 0 and n [ 0 the system is underdetermined and an additional condition is required for Consider also imposing the partition of nullity on (63) in the case of m = 0 and n [ 0, that is, in addition to (65), the following is imposed: The system (65) can then be recast with (67) in hand to yield a determined system for all free variables: where A unified approximation is finally obtained by solving for b ðaÞ from (68) and substituting into (63): The approximation in (70) for derivatives is termed the continuous reproducing kernel peridynamic approximation herein. The selection of possible values of x l , jaj, n and m, yield the implicit gradient approximation, as well as the manner in which the deformation gradient is approximated by peridynamics. In addition, two more approximations can be obtained which are termed the continuous nth order non-local deformation gradient, and continuous nth order non-local higher-order derivatives, enumerated as follows: 1. When jaj [ 0, m = 0, x l ¼ U a and n is a free variable, (70) yields the implicit gradient approximation (20), since Q ½0;n ðxÞ ¼ HðxÞ and Q ðaÞ ½0;n ðxÞ ¼ H ðaÞ , and further if m = 0, then partition of nullity is satisfied by the imposition of (67) in (68), and the approximation (70) yields: 2. When m = 1, n = 1, and jaj ¼ 1 (approximating first order derivatives only), and choosing x l ¼ w d , (70) yields the peridynamic deformation gradient (61) 3. When m = 1, jaj ¼ 1, x l ¼ w d , and n [ 1 is a free parameter, when taking the gradient of the displacement u i with respect to X j , (70) yields a nth order accurate non-local deformation gradient which is denoted F ½n herein: where In the above, M ½1;n can be interpreted as a highorder reference shape tensor, while Q ½1;n can also be understood in the context of states. The arbitrarily high-order versions of non-local deformation gradients can also be understood in terms of reduction and expansion of states to tensors, and tensors to states, respectively. Details are given in ''Appendix A''. 4. When m = 1, jaj [ 1, and n is a free variable, choosing x l ¼ w d in (70) yields nth order accurate non-local higher order derivatives: The terminology adopted herein is that when m = 1 and x l ¼ w d , the derivative approximations are termed nonlocal or peridynamic since they approximate derivatives in the same manner as the non-local deformation gradient, they embed the non-local length scale d, and also perform differentiation by integration. The generalization occurs with jaj [ 1 and/or n [ 1. That is, the original non-local deformation gradient is recovered in the non-local derivatives when jaj ¼ 1 and n ¼ 1. Thus the general expression for nth order accurate non-local derivatives in the continuous case is D ðaÞ Therefore, the present formulation can be regarded as a generalization of the way in which state-based peridynamics approximates derivatives using the non-local technique.

Continuous nth order non-local gradient and divergence operations
High-order non-local gradient and divergence operations can be derived from the general formulation (70) with m = 1, x l ¼ w d as in (77). Examining (70) with jaj ¼ 1 (for first order derivatives) and casting it as an operator to approximate a non-local derivative of a function u with respect to X j , one obtains: where Q r j is the same vector in (75). In vector form, (78) can be cast as wherẽ and here ru is a column vector. Thus the nth order nonlocal gradient operation on a vector field f can be expressed as: where on the right hand side, f is represented as a column vector. Likewise, the nth order non-local divergence on a vector field f can be expressed as: 4.3 High-order force density via non-local divergence of stress: continuous case According to the correspondence principle [20], the force state T can be calculated from the 1st PK stress r as a function of the non-local deformation gradient (61) for state-based peridynamics as: In the above, and in the following text, r denotes the matrix form of the 1st PK stress according to the context.
Comparing (1), (2) and (83), the integration of the force state as a function of stress can be interpreted as a type of non-local divergence operation on the nominal stress r T as: where the fact that Þwas employed. A force state consistent with the high order non-local deformation gradient (73) can be obtained using the same correspondence principle as (see ''Appendix B'' for the derivation): The corresponding non-local divergence operation that results from this force state can be expressed as: However, despite being derived from the high-order accurate deformation gradient, numerical testing shows that the discretized form of the force density (86) viewed as a mathematical operator on r T does not guarantee even 0th order accuracy in the general case. That is, when the nominal stress r T is constant, the total force density contribution to a point X may be non-zero. Notably, this is also true of (84), where under a constant state of stress, the corresponding force-density could be non-zero, and the original formulation using (84) also does not possess even 0th order accuracy in the general case. Thus, despite whatever the accuracy the non-local deformation gradient possesses (using the original formulation (47), with linear accuracy, or the high-order formulation (73), with nth order accuracy), the final solution computed using either (84) or (86) may not even have 0th accuracy because of the operation of computing the force density. This assertion is verified numerically in Sect. 6.5. Accordingly, a non-local divergence of the nominal stress r T is proposed, which can be computed from the generalization of (82) from vectors to tensors as For the linear case (n = 1), the force density computed by the proposed non-local divergence can be expressed as: where the fact that P X 0 À X ð Þ¼ÀP X À X 0 ð Þ was employed. Comparing to the standard force density by peridynamics (84), which can also be simplified as it can be seen that constant and linear accuracy can be introduced easily into state-based peridynamics via a small modification of (89) to (88). In Sect. 6.5, it will be demonstrated that any order of accuracy desired is maintained in the discrete case.

Discrete reproducing kernel approximation and peridynamic deformation gradient
In this section, the discretized versions of the continuous approximations discussed in Sect. 3 are analysed and compared.

Discrete reproducing kernel approximation
A discrete version of the reproducing kernel approximation (4) can be obtained by performing numerical integration on both (11) and (12) at a set of NP nodes x J jx J 2 X f g NP J¼1 that discretize a domain X: Computational Particle Mechanics where u J uðx J Þ are nodal coefficients, W J ðxÞ is the reproducing kernel shape function, V J is the nodal quadrature weight for point x J , and MðxÞ is the discrete moment matrix. For notational simplicity, it should be understood that MðxÞ and other quantities denote the discrete or continuous counterparts depending on the context in which they are employed. Unlike the continuous case, the moment matrix (91) is conditionally invertible, which requires ðn þ d)!=ðn!d!) nodal kernels covering x which are non-colinear (in 2D) or non-coplanar (in 3D) [29]. The selection of a kernel value of a ¼ hðn + 1) with h the nodal spacing generally suffices as a rule of thumb.
It is important to note that when the quadrature in (90) is the same as in (91), nth order accuracy is maintained [5]. The construction in (90)-(91) can also be derived by using the Taylor expansion procedure in Sect. 3 (the derivation is omitted here), or can be derived by enforcing the reproducing conditions directly on W J , both of which can also demonstrate that the discretized form possesses nth order accuracy. The interested reader is referred to the literature for these procedures, e.g., [1], and for discussions on quadrature, see [5].
The set of shape functions W J x ð Þ satisfy the so-called nth order reproducing conditions, i.e., possess nth order completeness: or equivalently, as it is often expressed and employed for better conditioning of the moment matrix: Alternative to the construction in (90)-(91), determination of quadrature weights may be avoided by constructing the so-called discrete RK approximation [62], which directly imposes (92) on a corrected kernel function.

Discrete implicit gradient approximation
Analogous to the case of the discrete RK approximation (90), a discrete implicit gradient approximation can be obtained by employing quadrature on (20) at nodal locations: where MðxÞ is the same discrete moment matrix in (91). It is apparent that as in the continuous case, as a special case of (94) with jaj ¼ 0 we obtain (90), that is ð Þ. So long as the moment matrix is discretized with the same quadrature in (94), the gradients also enjoy nth order accuracy, which can be confirmed by deriving (94) from a Taylor expansion point of view (again the derivation is omitted here), as in the RK approximation, or by directly enforcing gradient reproducing conditions on W ðaÞ J . It can be shown that the derivative approximations (94) enjoy gradient completeness [55]: Or again, analogous to (22), If desired, quadrature weights in implicit gradients may also be avoided by employing a discrete implicit gradient approximation [55], which directly imposes (95) on a corrected kernel. One final point to note which will be revisited, is that given a discrete set of scattered data fu J g NP J¼1 , an approximation to derivatives o a u can be obtained at any given point x of interest using the implicit gradient approximation, providing a smooth field of derivative estimates in the entire domain.

Deformation gradient under discrete implicit gradient approximation
The discretization of the deformation gradient under the implicit gradient approximation (45) is calculated in the Lagrangian coordinates in a similar fashion as the continuous case: where H r j is the same vector in (46).

Computational Particle Mechanics
Since the deformation gradient is constructed by discrete implicit gradients with nth order accuracy (or nth order consistency), it can again be directly inferred that the discrete deformation gradient (97) possesses nth order accuracy without additional analysis needed.

The discrete deformation gradient under peridynamics
To discretize (47)- (49), nodal quadrature at NP nodes is employed in the meshfree implementation [24], similar to the discrete RK approximation and discrete implicit gradient approximation: The discretized non-local deformation gradient is calculated from the above quantities as Again, for notational simplicity, depending on the context it should be understood whether K, S or F is the continuous or discrete form.

Analysis of the discrete deformation gradient under peridynamics
Following the procedures in the continuous case in Sect. 3.6, a Taylor expansion on the displacement in (98)-(99) obtained from xðXÞ ¼ X þ uðXÞ yields the discrete expression: In order to interpret the implications of (100), first consider the case of a uniform discretization, away from the boundary, with symmetric influence functions. The fourth term on the right hand side is an even function (centred around X) and will be non-zero due to the symmetry and constant nodal weights. Under these conditions, the following is obtained: In the same situation, the third term disappears, but this time only at nodal locations due to the discrete quadrature, as the terms being summed will only cancel when they are ''anti-symmetric'' (which only occurs at nodal locations): This also indicates that in the general case, it will be nonzero unless a careful selection of the combination of influence functions, nodal quadrature weights, and so on, are selected, as was performed in [52], but is difficult to generalize to non-uniform discretizations. Finally, in the case of a non-uniform discretizations, or even in a uniform discretization near the boundary, and also away from nodal locations, the third term in (100) will not disappear and the expression reduces to Thus the discrete form and the continuous form share the similar order of accuracy and behavior in accuracy; in the best-case they are both second order accurate, and in the general case, they are first-order accurate.
In summary, the same order of accuracy is attained for both continuous and discrete versions of implicit gradients, yet the continuous and discrete versions of the non-local deformation gradient in peridynamics slightly differ. That is, the accuracy is second-order in the best case (away from the boundary, uniform discretizations, symmetric influence functions), but in the discrete case second-order accuracy can only be obtained at the nodes in this situation. In the general case, the constructions are first-order accurate, for both integral and discrete forms. These situations can be rectified with the proposed generalized discrete formulation given in Sect. 6.1.

Comparison between discrete implicit gradients and peridynamics
In order to facilitate a comparison between the discrete non-local deformation gradient by peridynamics and the discrete RK approximation, we first express the discrete shape tensors (98) in matrix form: As before, the undeformed shape tensor is coincident with the Lagrangian RK discrete moment matrix (91) with a ¼ d, U a ¼ w d , and linear basis, but omitting the unity term in the vector HðXÞ.
Following the procedures in Sect. 3.7, the discrete nonlocal deformation gradient can be rearranged and expressed as: For comparison, the local deformation gradient F ij calculated by discrete implicit gradients (94) with a ¼ d, and U a ¼ w d can be expressed as: It can be seen that in the discrete case, if PðXÞ were to be employed in the implicit gradient approximation, partition of nullity would also not be able to be satisfied. This fact is again ''compensated for'' in the discrete peridynamic gradient by the summation with u i ðX J Þ À u i ðXÞ ½ rather than u i ðX J Þ alone. That is, if u ¼ constant then the discrete nonlocal deformation gradient still yields the correct result of Another interesting point is that the gradient approximation (105) does not allow a continuous gradient field representation from a discrete set of nodal data, since in computing the quantity away from the discrete points with known solutions (e.g., from scattered data or a PDE) u i ðX J Þ, the quantity u i ðXÞ is unknown. This places a serious drawback on the approximation, as it can only yield gradient estimations at the scattered data points themselves, but does not provide for an interpolation function for the data. Thus formally, the peridynamic approximation cannot provide a smooth field at all points X in the domain given a finite set of nodal coefficients fu J g NP J¼1 . It can however, given a function uðXÞ defined in the entire domain, provide an estimate of derivatives. Finally, it can be noted that it would be possible to interpolate the derivative estimates at nodes, although this yields some additional complexity.

Reproducing kernel peridynamic approximation
In this section, the discrete form of the reproducing kernel peridynamic approximation is given. High-order non-local discrete deformation gradients and non-local divergence operations are derived, as well as several other discrete approximations. The order of accuracy of these operators is also verified numerically.

Reproducing kernel peridynamic approximation
Similar to the continuous case, the operations for the discrete approximation of the gradient by peridynamics and implicit gradients can be unified as follows. First, consider a discrete approximation of the type (94) with a basis of monomials from order m to n to estimate gradients of a scalar field u X ð Þ, with use of uðX J Þ À uðXÞ rather than uðX J Þ as in (105) where Q ½m;n ðXÞ is the same column vector of the set of monomials fX b g n jbj¼m as the continuous case. To facilitate nth order accuracy in this approximation, taking the Taylor expansion on u X J ð Þ in (7) yields: where DðXÞ and J are the same vectors and matrices in the continuous case (64). Examining (108), it is apparent that in order to reproduce gradients up to nth order accuracy, we have the following discrete vanishing moment conditions: where Q ðaÞ ½m;n is a again a column vector of fa!d ab g n jbj¼m . As before, when m = 0 and n [ 0 the system is underdetermined and an additional condition is required for determining b ðaÞ . Imposing the partition of nullity on (107) in the case of m = 0 and n [ 0: the system (109) can then be recast with (110) in hand to yield a determined system for all free variables: where A unified discrete approximation is obtained by solving for b ðaÞ from (111) and substituting into (107): The approximation in (113) for derivatives is termed the reproducing kernel peridynamic approximation herein. The selection of possible values of x l , jaj, n and m, yield both the discrete implicit gradient approximation, as well as the manner in which the deformation gradient is approximated by the discretized version of peridynamics, and as before, two additional approximations can be obtained which are termed the nth order non-local deformation gradient, and nth order non-local higher order derivatives herein: 1. When jaj [ 0, m = 0, x l ¼ U a , and n is a free variable, (113) yields the discrete implicit gradient approximation (94), since Q ½0;n ðxÞ ¼ HðxÞ and Q ðaÞ ½0;n ðxÞ ¼ H ðaÞ , and further if m = 0, then the discrete partition of nullity is satisfied by the imposition of (110), and the approximation yields: 2. When m = 1, n = 1, 3. When m = 1, jaj ¼ 1, x l ¼ w d , and n [ 1 is a free parameter, when taking the gradient of the displacement u i , (113) yields a nth order accurate non-local deformation gradient denoted F ½n herein: and Q r j is the vector in (75). 4. When m = 1, jaj [ 1, and n is a free variable, choosing x l ¼ w d , (113) yields nth order accurate non-local higher order derivatives: Computational Particle Mechanics As before, a general expression for discrete nth order accurate non-local derivatives can be found by setting m = 1 and It should be noted that similar to the discrete deformation gradient in (99), this approximation cannot yield a continuous field of derivative estimates from a set of data fu J g NP J¼1 , since u X ð Þ would be required at other points aside from the nodal positions. This again places a limitation on the approximation. However derivative approximations can be obtained at the nodal positions X I themselves as: with u J u X J ð Þ.
In contrast, the implicit gradient approximation (94), as well as the direct derivative of the RK approximation (90) can yield an approximation to derivatives of a field at every point, constructed from a set of scattered nodal data fu J g NP J¼1 . For the construction in (119), a form of postprocessing or reinterpolation would be involved.
Finally, it is notable that a discrete version of (113) can be employed that omits quadrature weights entirely. Following the procedures common in the RK literature [1], the following expression for the discrete reproducing kernel peridynamic approximation without quadrature can be shown to maintain nth order accuracy:

Deformation gradient test
Here the accuracy of the non-local deformation gradient by the standard technique (99) and the proposed high-order technique (116) is assessed and verified. The following displacement fields are considered, with each component (i = 1,2, two dimensions) prescribed the same value: 1. First order: u Here it can be seen that in uniform discretizations, the original non-local deformation gradient has second-order accuracy away from the boundary, but only first order accuracy in the presence of the influence of the boundary. For the non-uniform discretization, the deformation gradient is first-order in all cases. This confirms the analysis and discussions in Sect. 5.5. In contrast, the proposed highorder formulation yields the desired order of accuracy in all cases, in both uniform and non-uniform discretizations, as well as in the presence of a boundary, and verifies the ability of the formulation to provide arbitrary-order accuracy in the discrete case.
Computational Particle Mechanics

nth order non-local gradient and divergence operations
High-order non-local gradient and divergence operations can be derived from (113) analogous to the continuous case. Setting m = 1, jaj ¼ 1 and x l ¼ w d in (113) and casting it as an operator for the non-local first order derivative with respect to X j , one obtains: where Q r j is the same vector in (75). In column vector form, (123) can be expressed as whereQ r is the same matrix in (80).
it can be seen that constant and linear accuracy (confirmed in the discrete case now, since the proposed operator is inherently nth order in the discrete case as well) can be introduced easily into the existing peridynamics formulation via a small modification.

Force density test
Consider again the discretizations used in Sect. 6.2, with the uniform and non-uniform node distributions shown in Fig. 3. The interior and boundary points are again considered for testing the discrete force density by peridynamics (127) and the proposed non-local divergence of stress (129). Stresses are computed from the displacements fields listed in Sect. 6.2, yielding one order lower polynomial in the stress solution. The Lamé constants for this computation are chosen as Young's modulus E = 100 and Poisson's ratio t = 0.3. The force density computed from the standard peridynamic technique is denotedT PD i , while the force density from the high-order non-local divergence operation is denotedT HOPD i . The order chosen in the higher-order accurate formulation is chosen to be one order higher than the stress, corresponding to the case that the same order of accuracy is used in both the deformation gradient and the force density calculation. Tables 3 and 4 show the absolute error in the solution at the two points of interest, computed using these two methods for various orders of solutions, both for uniform and non-uniform discretizations, respectively. It can be seen that in uniform discretizations, the traditional technique can compute the correct solution in the interior of the domain for stress fields up to second-order. On the other hand, near the boundary, this method does not produce even 0th order accuracy. That is, when the stress is constant, the force density will be non-zero. For non-uniform discretizations, the traditional approach does not ever yield 0th order accuracy.
In contrast, the proposed formulation can compute the exact force density for constant, linear, and quadratic stress fields, regardless of the nodal arrangement (uniform, nonuniform), and also in the presence of a boundary. These numerical tests verify the necessity of the introduction of the non-local divergence operation introduced in Sect. 6.4.
The numerical examples in Sect. 8 further demonstrate the consequences of the choice of formulation.

The reproducing kernel peridynamic method
In this section, a short summary of the proposed reproducing kernel peridynamic method is given. First, the discrete high-order deformation gradient is given in compact form, with some examples of implementation. The high-order non-local divergence of stress is then introduced, and finally, the node-based collocation formulation is given.

High-order non-local deformation gradient
The high-order accurate non-local deformation gradient in (116) can be expressed succinctly at a node X I as: Or, in matrix form, (132) can be expressed as: with To take an example of the terms involved, selecting n = 2, one obtains the following: where we recall here for convinceQ r ¼ ½Q r 1 ; Q r 2 ; Q r 3 .

High-order non-local divergence of stress
With (132) in hand, classical techniques can be employed to compute the 1st PK stress r at each node. The non-local divergence of the nominal stress (r T ) can then be computed at node X I as: or, simply, where U IJ is the same term in (135). Finally, the non-local equation of motion is solved for the proposed formation under the nodal collocation framework, for all nodes X I which are not boundary nodes: Procedures for construction of an elastic stiffness matrix based on peridynamic-type approximations (134) and (140) can be found in [63]. For dynamic problems, standard time integration techniques can be employed to solve the semidiscrete Eq. (141).

Numerical examples
In this section, the order of solution exactness of the collocation method using (134) and (140), described in Sect. 7, and the existing peridynamic method using (99) and (127) are tested, along with their associated convergence rates, with the resulting formulations denoted as PD for standard state-based peridynamics, and RKPD for the reproducing kernel peridynamic method, respectively. Permutations of their operators are also tested to examine the effect. Nodal quadrature is employed in all approximations of integrals.
In order to enforce essential boundary conditions, a ghost boundary layer is considered in the examples, along with direct enforcement on boundary nodes without a ghost layer to examine the effect. For ghost nodes, a layer of uniform ghost nodes is generated with thickness based on the horizon, sufficient to eliminate effects of a finite boundary. All boundary conditions are pure Dirichlet, to set aside any complications with enforcing natural boundary conditions. Finally, in addition to other nodes, ghost nodes are also employed as collocation points.
As discussed previously, convergence in peridynamics can be interpreted in a few ways [48]. For one, a non-local solution can be obtained with d fixed, and one may examine the error as the nodal spacing goes to zero. The numerical solution in this case converges to the non-local solution [46,48]. On the other hand, as the non-local length scale in peridynamics goes to zero concurrently with the nodal spacing in the discretization, the numerical solution converges to the local solution [48]. In the following examples, the latter is chosen to be tested. Accordingly, horizon sizes with a fixed ratio to nodal spacing are chosen as 1.75, 2.5, and 3.5 for the linear, quadratic, and cubic order formulations, respectively. All problems are solved with influence functions chosen as the cubic B-spline in (5).

Patch tests
In the following set of patch tests, ''PD'' denotes the statebased peridynamic force density and deformation gradient, while ''RKPD'' indicates the proposed high-order peridynamic formulation for force density and deformation gradient. Permutations of these two formulations are tested in order to assess and verify the order of accuracy of the operators involved.
Consider a two-dimensional linear patch test, which requires recovering the exact solution by a numerical method when the solution to a boundary value problem is linear. Zero body force is prescribed with the following Dirichlet boundary conditions: With zero body force and the prescribed displacement in (142), the solution is the same linear displacement (142). Both the uniform and non-uniform discretizations shown in Fig. 4 are examined, in order to test the effect of non- Fig. 4 Nodal discretizations for patch tests: a uniform, b uniform with ghost nodes with two layers (used for the quadratic case), c nonuniform; and d nonuniform with ghost nodes with two layers (used for the quadratic case)

Rows of text in boldface indicates the method passes the patch test
Computational Particle Mechanics uniformity of node distributions. Here, the case with two layers of ghost nodes necessary to eliminate the boundary effect for the quadratic test is depicted for illustration. For the non-uniform case, the nodes in the uniform discretization are perturbed away from their original position in a random fashion, as shown in Fig. 4. Table 5 shows the results for the linear patch test under various conditions. Note that the deformation gradient with linear accuracy under RKPD is coincident with the standard peridynamic deformation gradient, although the force density is not. First, it can be seen that the original formulation can only pass the linear patch test under a uniform discretization with ghost boundary nodes. On the other hand, RKPD is able to pass the patch test without ghost nodes, and in both uniform and non-uniform discretizations.
The deformation gradient analysis and tests show that the standard peridynamic deformation gradient exhibits at least linear accuracy in all situations: for a linear displacement field, such as the one in this problem, the associated constant deformation gradient is calculated exactly. Thus in this problem, the peridynamic force density is operating on a constant stress field. In the case of a non-uniform discretization and/or without ghost nodes, the standard peridynamic method fails the patch test, indicating that the peridynamic force density, considered as a mathematical operator on the stress, does not possess 0th order accuracy in non-uniform discretizations or near the boundary of the domain, which is confirmed also by the tests made in Sect. 6.5. On the other hand, these results confirm that the force density is at least 0th order accurate in uniform discretizations and away from the boundary, since in this problem it operates on a constant stress, and the patch test is passed in a uniform discretization with ghost nodes. Now consider a quadratic patch test with the exact solution: The procedure to design such a higher-order patch test (boundary conditions and other prescribed data) has been described in various references (cf. [41]). Table 6 shows the results using both linear accuracy in the deformation gradient (PD, RKPD coincident) and quadratic accuracy using RKPD, under various conditions. Once again the standard peridynamic method exhibits quadratic exactness in uniform discretizations with ghost nodes, which agrees with the analysis and other numerical tests herein. However when one or more of these conditions is violated, the method fails to pass the patch test. The implications are discussed as follows.
Peridynamics in a uniform discretization with ghost nodes will calculate the exact linear deformation gradient for quadratic displacements; thus the peridynamic force density in this problem operates on a linear stress. And since the quadratic patch test is passed with peridynamics only in uniform discretizations with ghost nodes, this indicates that the peridynamic force density operation on the stress is at least first-order accurate in uniform discretizations, away from a boundary. And otherwise, the operator is again not even 0th order accurate, as in these cases the method fails the patch test.
Finally, the proposed formulation with second-order deformation gradient and second-order non-local divergence of stress can pass the quadratic patch test under both uniform and non-uniform discretizations, and does not need ghost nodes.
Next, a cubic patch test is considered. The boundary conditions and body force are prescribed according to the exact displacement solution: The results under the various conditions and permutations of formulations are shown in Table 7. With cubic accuracy in the deformation gradient, the exact quadratic deformation gradient is obtained under RKPD. Thus the standard peridynamic force density in this problem operates in some cases on a quadratic stress, and from the results in the patch test, the correct result can be obtained with uniform discretization and using a boundary layer. With this test, it is apparent that the original deformation gradient and the original force density can both yield quadratic accuracy in a uniform discretization away from the boundary, but in the general case, the force density operator does not possess even 0th order accuracy.
On the other hand, RKPD is seen to pass the patch test when both the deformation gradient and force density are calculated using the proposed formulation with cubic accuracy. Ghost nodes are not necessary, and it can pass the cubic patch test in a non-uniform discretization.
One final interesting note, is that formulations using two operators with quadratic accuracy can pass the cubic patch test in uniform discretizations with ghost nodes (i.e., the standard technique in uniform discretizations with ghost nodes, and RKPD). That is, when two operators of a certain order are used in conjunction with one another (such as a Computational Particle Mechanics deformation gradient paired with the force density here), the combined formulation possesses one order of accuracy higher, at least in uniform discretizations, away from the boundary. This observation is consistent with the results of the recursive gradient technique of using two combined operators [37].

Rows of text in boldface indicates the method passes the patch test
over the domain [-1, 1], with boundary conditions: This problem can be considered a one-dimensional elastic bar with unit Young's modulus, with a high-order body force term. The solution to this problem is u ¼ e x . Two cases are considered of uniform and non-uniform discretizations. The uniform case is discretized with nodal spacing h ¼ [1=4, 1=8, 1=16, 1=32], as shown in Fig. 5a. For the non-uniform case, shown in Fig. 5b, the first refinement of the uniform case is perturbed, with each subsequent refinement adding nodes at the halfway point between two nodes.
The first case tested is the standard peridynamic formulation (PD) and higher-order peridynamics (RKPD) with linear accuracy. In this case, the deformation gradient calculations are the same, but the force densities differ. The convergence of the solutions in the L 2 norm for uniform, non-uniform, ghost nodes, and no ghost nodes are shown in Fig. 6. Here it can be seen that the only case in which the standard peridynamic formulation converges is in the uniform case with ghost nodes. This could be attributed to the fact that in non-uniform discretizations, or in the presence of a finite boundary, the force density operation on the stress does not possess 0th order consistency. On the other hand, the proposed RKPD formulation converges at the rate of two (n ? 1) for both uniform and non-uniform discretizations, with and without ghost nodes. These results confirm the necessity of a high-order non-local divergence for a high-order numerical solution, that is, a high-order non-local deformation gradient alone is insufficient for high-order accuracy.
Finally, it should be noted that the rate of n ? 1 for this odd-order approximation indicates superconvergence of the solution. That is, in nodal collocation approaches, the rate of n -1 is observed for odd orders of accuracy [17,36,37,40], with the further implication that linear completeness does not yield convergence. Meanwhile, the present formulation can converge with linear accuracy. This could be attributed to the fact that in essence, the combination of the two linear operators here produce something akin to higher-order differentiation, which is similar to the recursive gradient technique recently introduced [37], where linear basis was shown to be sufficient for convergence in strong form collocation, while also exhibiting the superconvergence phenomenon. The second case considered is a quadratic-order deformation gradient paired with the standard force density, denoted ''PD'' and the high-order force density, denoted ''RKPD'', respectively, in order to test the effect of the different permutations, as well as the proposed RKPD formulation. Figure 7 shows the convergence plots for the cases tested. Here it can again be seen that using both the high-order deformation gradient and the high-order force density yields consistently convergent results. In addition, the proposed RKPD formulation can provide much lower errors than PD, and yields convergent solutions in both uniform and non-uniform cases tests, as well as with and without ghost nodes. On the other hand, the other cases either slowly converge, or do not converge at all. A slope of two (n) is indicated in the figures, where it can be seen that RKPD obtains rates consistent with nodal collocation of the strong form [17,[36][37][38]. That is, for even orders of accuracy, a rate of n is expected, and is obtained with the current formulation.
Finally, the cubic RKPD deformation gradient is tested, paired with the traditional force density, denoted ''PD'', and the high-order force density, denoted ''RKPD''. Figure 8 shows that the proposed formulation converges at the rate of n ? 1, while PD does not converge at all. This again indicates superconvergence for odd orders of approximation (two orders higher than n -1 for standard nodal collocation formulations). In addition, the high-order accuracy can be obtained in uniform and non-uniform discretizations, with and without ghost nodes.

Convergence in a manufactured 2-D solution
A two dimensional manufactured plain-strain elasticity problem is considered over a domain [-1, 1] 9 [-1, 1] with the following exact solution: Accordingly, the body force density components associated with this solution, given in components are where E is Young's modulus, with E = 100,000 and t is Poisson's ratio, with t = 0.3 in this example. The essential boundary conditions on all four edges of the domain are given by the displacement in (147). As before, uniform and non-uniform discretizations are considered, shown in Fig. 9, with and without ghost boundary nodes, resulting in four test cases. The ghost nodes are uniform and the thickness of the ghost layer is selected to be sufficient to eliminate any effect of the boundary on the approximations, as in the previous examples.
First, linear accuracy in the deformation gradient (PD, RKPD coincident) paired with both the standard force density (denoted PD) and the proposed high-order force density (denoted RKPD) is considered. Figure 10 shows the convergence of the solution for both uniform and nonuniform discretizations. PD without ghost nodes converges Computational Particle Mechanics at a rate of approximately one in the uniform case without ghost nodes, and approximately two in the case of a uniform discretization with ghost nodes. PD is seen to essentially not converge at all in the non-uniform case, likely due to the lack of consistency in the operation on stress to produce the force density in the standard formulation. Finally, it is again seen that the proposed formulation can consistently give convergent results across all types of discretizations, with the superconvergent rate of n ? 1, rather than n -1 in standard direct nodal collocation techniques.
A second-order case is considered next, with RKPD denoting second-order accurate operators, and PD denoting the second-order accurate deformation gradient paired with the standard force density. Figure 11 shows the convergence behaviour of the various permutations in discretizations and solution techniques. The trend is similar to the one-dimensional case. Here, it is seen that secondorder convergence is obtained (order n), which is consistent with nodal collocation approaches with even-order accuracy. For the PD formulation, the solution is again seen to converge slowly or not at all, and in the case of PD without ghost nodes, the solution is seen to diverge in non-uniform discretizations.
Finally, the cubic case is considered. RKPD again denotes the pairing of the high-order deformation gradient and non-local divergence, while PD denotes the high-order deformation gradient paired with the standard force density. Figure 12 shows that the proposed RKPD formulation converges in all situations with the superconvergent rate of n ? 1 (rather than n -1), both with or without ghost Computational Particle Mechanics nodes, and in uniform and non-uniform discretizations. Similar to other cases, the PD approach either converges incrementally, or not at all.

Conclusion
A unification of local and non-local meshfree approximations is presented and is termed the reproducing kernel peridynamic approximation. The continuous, or integral form is presented, as well as the discrete form. With the selection of free parameters, the generalized approximation can yield both the implicit gradient approximation, and the way in which state-based peridynamics under the correspondence principle approximates derivatives via the nonlocal deformation gradient. Perhaps just as important, the generalization yields the ability to obtain arbitrary-order accurate non-local gradient and divergence operations, as well as higher-order non-local derivatives, with arbitrary accuracy as well. Thus, the formulation generalizes and formalizes the concept of non-local derivatives beyond first-order.
The framework also demonstrates that the non-local deformation gradient is not in fact equivalent to the implicit gradient approach, which has been speculated previously. In addition, it has been shown that the peridynamic differential operator is equivalent to the implicit gradient approximation, which by analogy is not equivalent to the generalization of the non-local approximation of derivatives.
The analysis presented demonstrates that the continuous form of the non-local deformation gradient is at best second-order accurate, but in the general case near the influence of a boundary, is first-order. In the discrete case, this deformation gradient is also at best second-order accurate, in the case of uniform nodal distributions, away from the influence of the boundary, but only at nodal locations. In all other situations, it is first-order accurate. Several numerical examples are provided to confirm this analysis.
The force density in terms of the stress, obtained via the standard peridynamic formulation is also found to be at best, second-order accurate, and at worst, without any order of consistency, which has been confirmed numerically using several tests. This work proposes a high-order accurate non-local divergence operation on the stress to replace this force density, in order to obtain globally highorder accurate numerical solutions. The pairing of the highorder non-local deformation gradient, along with the highorder non-local divergence of stress is termed the reproducing kernel peridynamic (RKPD) method. The method is tested under the collocation framework, although weak formulations are also certainly possible. It is suspected however, that in addition to the computational expense of double integrals, numerical integration may play a key role on the convergence of the method, as in local Galerkin meshfree methods [1].
The numerical examples demonstrate that the existing peridynamic formulation can pass the linear, quadratic, and cubic patch tests, but only in uniform discretizations, with ghost nodes. In all other situations, the method fails to pass any patch test. On the other hand, the proposed formulation is able to pass arbitrary-order patch tests in both uniform and nonuniform discretizations, with and without ghost nodes.
The examples further show that the existing peridynamic formulation can generally converge at rates of two and one for uniform discretizations, with and without ghost nodes, respectively. The examples further demonstrate that in non-uniform discretizations, the standard deformationgradient based peridynamic formulation essentially does not converge at all. This situation can be rectified by the proposed formulation, which shows consistent convergence behavior in both uniform and non-uniform discretizations, with and without ghost nodes. One noteworthy aspect of this approach is that linear accuracy will still yield convergent results, in contrast to the standard collocation approach. In addition, the proposed formulation exhibits superconvergence: a rate of n is obtained for even orders of accuracy, and n ? 1 for odd orders. This is in contrast to typical nodal collocation approaches, where the rates observed are n and n -1 for even and odd orders, respectively.
So far, the convergence behavior has been tested with respect to local solutions. Quadrature weights were also employed in the approximation, which are not necessary to obtain nth order accuracy. In addition, all problems tested were with pure Dirichlet boundary conditions. Finally, the performance against some existing similar methods, such as collocation with implicit gradients, or collocation with explicit gradients, was not tested, although the present generalized formulation can encapsulate these methods. Future work could examine the convergence to non-local solutions, natural boundary conditions, the approximation accuracy of the non-local approximation with and without quadrature weights, comparisons with similar numerical formulations, as well as formulation of the RKPD method under the Galerkin framework.
Other future work could possibly leverage the unification presented for techniques that have been extensively developed for local methods over the past several decades, such as the use of enrichment functions. Finally, the higher order non-local derivatives that can be obtained by the formulation have not been investigated for any particular use yet. Two possibilities would be to employ the non-local derivatives to directly discretize partial differential equations, or employ them for strain regularization in order to avoid ambiguous boundary conditions. The direct discretization of derivatives however, is suspected to yield the typical collocation convergence rates, rather than the superconvergent rates of the present formulation.
One last point, the importance of which was emphasized by feedback from a talk on this method by the late Steve Attaway, is that under this unification, state-based peridynamic codes can be converted to use local meshfree approximation methods, while local meshfree codes can be converted to use non-local peridynamic approximations.
The nth order reduction of the deformation state Y yields an nth order non-local deformation gradient tensor F ½n : where S ½n ðXÞ is a high-order deformed shape tensor and K ½n ðXÞ is the high order reference shape tensor in (156) The high-order non-local deformation gradient (160) is identical to the one expressed in (73).
where e ijk is the Levi-Civita symbol.