Peridynamic Galerkin method: an attractive alternative to finite elements

This work presents a meshfree particle scheme designed for arbitrary deformations that possess the accuracy and properties of the Finite-Element-Method. The accuracy is maintained even with arbitrary particle distributions. Mesh-based methods mostly fail if requirements on the location of evaluation points are not satisfied. Hence, with this new scheme not only the range of loadings can be increased but also the pre-processing step can be facilitated compared to the FEM. The key to this new meshfree method lies in the fulfillment of essential requirements for spatial discretization schemes. The new approach is based on the correspondence theory of Peridynamics. Some modifications of this framework allows for a consistent and stable formulation. By applying the peridynamic differentiation concept, it is also shown that the equations of the correspondence theory can be derived from the weak form. Likewise, it is demonstrated that special moving least square shape functions possess the Kronecker-δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} property. Thus, Dirichlet boundary conditions can be directly applied. The positive performance of this new meshfree method, especially in comparison to the Finite-Element-Method, is shown in the calculation of several test cases. In order to guarantee a fair comparison enhanced finite element formulations are also used. The test cases include the patch test, an eigenmode analysis as well as the investigation of loadings in the context of large deformations.


Introduction
Meshfree methods have represented a big promise in computational mechanics for years. Due to the absence of a sorted connectivity and requirements for the distribution of points, even large deformations with free surfaces can be represented in a Lagrangian description. The simulation of crack formation or the fusion of materials in Additive Manufacturing can be facilitated with these methods. On the other hand, additional algorithms or stabilization schemes are often needed to achieve the required accuracy in solving initial boundary value problems [47]. Also for this reason, mesh-B C. Weißenfels christian.weissenfels@mrm.uni-augsburg.de 1 based methods, such as the Finite-Element-Method (FEM), are mainly used in engineering applications until today. The FEM convinces by an accurate solution as long as the underlying elements are not too distorted. This places additional demands on the algorithms to generate a proper FE-mesh. In addition, very large deformations cannot be represented in a Lagrangian description without constant remeshing. Therefore, the simulation of fluid flows is usually based on an Euler description of the differential equation. To represent free surfaces, additional algorithms, such as the volume-offluid method [19] or the level set approach [33] are often employed in this case.
Peridynamics [38] represents an alternative description of the balance laws in continuum mechanics motivated from molecular dynamics. In its discrete form, this concept can be classified as a meshfree particle method. The correspondence formulation allows the incorporation of continnuum mechanical laws into state-based Peridynamics [42]. However, this approach also has shortcomings. One lies in the observation of low energy modes. These negative effects occur in the FEM in the case of under-integration. To remedy this defi-ciency, a stabilization term is often added [27], [39]. This correction is also found in Smoothed Particle Hydrodynamics [15] or in the Optimal Transportation Meshfree (OTM) method [48]. The disadvantage of this concept is the use of an additional parameter whose optimal value is not known before the start of the simulation. Another possibility is to divide the neighborhood into individual smaller areas where the equations are evaluated separately [14]. As shown in [17], this partitioning can lead to incorrect results in some cases. In addition, other modifications can be found that rely on a higher order approach in the approximation of the deformation gradient [52] or introduce additional stress points to increase the number of evaluation points [30].
Another possibility, which does not require numerical parameters, additional evaluation points, or a higher order formulation, but still remedies the negative effects from under-integration, can be found in [6]. In contrast to the classical correspondence theory, the kinematic and kinetic quantities in the neighborhood are assumed to be nonconstant and evaluated at every point in the neighborhood. An extension for weakly compressible or incompressible material behavior can be found in [7].
Another negative aspect of the classical correspondence theory is the failure of the patch test even when the negative effects from under-integration have been removed. A suitable approach to overcome this deficiency in the peridynamic framework can be found in [8]. This correction also allows optimal convergence rates to be achieved.
The scheme in [6] and [8] is based on a Petrov Galerkin formulation. Thus, the resulting stiffness matrix is not symmetric even for conservative systems.
This work presents an approach that satisfies essential requirements for spatial discretization schemes. Moreover, the stiffness matrix is symmetric for conservative systems.
In addition, there exist approaches in which Peridynamics is related to FEM in order to exploit the advantages of both methods, especially in the approximation of fracture formations. In [31], the bond-based variant is directly integrated into a FEM software. In addition, algorithms for coupling both discretization schemes are proposed. In [13] as well as [36], Peridynamics is expressed in the weak form and classical FE shape functions are used for the approximation. A convergence analysis for this approach in the context of an ordinary state-based formulation can be found in [22].
In contrast, the approach presented here does not use FE shape functions. Likewise, a unique partitioning of the domain into elements can be omitted. More specifically, a guideline is presented on how to rewrite the equations of classical peridynamic correspondence theory to satisfy the necessary conditions for spatial discretization schemes. This allows the direct imposition of boundary conditions despite the non-local nature of Peridynamics. In contrast to the FEM, the individual particles can be distributed arbitrarily in the domain.
To facilitate the readability, the correspondence formulation of state-based Peridynamics is first introduced in Sect. 2. In addition, it is shown that the resulting equations can be derived directly from the weak form. The theoretical foundations for the development of an alternative formulation in the framework of Peridynamics are shown in Sect. 3. After listing the conditions on spatial discretization schemes in Sect. 5, the following section is devoted to the introduction of the Peridynamic-Galerkin-Method. Finally, the performance of the new meshfree scheme in comparison to the FEM is demonstrated by means of several demonstrative examples including an eigenmode analysis, the patch test, and cases with loadings leading to very large deformations.

State based peridynamic theory
To distinguish the presented alternative formulation from the standard approaches, the main equations in the framework of state-based Peridynamics are introduced first. Further information on the classical theory can be found, for example, in [42] or [41].
The differential equation for calculating the response of structures under mechanical loading is based on the balance of linear momentum. In Peridynamics, the equations are presented in the Lagrangian description with respect to the initial configuration. In the local form, the divergence of stress together with the density of volume forces leads to an accelerationü of the considered mass point ρ 0ü (X, t) = Div P(X) + ρ 0b . (1) The mass point is identified by its position vector in the initial configuration X. In the balance of momentum, ρ 0 denotes the initial density, P the 1st Piola Kirchhoff stress tensor, andb contains the gravitational acceleration. The initial boundary value problem requires the specification of the displacement vectorū on the Dirichlet boundary ∂ B D as well as the derivative of u on the Neumann boundary ∂ B N . For the latter, the stress vectorT is normally chosen In addition, the initial values for the displacement vector u(X, t = t 0 ) = u 0 and for the velocity vector v(X, t = t 0 ) = v 0 must be present for each mass point. The state based Peridynamics [42] replaces the divergence of the stress by the difference of force vector state fields All mass points within the spherical environment of a mass point X defines its neighborhood H X . The subscript vector in the volume differential indicates the coordinate to be integrated. The radius of the spherical neighborhood is called the horizon δ. A bond corresponds to the difference X − X between the position vector of the mass point in the neighborhood X and X. The sum of all bonds within the neighborhood denotes its family. The term state describes the mapping of bonds within the neighborhood and is denoted by an underscore. Different kinds of brackets are used, each type having a special meaning. In Peridynamics, angle brackets are usually employed to indicate a vector on which a state operates [42]. Square brackets contain differences of vectors. Arguments of functions are given in round brackets. The bond to be mapped is written inside an angle bracket. The force vector state field T X corresponds to the mapping that assigns a bond in the family of the mass point X a vector with unit [N/m 6 ]. The superscript vector indicates the family in which a quantity is evaluated. The second force vector state field T X , on the other hand, maps the same bond from the family of the mass point X .

Correspondence formulation
In [42] a description is presented which allows an application of the material models of the continuum theory in Peridynamics and is called correspondence formulation. This relation is based on the deformation vector state field Y, which maps a bond into the corresponding position vectors of the current configuration x Equivalence with the deformation gradient is given if both mappings lead to the same result In this case, the strain energy density of the peridynamic formulation V corresponds to the continuum mechanical equivalence W . A reduction operator R can be defined which approximates the deformation gradient F from the deformation vector state field [42] The operator * indicates the tensor product of two states and ω corresponds to the influence function. A precise definition and a list of the most important mathematical calculation rules for states can be found, for example, in [41]. To facilitate a distinction with the arguments of a function, differences are formulated in square brackets. The shape tensor K is determined from the tensor product of two reference position vector states X This product is complete and the inverse exists. The 1st Piola Kirchhoff stress tensor is calculated for hyperelastic material behavior from the derivative of the strain energy density W with respect to the deformation gradienṫ The equivalence (5) also holds for the time derivative. With the definition of Y from (4), a connection can be established between the force vector state field and the 1st Piola-Kirchhoff stress tensoṙ The operator • denotes the scalar product of two states. Considering (6) and the symmetry of the shape tensor K, this relation can be represented in an alternative form The comparison between (9) and (10) provides the relation between the 1st Piola Krichhoff stress tensor and the force vector state field Remark To show the analogy and the differences between the continuum mechanical and the peridynamic formulation, the usual abbreviations of continuum mechanics are preferred.
In [42] the peridynamic strain energy density is abbreviated by the letter W and the position vectors of the initial configuration by x and of the current configuration by y.

Equivalence with weak form
The peridynamic balance of linear momentum (3) can also be derived from the weak form of (1). This gives also some information about the imposition of Neumann boundary conditions. The weak form of the balance of linear momentum results from multiplying (1) by a test function η together with an integration over the volume of the body B ⊂ R d with dimension d. By applying the Gauss integral theorem, the order of the differential equation is reduced by one The Neumann boundary conditionT is part of the balance equation. The test function must be zero at the Dirichlet boundary, i.e. η = 0 on ∂ B u , sinceT is not known at these mass points. The second term on the left-hand side corresponds to the virtual internal work The transition to the peridynamic formulation results from the Taylor series expansion of the test function around the mass point X The gradient can now be determined using the method of least squares. Only the linear part in the Taylor series expansion is considered and the error in the integral over the whole neighborhood is minimized Therein, the residual in the integral is additionally weighted by the influence function ω. The operator we are looking for is determined by the derivative of the functional J with respect to the gradient of the test function The tensor K corresponds exactly to the shape tensor (7). Inserting (16) into the weak form and taking advantage of the symmetry of K results in an alternative representation of the virtual internal work The second factor in the integral corresponds to the mapping by the force vector state field (11), and the virtual internal work can also be formulated in terms of T Since T X X −X = 0 if X / ∈ H X , the integral can be extended over the whole domain. The indices can then be interchanged and the right-hand side can be reformulated Thus, the weak form of the balance of momentum (12) can also be approximated as the difference of force vector state fields Since the test function can be arbitrary and thus non-zero, the usual equation of Peridynamics (3) is obtained for inner mass points For mass points at the loaded Neumann boundary (T = 0) this relation does not hold. In this case, the last term on the right-hand side in (20) must be included.

Remark
The representation in (15) corresponds to the initial form for deriving peridynamic derivatives, as given, for example, in [32].

Relation to weighting functions
Smoothed Particle Hydrodynamics is based on the discretization by radial weighting functions W , see e.g. [47]. Also the equations within Peridynamics can be formulated by means of W using the definition This representation is also drawn in [16]. The superscript indicates the neighborhood in which the function is evaluated. This label is equivalent to the indication of the element in the FEM. From (21) a useful property of radial weighting functions follows Also the shape tensor (7) can be formulated in terms of W Taking advantage of (22), the shape tensor corresponds exactly to the negative value of the correction tensor L from [35], which ensures the 1st order reproducing condition in the derivatives Thus, a corrected weighting functionW can be defined The peridynamic derivative of the test function (16) can thus be specified in an alternative form With this interpretation, the virtual internal work from (17) can be reformulated accordingly Thus, an alternative form for the relation between the force vector state field and stress tensor can be drawn Using (21) and applying the concept of averaged gradients (26) the calculation of the deformation gradient can also be transformed

Averaged correspondence theory
In satisfying the essential conditions on spatial discretization schemes, care must also be taken to ensure that the stiffness matrix is always symmetric in a conservative system. This can be achieved in the framework of Peridynamics if the equations are formulated on the basis of an averaging. For this, the product of the gradient of the test function and the stress tensor is integrated over the neighborhood and then divided by the corresponding volume V H X [5] Compared to the classical correspondence formulation, the stress within the integral can now take different values. For a better identification, the neighborhood in which the stress is calculated is indicated by a superscript vector. Considering the derivative rule (26), the gradient can be formulated as a function of a difference. It should be noted that the gradient is determined with respect to the neighborhood of X. Thus, compared to (26), the first factor in the integral does not change Substituting (31) into (30) leads to a formulation of the virtual internal work G int as a function of differences From the comparison with (18), a new relation between the stress tensor and the force vector state field can be obtained, which now results from averaging The corresponding quantity from (19) with respect to the neighborhood of X is calculated in an analogous way For hyperelastic material behavior, the stress at the point Y is determined from the deformation gradient at the corresponding point . For its calculation the derivative rule (31) is used again

Discretized averaged Peridynamics
In the FEM, the body B ⊂ R d is divided into individual elements. The discretization concept of Peridynamics can be assigned to the one-point methods [47]. In this case the body is represented by a finite number of particles which can be identified by their respective position vector. The evaluation of all quantities is limited solely to these points. In the discrete form the integral is replaced by a sum. The balance of linear momentum in the state-based peridynamic framework derived from the weak form (20) is thereby evaluated at all particles n in the domain The right-hand side is only taken into account for particles at the loaded Neumann boundary. The abbreviation n s denotes the total number of boundary particles, N K the number of neighboring points of the particle X K ∈ R d , and V K ∈ R its corresponding volume. The number of particles in the corresponding neighborhood is usually obtained by a search algorithm, where the horizon δ defines the radius of a sphere around the corresponding point. More details can be found in Sect. 5.6. In a pure non-local formulation, the horizon is assumed to be constant, which offers advantages especially for the computation of cracks [36]. Alternatively, the radius can be defined as a function of the distance between individual particles or, if a subdivision of the domain into individual elements is present, the neighborhood can also be determined from the given connectivity. In the context of this work, only the second variant is considered in order to allow a comparison with the FEM. Nevertheless, with this formulation the radius can also be assumed to be constant. The discrete force vector state field is determined from (33). The weighting functionW multiplied by the volume can also be interpreted as a shape function. Thus, for the discrete derivative of (25) it follows To simplify the representation, the explicit specification of the coordinate is omitted except for the shape function, i.e., e.g., V (X K ) = V K . The superscript indicates the neighborhood of the corresponding particle. The formulation by means of shape functions allows a generalization, since the functions used need no longer be radial [5]. The discrete mapping of the force vector state field from (33) multiplied by the volume can hence be expressed alternatively Similarly, the mapping of the force vector state field at particle I is calculated using the information of its neighborhood The discrete form (36) can also be interpreted as a balance of forces that must be satisfied at each particle. The internal force due to the state of stress inside the material f s , the fraction from self-weight f g , and the inertia term f t must be equal to the force f n for particles at the loaded Neumann boundary and zero at all other particles The first term on the left-hand side corresponds to the difference of T from (38) and (39) integrated over the neigh- where N K I and N I K correspond to the test shape functions. The proportions of the self-weight and the inertia term are determined only from quantities at the corresponding particle Assuming hyperelastic material behavior, the 1st Piola Kirchhoff stress tensor is calculated from the derivative of the strain energy density W with respect to the deformation gradient The discrete deformation gradient follows from (35) where N K I now correspond to the approximate or trial shape function

Requirements on spatial discretization schemes
A prerequisite for an accurate solution of initial boundary value problems is the fulfillment of essential criteria. The solution of the differential equation must exist and be unique. In addition, the available approximates must contain the true solution. The verification of these requirements belongs to the field of functional analysis. These conditions are assumed to be satisfied. Furthermore, the discrete equation must fulfill several criteria, so that the true solution can be approximated as accurately as possible. Since only quasi-static cases are investigated, the requirements on the temporal integration procedure are neglected. In the following section, criteria on spatial discretization schemes are presented. A detailed treatment and comparison between different spatial discretization concepts can be found in [47]. The finite element method (FEM) fulfills all conditions, although many of the criteria presented here are not explicitly mentioned in the textbooks on the FEM.

Reproducing conditions
A central requirement on numerical solution schemes is to ensure convergence. The approximation must tend to the true solution at a finer resolution. For this, the method must be consistent and stable. Since consistency is difficult to prove for most numerical solution schemes, polynomial completeness or reproducibility of the shape function is one requirement formulated for the FEM. This condition must solely be met by the trial shape functions used to calculate the deformation gradient in (44). The derivation can be based on either the Weierstrass theorem or the Taylor series, as described in detail in [47]. An illustrative explanation of the meaning of this criterion can be found in the context of FEM in [21]. The degree of reproducibility which has to be met corresponds to the highest occurring derivative. As can be seen from (30) and (35), at most the first derivative occurs in the peridynamic formulation. The 0th order reproducing condition corresponds to the socalled partition of unity. This criterion can also be represented in the form of the derivative. For the shape functions, which are used in the neighborhood of X K , the following must be true Compliance with this condition ensures that a constant quantity is properly approximated and its derivative is zero. Due to the difference-based formulation, the second requirement is always satisfied, see (44). The 1st order reproducing con-dition is ensured if the approximation of the position vector corresponds exactly to its value. The condition on the derivative results in the requirement that the sum corresponds to the unit tensor For differences, the equivalent form reads If the equations in (45) are satisfied, then (47) equals (46).

Integration constraint
In [25] a condition on the test shape function is formulated.
In [11] this criterion was further specified in the context of meshfree Galerkin methods and called the integration constraint. A generalization can be found in [12]. This condition addresses the Gaussian integral theorem, which must also be satisfied in the discrete case. Since the equations are derived from the weak form, Peridynamics must also fulfill this requirement. This theorem reads In the Peridynamics framework, the second term on the righthand side is formulated in terms of a difference, see (20) For the classical peridynamic formulation, the relation between force vector state fields and stress tensor is given in (11). When using the averaged approach, the individual contributions from (33) and (34) must be used. Different orders also exist for the integration constraint. If only linear shape functions are used, it is sufficient if the integration constraint is satisfied for a constant stress in the considered area. For the averaged correspondence formulation the requirement in (49) reduces to The discrete form of this condition results in a criterion that must be satisfied at each particle. The left-hand side must vanish for particles in the interior and be equal to a resulting normal vectorN at the boundary. This quantity is calculated from the normal vector N multiplied by the area fraction of the particle, i.e.

Discrete conservation properties
Balance equations must also be valid in the discrete case. Only external forces lead to a change of momentum. If these are not present, the total part from f t in the body must be zero [3]. Since there is no loading from self-weight and Neumann boundary conditions, it follows from (40) that the internal forces from the state of stress must vanish in sum. Similarly, it follows from the discrete conservation of angular momentum that the cross product of internal force and position vector must be zero Due to the peridynamic formulation based on differences, the conservation of linear momentum is always satisfied. Since the force vector state fields outside the neighborhood are zero, the sum can be extended over the whole domain. Thus, the indices can be interchanged As pointed out in [42], the requirement for conservation of angular momentum from (52) can also be represented in terms of a difference of position vectors n K =1 In the averaged approach, substituting (38) into (54) and taking advantage of the 3rd order permutation tensor E, the condition can be further transformed by considering (44) n K =1 Due to the symmetry of PF T , the proposed new peridynamic formulation per se satisfies the angular momentum conservation requirement [5]. An exact specification for the tensor E can be found, e.g., in [20].

Configurational consistency
Another important property that must be satisfied during the simulation is the preservation of pull-back or push-forward operations [47]. For instance the current coordinates obtained from the solution of the initial boundary value problem must also be determined from the update using the deformation gradient. Thus, for each bond in the corresponding neighborhood of the particle X K it must hold The left-hand side corresponds to the calculated coordinates from the solution of (36). This condition often serves as a starting point for the development of stabilizations aimed at correcting the negative effects due to under-integration, see e.g. [27], [39], or [10]. However, the above criterion must be satisfied separately, as shown, for example, by simple studies using the Optimal Transportation Meshfree method [47].

Stability
In spatial discretization schemes, either a tensile or a rank instability can occur. The first case can be circumvented if the equations are formulated with respect to the coordinates of the initial configuration [4]. This is given in the approach presented here. Thus, no tensile instability occurs.
In one-point methods, like Peridynamics, the evaluation takes place only at the particles. Thus, the number of evaluation points is fixed. In the classical correspondence formulation of Peridynamics, effects occur that can be observed in the case of under-integration within the FEM [45], [6]. For a detailed investigation of this unphysical behavior spectral analysis is well suited. From the condition the eigenvalues λ and the eigenmodes N of the stiffness matrix K can be determined. The eigenvalues belonging to the rigid body modes are identically zero. An inaccurate evaluation or under-integration is present if further eigenvalues are zero or very small. The eigenvectors indicate the deformation mode of the associated eigenvalue.

Search algorithm
In meshfree methods, a search algorithm determines the neighborhood of each point. In Peridynamics, the search region is usually spherical and the radius is called the horizon δ. If the distance between the particles X J and X K is smaller than the horizon, then X J is a part of the neighborhood H K of the point X K The search algorithm must also meet certain requirements. For example, it must be possible to calculate the inverse of the shape tensor. For this, the point distribution in the neighborhood should be homogeneous and isotropic on statistical average [29]. In addition, there must be no pathological distributions, i.e., in the 3-dimensional continuum there must be 4 points in the neighborhood spanning a tetrahedron. Furthermore, the search algorithm must ensure that the requirements of Sect. 5.3 are not violated. Equations (53) and (55) are automatically satisfied if pairs of forces are possible. If the particle J is part of the neighborhood of the particle K , then K should also be part of the neighborhood of J . If this is not given, unphysical forces can occur, which are called ghost forces in [37]. If the horizons in all neighborhoods are equal, this requirement is always satisfied.
If a subdivision of the domain into individual finite elements is available, the connectivity can also be used to define the neighborhood of a particle. This automatically results in the required pairs of forces.
For the description of solid deformations due to external loadings, the equations are usually expressed in the form of a total Lagrangian description. In this case, the computational effort is limited to the determination of the neighborhoods before the simulation. The connectivity is assumed to be constant later on.

Kronecker-ı property
In standard discretization approaches, such as in the FEM, the solution is approximated by the product of base or shape functions and unknown coefficients a If the shape functions possess the Kronecker-δ property, the coefficients correspond to the values of the approximate solution at the corresponding particle This property allows the direct imposition of Dirichlet boundary conditions in computational solution schemes.

The Peridynamic-Galerkin-Method
Starting from the averaged correspondence formulation from Sect. 3, a meshfree formulation is presented that satisfies all the requirements for a spatial discretization scheme. Therefore, the right choice of shape functions is of central importance. Due to the different requirements on test and trial functions a Petrov-Galerkin approach seems to be preferable. An application of this concept can be found in the Peridynamics framework in [6], [7] and [8]. On the other hand, a conservative system must lead to a symmetric stiffness matrix, see e.g. [49]. This can be given if the test and the trial shape functions are identical. In order to fulfill the requirement for conservative systems, the same shape functions are used for η and u in this approach. The presented averaged correspondence formulation satisfies per se the requirements for the discrete conservation equations. In the examples, the particle distribution can be considered statistically homogeneous and the conditions from the search algorithm are satisfied. Hence, no ghost forces occur. Since the formulation is based on a total Lagrangian description, no tensile instability emerges. As shown in Sect. 7, no negative effects due to under-integration can be detected due to the formulation within the framework of the averaged correspondence theory. An evaluation of the quantities at the individual bonds, as described in this formulation, counteracts these effects, see also [6].
Thus, the requirements that remain to be satisfied are the reproducing conditions, the integration constraint, and the configurational consistency. In addition, due to a modification, the shape functions possess the Kronecker-δ property.

Reproducing conditions and shape functions
The fulfillment of the reproducing conditions up to the highest occurring derivative represents an advantageous criterion to enable a convergent behavior of the computational solution scheme. Weighted or moving least square functions have the big advantage that the reproducing conditions up to the desired order can be guaranteed even for an arbitrary number of particles and distribution in the neighborhood. In [29] an approach can be found which is also suitable for Peridynamics, since it is based on differences and an integral formulation. According to the Weierstrass theorem, any function in a closed domain can be approximated with sufficient accuracy by a polynomial The argument of the base function vector p contains the variables of the polynomial. In this case, the argument is based on the difference between X and X. In the approach of [29], the unknown coefficients d are calculated from the minimization of a functional J . The connection of d to the coefficients a from (59) follows from the Kronecker-δ property. The error in (60) should be as small as possible in the weighted average over the integral of the neighborhood The unknown coefficients d are calculated from the requirement that the derivative of the functional must be equal to zero The moment tensor M corresponds to a complete dyadic product. Thus, the inverse exists provided the particle distribution is not pathological Substituting (63) into (61) the displacement vector can be represented as a function of the coefficients a In the discrete approach, weighted least-square shape functions can be defined Moving least-square shape functions result from the transition This transition must also be considered in the moment tensor In order to satisfy the 1st order reproducing conditions, a linear approach is sufficient in the weighted or moving least square shape functions, i.e. for the 3-dimensional case the base function vector reads

Kronecker-ı property
The discrete moving least square shape functions generally do not satisfy the Kronecker-δ property. To ensure interpolation [26], an influence function can be chosen as Since the influence function in a one-point method must also be evaluated at the central particle, the denominator can be zero. To avoid infinity, an augmentation can be used. A disadvantage of this approach is a poor condition of the mass matrix. An alternative can be found in [8] or [5]. In the base function vector, the constant fraction is neglected, i.e., p (X I − X J ) = X I − X J . To compute the derivative of the shape function, first the weighted least-square formulation in (65) is differentiated with resect to X . Afterwards, the transition X → X is made leading in the discrete case to The influence function corresponds to Eq. (70) and also in the moment tensor the reduced base function vector is used The shape function derivative associated with the particle X J is computed so that the 0th order reproducing condition (45) is fulfilled. Using this modification the singularity can be bypassed and the requirements from Sect. 5.1 remain satisfied.

Configurational consistency
The correspondence formulation of Peridynamics does not guarantee the configurational consistency. In [10], a scheme is presented that can be used to satisfy the condition from (56). This approach can also be interpreted as a correction of the shape function derivative to compute the deformation gradient in (44) The modification is solely based on bonds On the right-hand side, the derivative of the shape function from (71) is employed.

Integration constraint
The fulfillment of the condition from Sect. 5.2 is not automatically met when using weighted or moving least-square shape functions. Therefore, a suitable correction scheme is needed. An approach in the framework of Peridynamics can be found in [8]. The idea of this correction goes back to [9], which is developed for Smoothed Particle Hydrodynamics. In this case, an additional term is added to the derivative of the shape function given in (74) The correction functions K ,d I are predefined. The number of unknown coefficients α d per node corresponds to the dimension n dim of the problem under consideration. The direction is defined by the Cartesian basis vectors e d . The individual values for α d are calculated from the satisfaction of the integration constraint from (51). The individual contributions of each particle can also be summarized in the form of a vector The symbol A corresponds to the assembly operator and assigns the correct location in the vector to each contribution. Furthermore, it must be ensured that the corrected shape functions do not violate the reproducing conditions and configurational consistency. For this, the additional part must satisfy (45) and (46) in the derivatives For a better representation, the individual contributions K ,d I are combined into a vector K I . This correction must not violate the configurational consistency. Thus, from (56) together with the definition of the discrete deformation gradient (73), the following must hold additionally A modification of the approach given in [34] allows all three requirements to be satisfied The shape function in (79) corresponds to the formulation given in (66). Equation (76) can also be represented as a linear system of equations with the unknown vector α. The solution requires proper boundary conditions. For this purpose, the part of α perpendicular to the boundary is set to zero [8]. These modifications yield a shape functionÑ that satisfies all the conditions on discretization schemes presented in Sect. 4.

Imposition of boundary conditions
Due to the Kronecker-δ property and the fulfillment of the reproducing conditions of the shape functions, Dirichlet boundary conditions can be imposed directly analogous to the FEM. Neumann boundary conditions require an integral over the surface, see (20). The surface is normally not known in meshfree methods. As shown in [46] for the FEM, stresses at the boundary can be directly transformed into volume loads using the integration constraint, provided that the condition (51) is satisfied. Therefore, no knowledge of the surface is required. The magnitude of the resulting normal vector corresponds exactly to the area fraction of the particle. For instance a pressure boundary condition can be imposed using (51) with the property This kind of application differs from the classical treatment in the framework of Peridynamics. These approaches mostly use ghost particles to impose Dirichlet boundary conditions. The number and the position of these particles are chosen such that the outer particles have a complete neighborhood [31], [28]. The Dirichlet boundary conditions are defined at the ghost particles. In [40] a linear distribution of these values over the position of the ghost particles is proposed. An alternative approach uses a correction of the shape function at the outer layer to apply the displacements directly at the boundaries [50]. So-called general meshfree approximation functions [51] are introduced, which possess a weak Kronecker-δ property. Neumann boundary conditions are converted into corresponding volume forces [32]. The major advantage, of the new peridynamic formulation presented here, is the application of boundary conditions analogous to the FEM, without the need for additional ghost particles and the imposition over a certain layer.
Using the modifications given in this section a meshfree method results that satisfies the essential conditions on spatial discretization schemes and simplifies the imposition of boundary conditions.

Examples
The Finite-Element-Method provides a good approximation in many cases. When solving the differential equation in a Lagrangian description, strong deformation limits the use of this spatial discretization scheme. The investigations presented in this section aims to show that the Peridynamic- Galerkin-Method (PGM) results in a very similar approximation as the FEM, but also allows for much stronger deformations. In addition, the regularity on particle distributions is drastically reduced. The selection of examples addresses different characteristics in the framework of spatial discretization schemes. First, the rank instability is investigated using an eigenmode analysis. Subsequently, the patch test is conducted, which gives an indication of the consistency of the method. A convergence study is presented in the context of a demanding test. Very large deformations usually lead to major challenges in mesh-based methods. Especially a poor initial mesh can lead to erroneous solutions or to the termination of the simulation. For this reason, three different applications are investigated in the context of finite deformations. The specimens are subjected to large pressures, a large torsional as well as large shear and bending loads. Especially, the case of poor initial mesh is evaluated.
In all examples, linear elements are used in the FEM analysis and linear base functions vectors in the PGM. To allow comparability, the nodes correspond to the particles in the peridynamic analysis. In the FEM, resulting normal vectors can be easily assigned to each node at the boundary, see [46].
Peridynamics is based on a non-local formulation where the horizon δ is usually assumed to be constant [31]. This enables convergent behavior to be obtained even in the calculations of fracture [36]. To allow comparability between the new peridynamic formulation and the FEM within a convergence analysis, in this section the horizon is adjusted upon refinement so that the number of particles in the neighborhood remains constant. Sometimes, the connectivity from the FEM is used to derive the neighborhood. However, to achieve convergent behavior in crack calculations, the horizon can also be assumed to be a material parameter having a constant value.

Eigen-mode analysis
As stated in Sect. 5.5, rank instability must be avoided in order to obtain physically meaningful results. Spectral analysis provides a suitable approach to determine if such a case exists. In the Finite-Element-Method, the stiffness matrix of an element is decomposed into eigenvalues and eigenvectors, see e.g. [49]. Due to the overlap of the individual neighborhoods, an investigation on the local level has only little significance for meshfree methods. Thus, negative effects due to under-integration may occur locally, but disappear due to the overlap at the global level. For this reason, the eigenvalues and eigenvectors of the entire system are considered.
For a better illustration, a 2-dimensional block with a length of 20 mm and a height of 10 mm is studied. The material behavior is assumed to behave linear elastic with a Young's modulus of 100 N/mm 2 and a Poisson's ratio of ν = 0.4 In the linear case, P corresponds to the Cauchy stress tensor.
To keep the representation consistent, the 1st Piola-Kirchhoff stress tensor is always used in this work. The formulas for converting E and ν to the Lame parameters λ and μ in the case of plane strain can be found, for example, in [49].
For the FEM analysis, the block is divided into 20 × 10 linear quadrilateral elements (Q1). The neighborhood is determined by the FE connectivity or is spherically spanned by a radius of δ = 3.01 , where corresponds to the particle spacing in each spatial direction.
Next to three rigid body modes, which are present to all formulations, the smallest 200 eigenvalues of the global tangent stiffness matrices are depicted in Fig. 1. The individual values of the Peridynamic-Galerkin-Method (PGM) agree well with the FEM Q1 eigenvalues in the range of low stiffness. Note that the enlargement of the family size to 3.01 times the nodal spacing in the PGM only shows a very small effect on its eigenvalues. The classical peridynamical formulation (PDC) is known to suffer from spurious oscillations. This is reflected in a large number of eigenmodes with low associated energy. In Table 1, the number of the lowest 11 eigenvalues which not belong to the rigid body modes are presented. The last two lines in Fig. 2 show the first three eigenmodes for the case of the PGM and the FEM, which are assigned to the corresponding eigenvalues from Table 1.
Here a good agreement can be found. The first row contains the modes of the PDC, which belong to the first two, as well as to the eleventh eigenvalue from Table 1. It is remarkable that the 11th eigenvalue of the PDC with the corresponding mode agrees with the 3rd eigenvalue and mode of the FEM  or the PGM, respectively. The shape of the modes of the first 10 eigenvalues of the PDC show a typical behavior which can be observed especially in the case of under-integration [6] or [49].

Patch test
To avoid a direct mathematical verification of consistency in the FEM the patch test introduced by [2] can be used instead.
A confirmation of equivalence to consistency is given in [43]. An extension to also study the stability of a formulation can be found in [44]. In addition, different versions exist. In this study the test given in [12] is used. The material behavior is assumed to be linear elastic. The true solution of the differential equation is given at the boundary. The results at the interior nodes or particles from the simulation must match the true solution in order to pass the patch test. For the linear . 3 a Deformed configuration of the bloc and boundary conditions. b Shear stress distribution σ xy at the particles using PDC. c Shear stress distribution σ xy at the particles using PGM. d Legend of σ xy case, the following solution is considered The stress is calculated according to Eq. (82) with a Young's modulus and a Poison's ratio of E = 100 kN/m 2 and ν = 0.3. The stress in case of plane strain is given by  (Fig. 3), which is derived from the stress tensor (84). At the other sides, Dirichlet boundary conditions are imposed according to the values from (83). Fig. 3 depicts the resulting shear stress field, which should be constant. 1345 irregularly distributed particles are used in this test. In order to be able to achieve a comparison, the neighborhood of the PGM corresponds to the connectivity of the FEM. The distance between the particles and hence also the radius is approx. 0.03 m in this case. In contrast to the classical peridynamic formulation PGM shows a constant stress field. The error of the approximate is evaluated using the Lebesque quadratic (L 2 ) norm. For the derivative the first Hilbert H 1 -seminorm is employed The first value in the difference corresponds to the true solution and the second term u h denotes the approximate. For more information on norms in the context of numerical solution procedures, see, for example, [1].
As can be observed in Table 2, the PGM satisfies the patch test compared to the classical peridynamic correspondence formulation. The PDC only fulfills this test, if the particles are regularly arranged and ghost particles are used to impose Dirichlet boundary conditions. Only in this case, the PDC satisfies the integration constraint. Effects due to underintegration are not so dominant in the patch test.

Convergence in a manufactured 2-D problem
To check the convergence behavior, the test from [18] is considered. The true solution in this case has a sinusoidal shape The material behavior is assumed to be elastic according to (82) with the parameters E = 100 kN/m 2 and ν = 0.3. The body is additionally loaded by volume forces For a set of refinements with N + 1 particles per meter in each dimension, the norms from (85) are evaluated. Fig. 4 depicts the undeformed configuration for N = 4, N = 16 and N = 64. The neighborhood is determined by the FEM connectivity. For this purpose, quadratic quadrilateral  In Fig. 6, the error norms with the averaged convergence rate are plotted. The PGM shows a convergence behavior with rate 2 in the L 2 norm. In the case of the H 1 semi-norm, a better rate can be obtained compared to the FEM using linear polynomial shape functions. Using the classical peridynamic correspondence formulation, on the other hand, no convergent behavior exist. During loading, the boundary gets curved which makes this test more challenging for imposing Dirichlet conditions.

Punch test
If the elements distort too much, mesh-based methods require remeshing, otherwise usually no or only a poor approximation can be achieved. Meshfree methods on the other hand offer a larger flexibility because the connectivity does not have to be sorted and the point distribution can be arbitrary. To test the new peridynamic method, a 2-dimensional punch test is conducted. As shown in Fig. 7, a block is compressed on the left side of the top surface by a displacement ofū = −1 mm.
The PGM results are compared with three different FEM formulations. In addition to linear quadrilaterals (Q1), two enhanced assumed strain elements based on a Taylor expansion of the shape function derivatives are also used which provide good results for more extreme loading cases. More information regarding the corresponding H1TS and the TSCG12 formulation can be found in [23] and [24], respectively. For the study with these approaches, the domain is divided into 3-dimensional elements. In order to represent the plane strain state, the displacements in the third direction are fixed. The domain is divided into 64 × 32 elements in x-and y-direction. In the 3-dimensional case, only one element is used in the z-direction. The depth is adjusted so that a cube-shaped geometry can be obtained. The material behavior is now non-linear elastic. The stress is calculated from a neo-Hookian approach. The strain energy density W is thereby decomposed into an isochoric and volumetric part The 1st Piola-Kirchhoff stress tensor is determined from the derivative with respect to the deformation gradient The Lame parameters are calculated from the Young's modulus E = 333 MPa and the Poisson's ratio ν = 0.2. The maximum displacement in the x-direction is plotted against the minimum displacement in the y-direction for different FEM formulations and the Peridynamic-Galerkin-Method   Fig. 7 can be found in Table 3.
The final deformed configurations of the PGM simulation is shown in Fig. 8.

Torsion of a square hyperelastic prism
Robustness to highly irregular discretizations is a major advantage of meshfree methods. To demonstrate the independence of the particle distribution, the torsion test according to Fig. 9 is examined. A squared prism of dimensions (length × width × height) 5 mm × 5 m × 125 m is clamped at the top as well as the bottom and rotated at the top end. The material behaves non-linearly elastic according to (89) with a Young's modulus of E = 12M Pa and a Poison's ratio of ν = 0.48.
The block (length × width × height) is divided into 10 × 10 × 125 cubes. The inner nodal coordinates are randomly perturbated in each spatial direction in an interval of (−d, d) times the initial nodal spacing, where d represents a dimensionless quantity whose magnitude can be at most 1. In total, 20 different discretizations are created using randomly distributed values for d.. In the case of the PGM, the nodes correspond to particles. All nodes of the adjacent elements define the neighborhood of the corresponding particle.
The rotation of the top end causes the prism to get twisted. At a critical rotation, an additional bending occurs as can be observed at rubber bands. All calculations are conducted up to this critical rotation angle. The calculation using the PGM shows a good agreement with experimental results (Fig. 9). In Fig. 10, the maximum possible angle of rotation is shown for different perturbations. For the calculation by FEM, the enhanced assumed strain formulation (TSCG12) from [24] leads to minimally improved results compared to the standard trilinear brick element (H1). However, the PGM again outperforms these formulations. As can be seen from Fig. 10, the failure angle of the Perdynamic Galerkin method is independent of the degree of perturbation. The H1TS element which shows improved results in the punch test allows only a much smaller angle of rotation. For all FEM formulations, at a certain perturbation a calculation beyond a small rotation angle is no longer possible if d exceeds a certain value. The outcomes also show good qualitative agreement with measured results. In the experiment, the critical torsion angle was about 2340 degrees (Fig. 9). The Young's modulus in the simulation were taken from the data sheet of the printed material. The Poison's ratio was estimated. However, the material can be characterized as weakly compressible. It was also assumed that the stress-strain behavior can be represented by a neo-Hooke's model.

Cantilever test
In the last test, the computational effort between the FEM and the PGM is compared. For this purpose, at the free end of a cantilever a surface load is applied. The dimensions, boundary conditions and material parameters can be found in Fig. 11.
The material behavior is assumed to be non-linear elastic. The stress is calculated according to (89). Due to the geometry, a bending and shear stress dominates in the body. For the FEM calculation, the cantilever is divided into linear tetrahedra (FEM O1) and as an additional variant again into the enhanced assumed strain elements from [24] (FEM TSCG12). The nodes of the adjacent tetrahedral elements define the neighborhood of the corresponding particle. The black horizontal line in Fig. 12 denote the overkill FE result in the convergence analysis. As expected, the solution based on the FEM with linear tetrahedral elements behaves excessively stiff whereas the TSCG12 formulation shows an improved behavior. After a fluctuation at a coarse resolution, the PGM reaches the asymptotic value even before.
This figure also shows the computational effort between the FEM and the PGM. Both algorithms are incorporated into the same program structure which allows a fair comparison. While in the range of a few degrees of freedom the Comparison of computation time between the PGM and the FEM assembly of the system plays the main role, for larger numbers of particles the solution of the system of linear equations dominates. Thus, the increase in the computational effort for smaller systems is caused by an increased number of evaluations of the constitutive equations and for larger systems by the less sparse tangent matrix resulting from overlapping families. For the finest resolution with 1.2 million degrees of freedom, the total computational cost of the PGM simulation is about 56% larger than for the corresponding FEM O1 simulation. The computational effort of the enhanced assumed strain element formulation is in the range of the PGM. On the other hand, the accurate solution can be obtained already at a coarser resolution of the domain by the PGM in comparison to FEM O1. Using linear tetrahedra a finer subdivision is required (Fig. 13).

Summary
The Finite-Element-Method offers a good approximation of the true solution for many engineering applications. In the field of solid mechanics, the differential equation is usually solved on the basis of a Lagrangian description. An important prerequisite is the subdivision of the domain into elements that are as regular as possible. However, complex structures often needs time-consuming manual post-processing to generate an acceptable regularity of the mesh. Very large deformations can also lead to strong distortions of the elements. In this case, remeshing of the domain is needed.
Meshfree methods promise the possibility of a good approximation of the solution with an almost arbitrary distribution of points in the domain. In this case, meshing is not necessary and even for very large deformations the mechanical differential equation can be solved in a Lagrangian description. On the other hand, often a good approximation of the true solution cannot be obtained. The reasons for this lie in the violation of important conditions on spatial discretization schemes.
The presented Peridynamic-Galerkin-Method (PGM) fulfills all these conditions. Moreover, the negative side effects due to a rank deficiency, which is usually a critical point in meshfree methods, are avoided. The results of the patch test as well as the manufactured test show that this method can approximate the true solution of the differential equation very well. Due to the meshfree character, loadings can be calculated that lead to very large deformations compared to the FEM. In addition, the requirements for the distribution of points in the domain are significantly reduced using the PGM.
In the context of this work, cases with additional constraints, such as incompressibility, were not investigated. In addition, the application is limited to linear basis functions. To ensure comparison with the FEM, the horizon was not identified as a material parameter, but was considered to depend on the particle spacing. Nevertheless, the equations can also be applied to the case of a constant horizon.