The computational framework for continuum-kinematics-inspired peridynamics

Peridynamics (PD) is a non-local continuum formulation. The original version of PD was restricted to bond-based interactions. Bond-based PD is geometrically exact and its kinematics are similar to classical continuum mechanics (CCM). However, it cannot capture the Poisson effect correctly. This shortcoming was addressed via state-based PD, but the kinematics are not accurately preserved. Continuum-kinematics-inspired peridynamics (CPD) provides a geometrically exact framework whose underlying kinematics coincide with that of CCM and captures the Poisson effect correctly. In CPD, one distinguishes between one-, two- and three-neighbour interactions. One-neighbour interactions are equivalent to the bond-based interactions of the original PD formalism. However, two- and three-neighbour interactions are fundamentally different from state-based interactions as the basic elements of continuum kinematics are preserved precisely. The objective of this contribution is to elaborate on computational aspects of CPD and present detailed derivations that are essential for its implementation. Key features of the resulting computational CPD are elucidated via a series of numerical examples. These include three-dimensional problems at large deformations. The proposed strategy is robust and the quadratic rate of convergence associated with the Newton--Raphson scheme is observed.


Introduction
Peridynamics is an alternative approach to formulate non-local continuum mechanics [1]; its roots can be traced back to the pioneering works of Piola [2,3]. However, it is fundamentally different from established non-local elasticity [see 4, 5, among others] as the concepts of stress and strain are absent. As a non-local theory, the behaviour of each material point is influenced by interactions with other material points in their finite vicinity. In contrast to CCM, the governing equations of PD are integro-differential equations appropriate for problems involving discontinuities Before defining the kinematic measures of CPD, we recall the three local kinematic measures of relative deformation in CCM, namely the deformation gradient F := Grad y, its cofactor K := Cof F and its determinant J := Det F.
In the spirit of these local measures, we introduce three non-local kinematic measures of relative deformation, namely ||| associated with CPD. The first relative deformation measure of CPD is ξ | . It mimics the linear map F from the infinitesimal line element dX in the material configuration to its spatial counterpart dx. In view of our proposed CPD formalism, the relative deformation measure is the main ingredient to describe one-neighbour interactions. The second relative deformation measure a | / || is reminiscent of the linear map K from the infinitesimal vectorial area element dA in the material configuration to its spatial counterpart da. This is essentially Nanson's formula from conventional continuum kinematics. In our proposed framework, the relative area measure is the main ingredient to describe two-neighbour interactions. The third relative deformation measure v | / || / ||| mimics the linear map J from the infinitesimal volume element dV in the material configuration to its spatial counterpart dv.
The relative volume measure is the main ingredient to describe three-neighbour interactions.

Governing equations
where b ext 0 denotes the external force density per volume in the material configuration, with units N/m 3 , and t ext 0 is the external traction on the boundary in the material configuration, with units N/m 2 . This format of the external loading is a particular sub-case of a more general case accounting for higher-gradient and non-local continua as detailed in [72,73] among others. The universal form of the quasi-static linear momentum balance (5) can alternatively be expressed in terms of volume integrals as The internal body force density in the material configuration b int 0 in CCM is the material divergence of the Piola stress P. In CPD however, b int 0 takes an integral form over the horizon. That is with p | the force density per volume squared, with units N/m 6 . Inserting the internal body force density (7) into the quasi-static linear momentum balance (6) yields, after localization, the linear momentum balance of CPD To derive the angular momentum balance, we start from the global form of the quasi-static moment balance which, after some mathematical steps and using the linear momentum balance (5), upon localization reduces to the angular momentum balance of CPD Table 1 gathers the key governing equations for both CCM and CPD for the case of quasi-statics. Table 1: Governing quasi-static equations of classical continuum mechanics (CCM) and continuum-kinematics-inspired peridynamics (CPD). The Piola stress is denoted by P and F · P t is symmetric due to angular momentum balance expressed using the third-order permutation tensor ε.
linear momentum balance angular momentum balance

Constitutive laws
Constitutive laws bridge the gap between the kinematics described in Section 2.1 and the kinetics in Section 2.2.
The constitutive laws of CPD, as in CCM, must be thermodynamically consistent and are thus derived via a Coleman-Noll-like procedure. Let Ψ denote the point-wise stored energy density per volume in the material configuration. The The interaction potentials (21) are generic examples of suitable one-, two-and three-neighbour interactions.
Finally, we propose specific examples of stored energy densities ψ | 1,2,3 that are both thermodynamically consistent and satisfy the angular momentum balance. These energy densities are employed in the numerical examples in Section 4 and their key characteristics are illustrated and discussed. For that it proves convenient to define the counterparts of l , a and v in the material configuration, denoted by L , A and V , respectively. That is An example of the stored energy density per volume squared in the material configuration for one-neighbour interactions ψ | 1 = ψ 1 (l ; L ) and in accordance with the original bond-based model of Silling [1] reads where C 1 is the one-neighbour elastic coefficient and can be viewed as the resistance against the change of length between a point and its neighbours, reminiscent of the elastic modulus in CCM. Note that the parameter S 1 is precisely the stretch of the bond from a point to its first neighbour. Motivated by the format of the stored energy density for oneneighbour interactions (23), we propose the stored energy density per volume cubed for two-neighbour interactions as with C 2 the two-neighbour elastic coefficient which can be interpreted as the resistance against the change of the area of the triangle formed by a point and a pair of its neighbours, analogous to Poisson-like effects in CCM for two-dimensional manifolds. The parameter S 2 is essentially the area stretch of this triangle. Similar to the stored energy densities for one-and two-neighbour interactions, we propose the stored energy density per volume to the fourth power for three-neighbour interactions as 15 and where C 3 is the three-neighbour elastic coefficient, which can be interpreted as the resistance against the change of the volume of the tetrahedron formed by each point and its triplet of neighbours, analogous to volumetric Poisson-like effects of CCM. The stored energy densities (23)- (25) are introduced in this fashion since the energy density (23) is identical to the common format used in bond-based PD [74].
Alternative formats for energy densities could be proposed and their consequences investigated. However, the main objective of this manuscript is to provide details of the computational implementation which remains largely independent of the specific format of the stored energy density. Table 2 summarises the discussion of constitutive laws and collects the fundamental relations and definitions of CPD together with generic and specific examples of the stored energy densities that sufficiently satisfy the angular momentum balance.
specific examples of interaction energy densities

Computational implementation
The computational implementation of CPD comprises three major steps. We begin by replacing the integral equations in Section 3.1 with quadrature relations using appropriate weighting coefficients. Next, in Section 3.2, a discretised form of the linear momentum balance (8) is derived. The discretised balance is a non-linear system of coupled equations of the form R = O that is solved using a Newton-Raphson scheme. To do so, we compute the tangent matrix K defined as the linearisation of the residual vector R. Finally, the force densities and their derivatives contributing to R and K, respectively, are calculated from the constitutive laws associated with the stored energy density given in Section 3.3.
Remark 1. The proposed computational framework corresponds to the three-dimensional setting. Nevertheless, both plane-strain and plane-stress assumptions can be recovered via applying appropriate boundary conditions on a 3D domain. Note that neither "stress" nor "strain" is present in the peridynamic formulation; they can only be computed through post-processing. Therefore, the notions of "plane strain" or "plane stress" become naturally less relevant since they are defined from a local view of continuum mechanics. Furthermore, it is straightforward to infer a fully two-dimensional counterpart from the provided discussion. Our formulation in 2D, however, corresponds to a purely two-dimensional case wherein both deformations and forces are absent in the third direction. This compares to the surface elasticity theory of Gurtin and Murdoch [75][76][77][78]. Obviously, three-neighbour interactions do not contribute in our 2D formulation. Both three-and two-dimensional numerical examples are provided in Section 4.

From continuous to discretise form
In general. The governing equations of CPD are expressed in integral form. Unlike CCM, even point-wise equations in CPD include integrals over the horizon. The points P a at which we evaluate the balance of momentum (8) are precisely the collocation points. At each collocation point P a , we use quadrature rules to evaluate the integrals over the horizon by employing quadrature points. In this contribution, the collocation points coincide identically with the quadrature points, henceforth, we refer to them collectively as grid points or simply points. This assumption is made for the sake of simplicity; alternatives will be investigated in a separate contribution.
To proceed, we first discretise the domain by a set of grid points that serve the double-purpose of being collocation and quadrature points. Every grid point in the present description represents continuum point as opposed to physical particles. Furthermore, each grid point defines a neighbourhood H 0 ⊂ B 0 . Each grid point P a is identified by its coordinates, X a and x a , in the material and spatial configurations, respectively. Therefore, the integral of an arbitrary quantity {•} over the domain B 0 is approximated by where V a is the volume of the support domain of the grid point P a and #P is the total number of grid points. The volume V a can be computed according to the discretisation strategy employed. For instance, if the points are chosen based on a Voronoi tessellation, the volume of each Voronoi cell can be assigned to the point P a at its centre, see [79].
It is also relatively straightforward to discretise the domain using the common discretisation tools of the finite element method and then associated each finite element with the point P a in its barycentre with V a equal to the volume of the respective element. Alternatively, for simple geometries, one can discretise the domain using a uniform grid for which V a = α 3 with the grid spacing and α a constant dimensionless correction factor accounting for the size of the support. If the grid spacing is non-uniform, the parameter can be replaced with a suitable average grid spacing avg . The summation on the right-hand side of Eq. (26) should correctly represent the volume of the domain itself in the sense that Next, we discretise the various multiple integrals that appear in the form Let #N denote the number of points within the horizon of the point P a . The effective volume of the neighbouring point i contributing to one-neighbour interactions at the point P a is denoted as V 1 , assuming that all neighbours equally contribute. The effective volumes V 1 can be defined by where V H denotes the volume of the neighbourhood of the collocation (continuum) point P a . For example if the point P a is completely inside the bulk and entirely surrounded (i.e. no part of its horizon extends outside of B 0 ) then V H = 4/3 π δ 3 0 . Otherwise, V H is modified by a dimensionless correction factor β < 1, as V H = β 4/3 π δ 3 0 where β accounts for the truncated support. In the definition of the effective volume for one-neighbour interactions (29), the total number of contributing neighbours within its horizon is denoted as #N 1 . A neighbour {i} of the point P a is identified as a contributing neighbour if the distance between the pair {a, i} is less than or equal to the horizon size of δ 0 . Clearly, all neighbours count as contributing neighbours for one-neighbour interactions and thus, #N 1 = #N. As will be made clear, this property does not necessarily hold for two-neighbour and three-neighbour interactions. The integral (28) 1 can be formally written in the discretised form as with #N being the total number of neighbours in the neighbourhood of the point P a . In a similar fashion, the inte-gral (28) 2 can be formally discretised as where V 2 is the the effective volume squared contributing to two-neighbour interactions defined by with #N 2 the total number of contributing pairs in the neighbourhood of the point P a .  : Schematic illustration of how a "contributing pair" is defined. The pair i, j in the upper half of the figure is a "contributing pair" and contributes to the energy due to two-neighbour interactions. However, the pair i, j in the lower half is not a "contributing pair" and does not contribute to energy.

Remark 2.
In view of the definition of contributing pairs, the first condition must hold since if the points {a, i, j} are collinear, the stiffness matrix can become singular. Furthermore, the derivations of the governing equations of CPD in [71] require that if a triangle { ai j} contributes to the energy, both triangles { jai} and { i ja} must also equally contribute to the energy which leads to the second condition. That is, if the distance between each pair is not less than or equal to the horizon size δ 0 , the stiffness matrix can lose its symmetry. Similar arguments hold for contributing triplets are defined next.
Finally, the discretised form of the integral (28) 3 reads with where V 3 is the effective volume cubed contributing to three-neighbour interactions in the neighbourhood of the point {k, j, i} contribute equally to the volume. To test the fidelity of the implementation, one can numerically compute the following integrals to ensure they hold exactly:

Discretised balance of linear momentum
The underlying governing equation of CPD is the linear momentum balance. The term b ext 0 in the linear momentum balance (8) corresponds to externally prescribed body forces and its incorporation into our framework is fairly straightforward and is omitted from this presentation in order to focus on the novel aspects of the computational implementation. The point-wise, non-local, form of the linear momentum balance (8) in the absence of body forces, i.e. the equilibrium equation, is thus given by We proceed with this reduced linear momentum balance and develop a discretised version based on the strategy presented in Section 3.1. Using the definitions of the force densities (17) at each collocation point, the integral over the horizon in (36) can be decomposed into three parts, corresponding to one-, two-and three-neighbour interactions, as follows This form of the point-wise balance of linear momentum can be expressed as where R is the point-wise residual vector and is decomposed into R 1 , R 2 and R 3 corresponding to one-neighbour, two-neighbour and three-neighbour contributions, respectively. That is, Next, we discretise the residual vector R. The global discretised residual vector R is composed of point-wise discretised residual vectors R a assembled into a global vector as follows The point-wise discretised residual vector R a of collocation point P a is comprised of the point-wise discretised residual vectors R a 1 , R a 2 and R a 3 corresponding to one-neighbour, two-neighbour and three-neighbour interactions, respectively. That is In the definitions of the point-wise discretised residual vectors (41), the bond vectors are related to the deformation via the relations with x a being the position vector of the collocation point P a in the deformed configuration. The deformation vector of the neighbours i, j and k are denoted by x i , x j and x k , respectively. The global residual vector R is energetically conjugate to the global deformation vector x that consists of the point-wise deformation vectors x a , that is As it is customary in nonlinear problems involving large deformations, the full deformation history can be divided into increments. The fully discrete nonlinear system of governing equations at each increment can be concisely stated as R · = O whose approximate solution is obtained via an iterative Newton-Raphson scheme. The consistent linearization of the resulting system at iteration k reads with the resulting deformation change ∆x k at iteration k given by where K k denotes the corresponding algorithmic tangent (stiffness) at iteration k. Finally, the deformation x is updated after each iteration by ∆x obtained from Eq. (45) according to Note that the tangent K is a matrix and composed of point-wise contributions K ab . That is Each contribution K ab itself can be further decomposed into the contributions from one-neighbour, two-neighbour and three-neighbour interactions. That is The next task is to compute the discretised point-wise residual R a and the discretised point-wise stiffness K ab from the stored energy densities (23), (24) and (25) corresponding to one-neighbour, two-neighbour and three-neighbour interactions, respectively.

Discretised residuals and tangents
The final steps in the computational implementation of the proposed scheme are (i) to express the discretised residual vectors (41) in terms of the deformation of the point P a and its neighbours, and (ii) to compute their associated tangents (48). We begin with the residual vectors. Before proceeding, recall the definitions in the material and spatial configuration, respectively. Furthermore, the following relations will prove useful throughout the forthcoming derivations: The point-wise discretised residual vector of point P a due to one-neighbour interactions reads where ∂ψ 1 | /∂ξ | can be expressed as Upon using Eq. (50) 1 and the relation S 1 = l /L = |ξ | |/|Ξ | |, one obtains Thus, the point-wise discretised residual vector of point P a due to one-neighbour interactions reads The point-wise discretised residual vector of point P a due to two-neighbour interactions is given by where and using Eq. (50) 2 and the relation Using the identity the point-wise discretised residual vector of point P a due to two-neighbour interactions reads Lastly, the point-wise discretised residual vector of point P a due to three-neighbour interactions is rewritten using the relation and using Eq. (50) 3 and the relation Inserting Eq. (60) into Eq. (59) yields the point-wise discretised residual vector of point P a due to three-neighbour interactions as Equipped with the discretised residuals (53), (58) and (61), the algorithmic tangents are derived next. Unlike the commonly accepted strategy in classical state-based peridynamics, we do not approximate the tangent stiffness using finite difference or other numerical differentiation schemes. A key feature of the proposed methodology is that we compute the tangent stiffness K directly. This is mainly possible since we do not rely on the notion of "state". Computing K directly has enormous advantages. For example, the associated decrease in the number of Newton iterations required to achieve convergence can reduce the computational time and significantly boost the accuracy of the calculations. Throughout our numerical simulations we observe the asymptotic quadratic convergence associated with the Newton-Raphson scheme. Note, for highly non-linear sets of equations, numerical differentiation is inaccurate and ultimately unstable [80]. We derive the tangents for pairs of points P a and P b due to one-neighbour, two-neighbour and three-neighbour interactions separately and combine them additively according to Eq. (48). The discretised tangent matrix for the points P a and P b due to one-neighbour interactions reads which after using Eq. (53), can be written as To proceed, we use the chain rule and the relation where i is the identity tensor. Thus, the discretised tangent for the points P a and P b due to one-neighbour interactions The discretised tangent matrix for the points P a and P b due to two-neighbour interactions is which, after using Eq. (59), can be written as To proceed, we use the chain rule and the identities the derivations of which are omitted for the sake of brevity. Inserting relations (70) into (69), the discretised tangent for the points P a and P b due to two-neighbour interactions reads One may attempt to rewrite the tangent (71) in a more compact form since similar terms appear on different lines.
However, for the sake of clarity, we retain this expanded version as it immediately follows from the previous derivations. Finally, the discretised tangent matrix for the points P a and P b due to three-neighbour interactions can be expressed using the relation To proceed, we once again use the chain rule and the identities the derivations of which are omitted for the sake of brevity. Inserting relations (75) into (73), the discretised tangent for the points P a and P b due to three-neighbour interactions reads Again, since similar terms appear on different lines of the tangent (76), it could be written in a more compact form.
However, for the sake of clarity, we retain this expanded version.
The stiffness components due to one-neighbour (66), two-neighbour (71) and three-neighbour (76) interactions, are all symmetric for the variationally consistent mechanical problem of interest here wherein the residuals derive from a potential. While in some cases, such as the first and second lines of Eq. (71), the symmetric structure of the stiffness is obvious, in other cases, such as the third and fourth lines of Eq. (71) the symmetry is not evident. The reason for this is that, for instance, the term ξ | ⊗ ξ || is calculated twice for any given set of {i , j } for which i = i and j = j for the first time but i = j and j = i for the second time. This property could indeed allow us to reduce the number of loops over the neighbours and avoid multiple counting. For instance, the index j instead of counting from 1 could start from i + 1 if the associated implications are accounted for. However, for the sake of brevity, we do not include these additional steps in the presentation. This procedure is essentially a technical programming detail also relevant in molecular dynamics simulations and associated methods. An efficient implementation to avoid multiple counting reduces the number of loops over the neighbours by a factor of 2 in two-dimensional simulations and by a factor of 6 in three-dimensional simulations.

Numerical examples
The

Convergence and non-locality in CPD
The main goal of this section is to investigate the convergence behaviour and the inherent non-locality of the proposed framework. This numerical example is carried out at small deformations in order to focus on the main objectives of this study. A schematic of the investigation and a depiction of the horizon size δ and grid-spacing ∆ is illustrated in Fig. 3. Furthermore, for all the computations in this section we utilize only one-neighbour interactions, for the sake of clarity. We have carried out similar studies including, in addition, two-neighbour interactions and have observed similar trends. The influence of two-neighbour interactions as well as three-neighbour interactions are studied separately in Section 4.3 and 4.4.
The details of the two boundary value problems studied are shown in Fig. 4 with an exaggerated deformation depicted. For both examples, the grid points are uniformly distributed with horizontal and vertical spacing ∆. In the Figure 3: Schematic of a convergence study versus a non-locality study. For the convergence study, the horizon size δ remains constant while the grid-spacing is decreased. For the non-locality study, the ratio of the horizon-over-grid size δ/∆ is fixed while the horizon is decreased to obtain a more local solution.
first row of Fig. 4, a solid unit square is extended by 0.1% in the horizontal direction. In the second row, 40% of the central region of the specimen is removed and the domain again extended by 0.1% in the horizontal direction.
For each of the two domains, we carry out two separate studies. For the first we fix the horizon size δ and decrease the grid-spacing ∆ resulting in more neighbours within the horizon of each point. This study demonstrates that, by increasing the number of grid points within the horizon, the solution of CPD converges, as expected. In other words, for a given horizon size δ, decreasing the grid-spacing ∆ results in a more accurate solution. In the second set of examples, the number of neighbours within the horizon remains constant but the horizon size itself varies. More precisely, we change δ for a given δ/∆ to study the non-locality associated with CPD. We expect to observe decreased non-local effects with diminishing horizon size δ. In the limit of δ → 0, we expect the solution of CPD to converge to the solution of CCM. We include the solution from CCM, obtained using the finite element method, for comparison. The vertical displacement of the full specimen decreases from zero to an extremum in the middle and are symmetric about both AA and BB . In addition, the vertical displacement over the line BB is greater than over AA , due to its increased distance from the centreline. The vertical displacement behaviour of the specimen with a square hole is more Figure 4: Schematic illustration of the specimens and the prescribed deformation. The first specimen is a full unit square and the second specimen is a unit square with a square hole at its centre. A lateral extension of 0.1% is applied to the specimens. The exaggerated deformed shapes are depicted schematically on the right.
complicated. Over the line BB , we observe a decrease from the edge to the middle of the domain, whereas the vertical displacement along AA increases at the beginning and then decreases to an extremum in the middle. This trend is not an artifice of our CPD formulation and is also observed in CCM. Moreover, the vertical displacements along both lines exceed those obtained for the full specimen, which again can be explained due to the more pronounced Poisson effect for a specimen with a hole at its centre.
In the second set of numerical studies, the horizon-over-grid size δ/∆ is fixed and the horizon size δ varies to highlight the non-local nature of CPD. Figure 6 shows both the vertical and horizontal displacements throughout the specimens. To facilitate comparison, the problem is designed to have the same features as the first numerical study depicted in Fig. 5.  we observe a decrease in non-local effects with diminishing horizon size δ and asymptotically, in the limit of δ → 0, the solution of CPD converges to the solution of CCM.

Properties of the stiffness matrix
This section elaborates on properties of the stiffness matrix K such as sparsity and symmetry. Figure 7 shows the sparsity pattern of the stiffness matrix for the full specimen and the specimen with a square hole. Two different gridspacings ∆ = 0.02 and ∆ = 0.05 are examined and for each case several horizon-over-grid ratios δ/∆ are considered.
The horizontal and vertical axis in each figure corresponds the columns and rows of the stiffness matrix, respectively.
As expected, the stiffness matrix is symmetric for all the cases. The colours in Fig. 7  There is a significant difference between the ranges for these two cases indicating that the stiffness matrix is considerably larger when ∆ = 0.02, which is intuitive since a smaller grid-spacing ∆ translates to more points, and hence more degrees of freedom.

Interplay between Poisson effect and the material parameters
This section provides a detailed comparison between CPD and CCM with a focus on the Poisson effect and its relation to the material parameters associated with each theory. The material parameters of CCM are the Lamé parameters λ and µ while for CPD they are C 1 , C 2 and C 3 . The parameters C 1 , C 2 and C 3 correspond to one-neighbour, two-neighbour and three-neighbour interactions, respectively. Conceptually, C 1 , C 2 and C 3 can be interpreted as resistance against the change of length, area and volume, respectively. Both two-and three-dimensional analyses are carried out and the significance of C 2 and C 3 for two-dimensional and three-dimensional problems explained. To perform this study, the specimen is subject to an extension at small deformations and the effective Poisson ratio is calculated. It should be emphasised that the effective Poisson ratio for both CPD and CCM is treated as a geometrical feature computed via numerical simulation. That is, the Poisson ratio is calculated by dividing the lateral contraction by the extension at the centre of the domain. Note that for the CPD analysis, the specimen is discretised into grid ν ν  Remark 3. It is important to recall that both CPD and CCM require three constants for isotropic elasticity at finite deformations. However, the linearisation process at small strains reduces the number of independent parameters from three to two in CCM. The CPD formalism accounts directly for finite deformations and is not a linearised theory, and hence the three constants are present. It would be feasible to establish a linear CPD theory in which only two independent parameters contribute to the material behaviour.   auxetic regime. For the CPD analysis, the horizontal axis represents the ratio of the two-neighbour elastic coefficient to the one-neighbour elastic coefficient C 2 /C 1 whereas for the CCM analysis it corresponds to the ratio of the first to the second Lamé parameter λ/µ. The parameters C 1 and µ are set to 1 and we vary C 2 and λ to generate the desired range for the ratios. Setting C 2 = 0 implies that two-neighbour interactions do not contribute. Thus only the one-neighbour interactions of CPD are active corresponding to bond-based peridynamics. It is observed that a Poisson ratio of 1/3 is obtained in CPD when C 2 = 0. This is well-known for bond-based PD [58]. A Poisson ratio of 1/3 is obtained in CCM when λ = µ as expected according to the relation ν = λ/[λ + 2µ] for two-dimensional local (linear) elasticity. Increasing C 2 or λ results in a higher resistance to a change of area, and hence we approach the incompressible limit. In the limits of C 2 → ∞ and λ → ∞, the two-dimensional incompressibility limit ν = 1 is obtained. It is observed that although there is a one-to-one relation between the graphs on the right and the ones on the left, for a given Poisson ratio, C 2 /C 1 does not coincide exactly with the same value for λ/µ. However, the Poisson ratio corresponding to C 2 /C 1 = 0 is identical to that associated with λ/µ = 0. Figure 9 shows the Poisson ratio versus the material parameters for a three-dimensional domain. The horizontal axis for CPD, in contrast to the two-dimensional problem presented previously, corresponds to the ratio of the threeneighbour to the one-neighbour elastic coefficient C 3 /C 1 . The horizontal axis for CCM, similar to the two-dimensional problem, is the ratio λ/µ. Similar to the previous case, the parameters C 1 and µ are set to 1 and C 3 and λ are calculated according to the desired range for the ratios on the horizontal axis. Note that in order to examine the effect of threeneighbour interactions, the two-neighbour elastic coefficient C 2 is set to zero here. In CPD, the Poisson ratio 1/4 is recovered when C 3 = 0; this corresponds exactly to the bond-based PD, as expected. In CCM, the same Poisson ratio is obtained when λ = µ which agrees directly with the classical relation ν = λ/[2λ + 2µ]. Similar to the twodimensional study, in the limits of C 3 → ∞ and λ → ∞, the incompressibility limit is obtained. The incompressibility limit now corresponds to ν = 0.5 due to the three-dimensional nature of the problem. As demonstrated in these two examples, our methodology can not only exactly recover bond-based peridynamics but it is also capable of covering the whole range of possible Poisson ratios including the auxetic regime. Figure 10: Schematic illustration of the two finite deformation example problems. The first specimen is a unit square and the second a unit square with a square hole at its centre. Both specimens are subject to an extension of 100%.

CPD at finite deformations
As shown in the previous section, CPD is capable of capturing any Poisson ratio by incorporating two-neighbour and three-neighbour interactions. In this section, we carry out a series of computational studies at finite deformations to compare the influence of multi-neighbour interactions on the material response for both two-dimensional and three-dimensional frameworks. Furthermore, we demonstrate the robustness of the proposed framework. For the twodimensional analysis, two different specimens are considered with the geometry and loading conditions illustrated in Fig. 10. Both specimens are subject to 100% extension in the horizontal direction and are free in the lateral direction.
We consider a unit square without and with a square hole at its centre as shown in the first and second rows of Fig. 10, respectively. The hole is in the shape of a square with the sides of length 0.4. For the three-dimensional analysis that follows, we mimic similar conditions and dimensions. First we detail the numerical simulations for a two-dimensional domain.
The discretised specimens used in the analysis are depicted in Fig. 11 together with their deformed shapes obtained Figure 11: Illustration of the original specimens and their deformed shapes for the two dimensional example at finite deformations. Zoom boxes are included to elucidate the distribution of the points throughout the specimens. Figure 12: Large deformations of a unit square without and with a square hole at its centre for different C 2 /C 1 ratio associated with different levels of incompressibility.
from computational simulations. The shading corresponds to the vertical displacements and zoom boxes are provided for further clarity. Figure 12 illustrates the deformed discrete domain for both specimens and for various values of C 2 /C 1 . The transparent shapes depicts the undeformed configuration. Prescribing displacement-type boundary conditions in CPD is similar to PD. That is, we prescribe the displacements for a few layers of points on the boundary as shown. The number of the layers over which the boundary conditions are imposed depends on the horizon size.
In the first row of Fig. 12, C 2 = 0 and thus, only one-neighbour interactions are active. The remaining rows involve both one-neighbour and two-neighbour interactions where C 2 /C 1 provides a measure of incompressibility and hence, given the loading and geometry, larger values of C 2 /C 1 lead to increased lateral contraction. Figure 13 presents a schematic illustration of the specimens and the prescribed deformations for the threedimensional analysis. Similar to the two-dimensional analysis, both specimens are subject to 100% extension in the horizontal direction and are free in both lateral directions. We consider a unit square without and with a cubic hole at its centre as shown in the first and second rows of Fig. 13, respectively. The hole is also in the shape of a cube with sides of dimension 0.4. Figures 14 and 15 show the deformed shapes as well as the distribution of the vertical displacement on both domains. Three different interaction types are examined. The first column is associated with one-neighbour interactions. Figure 13: Schematic illustration of the two study cases undergone and the applied deformation type. The first specimen is a full unit cube and the second specimen is a unit cube with a cubic hole at its centre. Both specimens are subject to 100% extension.
The results in the first column provide an excellent reference against which to compare the results in the second and third columns. The second column takes one-neighbour and two-neighbour interactions into account whereas the third column accounts for one-neighbour and three-neighbour interactions. These simulations are devised to distinguish between two-neighbour interactions and three-neighbour interactions. We emphasise that it is not possible to have only two-neighbour interactions or only three-neighbour interactions. This important observation will be contextualised later via a geometrical example. In the first row, the reference configuration is illustrated where the colours correspond to the vertical displacement. The deformed shapes are illustrated in the subsequent rows. On the second row, both deformed and undeformed configurations are shown. The smooth colour pattern is produced by interpolating the displacements between the points. Hence, the third row is more representative of what we obtain from the simulations as no interpolation between points is assumed. In the last row the deformed set of points on a vertical plane is extracted to highlight the difference between the three interaction types. Since only one-neighbour interactions are considered in the 1st column, the least contraction is obtained which leads to increased compressibility when compared to the second and third columns. It is important to note that although both two-neighbour interactions and three-neighbour interactions result in decreased compressibility, the solutions are different indicating that they indeed correspond to a different deformation.
As mentioned, one-neighbour interactions are necessary for the stability of the computation. In addition to the required one-neighbour interactions, two-and three-neighbour interactions can be included. However, the opposite cannot be done. More precisely, we cannot run simulations in which only two-neighbour or three-neighbour interactions are active and one-neighbour interactions are absent. That is, while two-neighbour and three-neighbour interactions are important to capture the Poisson effect correctly, one-neighbour interaction are necessary to avoid instabilities. The requirement for one-neighbour interactions is a physical one. Consider Fig. 16. In the upper part of the image three triangles are shown with equal area. In the lower half three pyramids are depicted with identical volumes.
Although the areas of the various triangles and the volumes of the various pyramids are the same, the length of their edges differ among them. Hence in the absence of one-neighbour interactions one could have multiple configurations with identical energies that would in turn cause instabilities. However, it is not possible to change a configuration without changing the distances between the pairs of points. Thus one-neighbour interactions guarantee that different configurations correspond to different energies, which does not hold for only two-neighbour and three-neighbour interactions.
Finally, to illustrate the robustness of the proposed framework and its ability to simulate large deformations, we study the convergence behaviour. Table 3 gathers the convergence of the normalized residual for the three-dimensional simulations at large deformations shown in Fig. 14 for various interactions at different increments. The domain in this   study is subject to 100% extension in 25 increments and the 2 -norm of the residual |R| recorded every fifth increment.
Despite the large deformation, consistent quadratic convergence is obtained consistently at every increment. A key feature of the proposed methodology is that it is fully implicit. Furthermore, the tangent stiffness is computed directly and not via numerical differentiation schemes, unlike the commonly accepted strategy in classical statebased peridynamics. The numerical implementation and solution procedure is robust and shows the asymptotically quadratic rate of convergence associated with the Newton-Raphson scheme. For the first time, specific constitutive laws for CPD have been presented. Their numerical implementations has been discussed in detail together with analytical forms for their associated tangents. The utility and reliability of the proposed strategy is illustrated via a broad range of numerical examples including three-dimensional problems at large deformations.
One logical extension of this work is to account for elasto-plasticity. Another interesting extension of CPD, and our next immediate plan, is to consider thermomechanical problems.