{'fraction_non_alphanumeric': 0.03447635409653379, 'fraction_numerical': 0.011584983632437955, 'mean_word_length': 4.399774351259872, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the treatment of boundary conditions for bond-based peridynamic models

In this paper, we propose two approaches to apply boundary conditions for bond-based peridynamic models. There has been in recent years a renewed interest in the class of so-called non-local models, which include peridynamic models, for the simulation of structural mechanics problems as an alternative approach to classical local continuum models. However, a major issue, which is often disregarded when dealing with this class of models, is concerned with the manner by which boundary conditions should be prescribed. Our point of view here is that classical boundary conditions, since applied on surfaces of solid bodies, are naturally associated with local models. The paper describes two methods to incorporate classical Dirichlet and Neumann boundary conditions into bond-based peridynamics. The first method consists in artificially extending the domain with a thin boundary layer over which the displacement field is required to behave as an odd function with respect to the boundary points. The second method resorts to the idea that peridynamic models and local models should be compatible in the limit that the so-called horizon vanishes. The approach consists then in decreasing the horizon from a constant value in the interior of the domain to zero at the boundary so that one can directly apply the classical boundary conditions. We present the continuous and discrete formulations of the two methods and assess their performance on several numerical experiments dealing with the simulation of a one-dimensional bar.


Introduction
In recent years, novel non-local modeling theories have been developed for the simulation of problems in structural and fracture mechanics. Peridynamic theory [1], which can be viewed as a generalization of classical continuum mechanics, provides for instance such a class of non-local models. Peridynamics does not involve differentials of displacement fields, which makes it an attractive framework for the modeling and simulation of fracture mechanics applications. Validation of the peridynamic theory has been undertaken using experimental data from a variety of wave propagation applications as well as crack initiation and propagation experiments, see e.g. [2]. However, major issues identified in that review include the socalled surface or skin effect [3,4,5] and the manner by which boundary conditions should be modeled within a non-local model. This was recognized early on during the development of peridynamics as Silling and Askari in [6], for instance, explained that displacement and traction boundary conditions should be prescribed within a layer of finite thickness under the surface of a solid body. Unfortunately, such boundary data cannot be assumed to be known in practical applications of interest. Several methods have since been proposed, see e.g. [7,8,9,10,11,12,13], in order to circumvent the lack of data and correct the spurious surface effects that are created as a result in the vicinity of the boundaries. Starting from the fact that data is only available at the surface of a solid body, we believe that boundary conditions for peridynamic models should be prescribed in a manner similar to that for local classical continuum models. We study in this paper two approaches to incorporate classical Dirichlet and Neumann boundary conditions into a one-dimensional linearized bond-based peridynamic model [1]. Our main objective in this paper is to mathematically analyze the resulting models, and more specifically, to verify that their solutions converge to that of the compatible linear elasticity model with respect to the so-called peridynamic horizon and discretization parameters.
The first approach, referred here to as the extended domain method (EDM), relies on the idea to artificially extend the domain by a thin layer of the size of peridynamic horizon. The idea is not new and is similar to that proposed by the authors in [14,15,16,17,18,19], among others. In [14,15,16] though, the thin layer was introduced within the solid body and subjected to Dirichlet-like non-local boundary conditions, that is, the displacement field in the layer was either arbitrarily set to zero or was constrained to match the values of the exact solution in the case of a manufactured problem. Zhou and Du [20] proposed to constrain the solutions to be either odd or even from either side of the boundary, in order to mimic homogeneous Dirichlet or Neumann boundary conditions, respectively. We propose here to constrain the displacement field in the artificial layer to satisfy odd extensions with respect to the boundary points so that one can handle both Dirichlet and Neumann non-homogeneous boundary conditions.
The second approach, referred here to as the variable horizon method (VHM), proposes to vary the peridynamic horizon from a constant value in the interior of the domain to zero as one approaches the boundary. The approach was considered in [3,21] for the treatment of boundary conditions along with adaptive mesh refinement. The method is motivated by the fact that the local model and non-local model are consistent and compatible in the limit that the horizon vanishes, so that one can justify the use of classical Dirichlet and Neumann boundary conditions to be prescribed at the boundary of the domain. Varying the horizon size necessitates however a scaling of the material parameter for the peridynamic model. We propose here an approach similar to [3,21] but with a scaling in 1D that appears to be different from that suggested in [3]. We also mention that the variable horizon approach was utilized in [15,22,23] for the treatment of interface problems in which different peridynamic models or a continuum model and a peridynamic model are coupled together.
Performance and comparison of the two methods are investigated on several numerical examples involving a one-dimensional linear elastic bar. In particular, we carry out a δ-convergence analysis [3,16] by considering four test cases in which the manufactured solutions are linear, quadratic, cubic, and quartic, in order to put in evidence the differences and similarities between the methods. The EDM and VHM solutions are also compared with the solution of the local linear elasticity model (LLEM). All models are discretized here using the finite difference method in order to allow for consistent comparison and are assessed using the maximum relative error with respect to the exact solution of the linear elastic local problem. In all cases, we verify that the numerical rates of convergence at which the peridynamic solutions converge to the linear elasticity solution are consistent with the theoretical rates. We also present additional examples in order to study the effect of a correcting factor for the extended domain method, to carry out a m-convergence analysis, and to consider the case of a solution with a very steep gradient in the vicinity of one boundary point.
The paper is organized as follows. Following the introduction, we present in Section 2 the model problem based on the classical linear elasticity theory and corresponding linearized bond-based peridynamic model for the simulation of a one-dimensional bar held fixed at one end and subjected to a traction at the other end. We proceed in Section 3 with the description of well-known spurious effects near boundaries when using non-local models and assess a correction method to alleviate the issue. We introduce in Section 4 and Section 5 the continuous and discrete formulations of the extended domain method and the variable horizon method, respectively, and assess their performance in Section 6 on several onedimensional numerical examples. We complete the paper in Section 7 with concluding remarks and perspectives.

Model problem and preliminaries
The model problem consists in finding the static equilibrium of a bar, of unit length and constant crosssectional area, held fixed at one end and subjected to a longitudinal traction at the other end. We shall suppose here that the deformations in the bar are infinitesimally small and can be adequately described by the theory of linear elasticity. We present below the local model based on classical elasticity theory and the non-local model counterpart built upon the linearized bond-based peridynamic theory [1].

Continuum local model
Let Ω = (0, 1). The continuum local problem consists in finding the displacement u in the bar such that: where E and A are the constant modulus of elasticity and cross-sectional area of the bar, respectively, g ∈ R is the traction force applied at end point x = 1, and the scalar function is the external body force density (per unit length). We suppose here that f b is chosen sufficiently smooth so that regularity in the solution is not an issue when comparing the solutions from the two models. Note that we use the notation u for the solution to the classical elasticity problem in order to emphasize that it may be different from the peridynamic solution u introduced below.

Peridynamic model
Let δ > 0 denote the so-called horizon of the peridynamic model and let H δ (x) = (x − δ, x + δ) be the subdomain of the neighboring particles within the horizon. We observe that for any given point x in the interval Ω δ = (δ, 1 − δ), we have that H δ (x) ⊂ Ω. We also introduce the boundary domains B δ 0 = (0, δ) and Figure 1. The general formulation, in one or higher dimension, of the linearized microelastic model [1] is given by: where κ is the parameter that characterizes the stiffness of the "bonds" between point x and the neighboring points y ∈ H δ (x), ξ is the vector between two material points in the reference configuration, i.e. ξ = y − x, ξ is the Euclidean norm of vector ξ, and u(x) is the displacement of x in the deformed configuration. In the case of the one-dimensional bar, the above integral at a point x ∈ Ω δ can be rewritten as: Ideally, the material parameter κ could be identified from experimental data. However, one may be interested in the value of κ such that the solution of the linearized microelastic brittle peridynamic model is the same as that of the continuum local model. One may do so by matching the strain energy of the two models [1]. Alternatively, one can try to recover equation (1) from (4) when taking the limit δ → 0. Supposing that the displacement u(x) is sufficiently Figure 1: Definition of the computational domains for the peridynamic model.
smooth, one can write, for all y ∈ Ω δ : so that the integral in (4) reduces to: Substituting the integral in (4) with the above result yields: By taking the limit when δ → 0, one then recovers the differential equation (1) pointwise whenever κ is chosen as: in agreement with [24]. In other words, the peridynamic model can be used to describe the behavior in the bar as one is able to identify the value of parameter κ that makes it compatible to the linear elasticity model. The main and well known issue with peridynamics is whether the model can also be employed in the region Ω\Ω δ , i.e. in the region of size δ along the boundary. A solution to this issue could be to transpose the Dirichlet condition (2) and Neumann condition (3), prescribed at the boundaries of the bar, into volumetric boundary conditions to constrain the material in B δ 0 and B δ 1 . However, it is still unclear how this can be achieved, and if it were, how to identify the corresponding boundary data. The point of view that we will follow here is that the notion of boundary conditions, which are prescribed on the surface of solid bodies, is naturally associated with local models. Our goal in the following section is to present some approaches to adapt the Dirichlet and Neumann conditions to the peridynamic problem.

Corrected forces near free surfaces
In this section, we study the integral in (4) at a point in B δ 0 , i.e. near the left boundary of the bar. The analysis is similar for a point in B δ 1 near the other boundary. Let x ∈ B δ 0 . Since H δ (x) extends beyond Ω, the integral becomes: (10) where we have used the change of variable ξ = y − x. When comparing the above result to (7), it is clear that truncation of the set H δ (x) close to the boundary induces forces involving the first derivative of u. This term vanishes in the case of free surfaces since zero traction, i.e. EAu = 0, is prescribed on the boundary. Moreover, the coefficient of u is different from that in (7). Since 0 ≤ x ≤ δ, the coefficient is actually smaller than for points in the interior of the domain, meaning that the material becomes softer as one approaches the free boundary, which is referred to as the so-called skin effect [3].
Several methods to correct the skin effect, such as the volume correction method [25,26], the force density method [27,18], the energy method [27,18], the force normalization method [28], or the positionaware method [29] have been compared in [5]. We just note that some of the methods seem somewhat heuristic and that the force normalization approach of Macek and Silling [28] appears to be valid only for two-and three-dimensional peridynamic models, as it involves an integral that is singular in the case of one-dimensional problems.
Looking at (10), a simple approach would be to replace the parameter κ by a position-dependent parameterκ(x) such that the coefficient of the second derivative matches that obtained for a point inside Ω δ , see (7), or equivalently, allows one to recover the linear elasticity model, that is: which leads to:κ It is notable that this result is equivalent to the correction factor that one obtains by the volume correction method. In fact, one would recover the same result by the energy method, at least in 1D. Indeed, the method consists in identifying the parameterκ(x) such that the strain energy computed at a point x near the boundary matches that computed at a point inside Ω δ for a given deformation. The strain energy per unit length of a "bond" is given by: It follows that the strain energy at a point x ∈ B δ 0 is provided by: For a point inside the domain, i.e. x ∈ Ω δ , the strain energy reads: Matching those two energies for a given deformation u naturally leads to (11).
Replacing κ byκ(x) in (10) and using (11), one obtains the new expression of the integral: Therefore, if u (x) = 0, such as near free surfaces, the peridynamic model leads to: In other words, due to the loss of symmetry of the integration domain, the peridynamic model provides an approximation of the classical linear elasticity theory of order one in δ in the vicinity of free surfaces, to be compared to O(δ 2 ) in the interior of the domain, see (8). Nevertheless, we emphasize here that the above correction approach is still unsuitable for problems involving Dirichlet boundary conditions or nonhomogeneous Neumann boundary conditions, since the first derivative of the solution does not necessarily vanish near the boundaries in these cases. We present in the next two sections alternative approaches to alleviate boundary effects, namely the extended domain method, if the horizon must be kept constant everywhere, and the variable horizon method, if the size of the horizon is allowed to vary in the domain.

Extended domain method 4.1 Continuous formulation
Let Ω δ = [−δ, 1 + δ] ⊃ Ω be the domain Ω augmented by the boundary domains D δ 0 = [−δ, 0) on the lefthand side and D δ 1 = (1, 1 + δ] on the right-hand side, see Figure 2. The goal here is to construct an extension of the solution u in Ω to the whole domain Ω δ . In the previous section, the integral in (10) can be viewed as the integral over [x − δ, x + δ] obtained by extending u to u(x) = 0, ∀x ∈ D δ 0 . The issue with such an extension is that it introduces a term involving the first derivative arising from the discontinuity in u at x = 0. We propose here, given a continuous function u in Ω, to consider odd extensions of u with respect to x = 0 and x = 1 as follows: We first note that the above conditions imply that the function u is continuous at x = 0 and x = 1 and satisfies the following properties: where x = 0 or 1, u (n) (x − ) and u (n) (x + ) denote the limits of u n when one approaches x from the left and from the right, respectively. In other words, the odd derivatives of u are continuous at x = 0 and x = 1 while the even derivatives are opposite at those points.
Using the above properties and Taylor expansions of u at x = 0, one can show, following lengthy calculations, that for x ∈ (0, δ): Figure 2: Definition of domains for the extended domain approach.
where we have used the facts that u (0 + ) = u (x) − xu (x) + . . . and 0 < x < δ. This is a remarkable result which shows that the leading term is now a term that involves the second derivative. The term containing the first derivative vanishes here thanks to the fact that the first derivative of the extended function is continuous at x = 0. If one had chosen an even extension, this term would not have cancelled.
Moreover, the peridynamic model converges to the classical linear elasticity model as δ tends to zero if the material property κ is replaced, for x ∈ (0, δ), by: and, similarly for x ∈ (1 − δ, 1), by: However, the convergence is only linear in δ. As an example, if We now formulate the peridynamic problem for the extended domain method. For the sake of simplicity in the presentation, we shall keep the value of κ constant over the whole domain, just keeping in mind that its value can be corrected in the intervals (0, δ) and (1−δ, 1) according to (18) and (19), respectively. The problem consists then in finding u ∈ Ω δ such that: where (9) has been used to modify the Neumann boundary condition (3) to (24). The fact that the displacement is extended in the boundary domains by odd extensions allows one to prescribe Dirichlet or Neumann boundary conditions. Indeed, if the Neumann condition at x = 1 is replaced by a non-homogeneous Dirichlet boundary condition, one would simply substitute (24) by: Although the solution to the problem (20)-(24) converges to the solution of the local linear elasticity model as the horizon δ approaches zero, the formulation nevertheless presents several drawbacks: A the problem is solved on a larger domain than the domain occupied by the bar, which increases the number of degrees of freedom for a given discretization parameter h. This is in particular true when δ/h is large and in higher dimensions; B odd extension of the displacement field outside of Ω may become difficult to carry on in two and three dimensions for domains with complex geometries, especially when boundaries involve corners and edges; C conceptually, there may still be an ambiguity in the application of the Neumann boundary condition, which is derived from the local model, while the behavior in the vicinity of the boundary is governed by the non-local model.
The approach proposed in Section 5 is an attempt at circumventing those issues.

Discrete formulation
We describe here the discretization of Problem (20)-(24) using a collocation approach. For a given δ, we introduce a uniform grid spacing h chosen such that, as it is customarily done in the literature [6,25], δ is a multiple of h, i.e. δ/h = m, with m a positive integer. For the sake of simplicity, but without loss of generality, we shall choose here m = 2. In this case, choosing h = 1/n, with n a positive integer, the discretization of domain Ω δ is determined by the grid points: where the points x −2 , x −1 , x 0 correspond to the grid in D δ 0 and the points x n , x n+1 , x n+2 to the one in D δ 1 . Moreover, let u i , i = −2, −1, . . . , n + 2, denote the discrete displacement at the grid points x i .
Approximation of the integral in (20) is obtained by classical quadrature formula using the grid points x i . Note that other more advanced quadrature rules, especially developed for non-local models, could alternatively be used, see e.g. [30,31]. We then have, using the second-order trapezoidal integration rule, or, in a more compact form, with α = κδ 2 /(16h 2 ), The constraints (21) and (22) to enforce odd extensions on the displacements read: while the Dirichlet boundary condition at x = 0 implies: There remains to discretize the Neumann boundary condition at x = 1. Several options based on the finite difference method are available. However, since the function u is only C 1 at x = 1, we cannot consider a second-order centered difference of the first derivative. We therefore opt for the second-order one-sided finite difference: that, in order to preserve the symmetric feature of the peridynamic model, can be rewritten using (30) and (31) in the centered form as: We emphasize that this is a second-order finite difference of u (1), whether u satisfies the odd extension or not. We propose to employ the latter formula to discretize the Neumann boundary condition, i.e.
The resulting system of linear equations associated with the discretization of the extended domain method is summarized in Figure 3. Upon inspection, it is worth noting that the use of the Dirichlet boundary condition at x = 0 allows one to eliminate the first equation of the system as the value of the degree of freedom u −2 , associated with the grid point x −2 in the extended domain, does not influence the other degrees of freedom. The same conclusion would also hold for u n+2 if one had considered a Dirichlet boundary condition at x = 1. We actually present in the next subsection an approach to further reduce the system of equations.

Reduced system of equations
In fact, it is possible in the one-dimensional case to implicitly enforce the odd extension constraints so as to reduce the system of equations and avoid solving for u −2 , u −1 , u n+1 , and u n+2 . Using (29) and (30), the fourth and (n − 1)-th equations in the system of Figure 3 become: whereas using the left-sided second-order approximation of u (1) to discretize the Neumann boundary condition gives: The reduced system of equations is shown in Figure 4. We emphasize here that the extended domain method and its reduced version provide exactly the same solutions for i = 0, . . . , n. We will therefore use the reduced system displayed in Figure 4 for solving the problems in the numerical examples. We observe that the matrix has in fact a similar structure as the one obtained using the second-order finite difference method for the discretization of the local model (see system of equations in Figure 7). We verify below that the equations associated with rows i = 1, . . . , n − 1 provide second-order approximations in h, and a fortiori in δ since h = δ/2 here, of u at the points x i . Starting with (26), we see that the expression in brackets, multiplied by the factor 1/h 2 , can be rewritten as: . . .  In other words, the second derivative u (x i ) is approximated by the weighted average of the secondorder centered differences of the second derivative using, on one hand, the nearest neighbors and, on the other hand, the next-to-nearest neighbors. We now turn our attention to (32). The term in brackets, using classical Taylor expansions, becomes: where the details of the calculation have been omitted. This result was to be expected in view of (17). Since x 1 /δ = 1/2 here, the factor 7/8 = 0.875 actually corresponds to an approximation of the fac-tor κ/κ(x 1 ) ≈ 0.9 that we obtained from (18). The discrepancy is simply an error due to numerical integration. One could multiply the expression by the factor 8/7 ≈ 1.14 in order to recover −u (x 1 ), or consider the material property at that point to bē κ(x 1 ) ≈ 1.1κ according to (18), but we note that the approach would still remain suboptimal due to the fact that one gets only a first-order approximation. Note that the same conclusion holds in the case of (33). We will illustrate the spurious perturbations near the boundaries in the numerical examples depending on the chosen correction method.

Continuous formulation
Our motivation in this approach is to address some of the ambiguities associated with the extended domain approach, the first one being that the problem is defined on a larger domain than the original domain Ω occupied by the bar (although we have shown than the problem from the extended domain method can be recast, in the one-dimensional case, as a reduced problem defined on Ω). Our objective is therefore to introduce a new problem restricted to the computational domain Ω. Secondly, in order to avoid constructing corrected forces near the boundaries, the set H δ (x) should always be a subset of domain Ω, suggesting that the horizon should tend to zero as one approaches the boundary. Finally, we suppose that one is interested in using a non-local model with constant horizon δ in the interior of Ω sufficiently far from the boundaries. In order to satisfy all these requirements, one needs to construct a non-local model with variable horizon δ v (x) such that, ∀x ∈ Ω: The definition of the function δ v (x) is therefore not unique. However, the simplest continuous function that fulfill the above requirements is the piecewise linear function defined as: The function δ v (x) is shown in Figure 5. Alternative functions, in particular smoother functions, could also be considered as long as they take values between 0 and δ v (x) (as defined in (37)) for all x ∈ Ω. In view of (9), the horizon being a function of x implies that the material parameter κ should also depends on x, i.e. κ =κ(x), if the non-local model is to be compatible with the local model. The new problem consists in searching for u(x) such that: In order for the model to be compatible with the linear elasticity model everywhere in the domain with a convergence of order O(δ 2 ), the productκ(x)δ 2 v (x) needs to remain constant for all x ∈ Ω, so that where κ satisfies (9). We note that this scaling of the material parameter appears to be different from that in [3].

Discrete formulation
As before, we shall use a collocation approach to discretize the problem (38)-(40) and choose the parameter h such that δ/h = m, with m = 2. The grid in Ω is defined by the grid points: and u i , i = 0, 1, . . . , n, denote the discrete displacements at the grid points x i . From (37), we observe that: The approximation of (38) at x 1 and x n−1 by a second-order integration scheme for the integral is given as: However, using relation (41) and the fact that, in those two cases, δ v (x i ) = h, one can rewrite the above equation as: or, using again the coefficient α = κδ 2 /(16h 2 ), In the case of the remaining interior grid points x i , i = 2, . . . , n − 2, the approximation of (38) is given by (27), that is: The Dirichlet boundary condition is prescribed as usual as: u 0 = 0.
Finally, it remains to discretize the Neumann boundary condition. A second-order approximation of the first derivative using a left-sided finite difference scheme is: Therefore, the approximation of the Neumann condition, using EA = κδ 2 /2 and parameter α, reads: 4αhu n−2 − 16αhu n−1 + 12αhu n = g.
The resulting system of linear equations, corresponding to the discretization of the variable horizon method, is summarized in Figure 6. We note here that the coefficient α could be factorized in front of the matrix. We conclude the section by emphasizing that the variable horizon method provides an approximation of the classical linear elasticity model of order two in the parameter h, and a fortiori, in δ. We will see that this rate of convergence is consistently observed in the numerical examples presented below.

Numerical examples
The objective of this section is to present several numerical examples in order to compare the solutions of the extended domain and variable horizon approaches for implementing boundary conditions in non-local peridynamic models. For completeness, we will also show the results obtained with the local linear elasticity model. The system of equations for the discretization of the problem (1)-(3) by a second-order finite difference method is shown in Figure 7, where we have assumed the same grid on Ω than the one defined in Section 5.2. For simplicity, but without loss of generality, we shall consider a one-dimensional bar for which EA = 1, so that κ = 2/δ 2 , α = 1/(8h 2 ), and β = 1/(2h 2 ). Moreover, the value of the boundary data for the Neumann boundary condition is set in all examples to unity, i.e. g = 1, unless explicitly said otherwise. The first series of examples will be concerned with δ-convergence (see e.g. [3,16]), for which we shall measure the error in the peridynamic solutions with respect to the classical local solution. In this case, we note that the discretization parameter h implicitly goes to zero as well since h = δ/2. We then provide an example to study the effect of the correcting factor in the material property for the extended domain method. The following example will treat the case of m-convergence, for which δ is fixed and h goes to zero, in order to show that the approximate solution actually converges to the corresponding peridynamic solution. Finally, we consider an example whose solution is given in terms of an exponential function and exhibits a very steep gradient at one end point. The objective is to show that, as expected, the rates of convergence are unaltered even in the presence of large variations in the solution.
. . . Figure 6: System of linear equations associated with the discretization of the variable horizon method where . . .

δ-convergence numerical study
We consider in this section four examples for which the manufactured solutions are polynomial functions of degree up to four. The objectives of this numerical study are to analyze the rates of convergence in δ as the peridynamic solutions converge to the linear elasticity solution and to highlight the differences between the methods. For comparison, we compute the maximum relative error: with n the total number of degrees of freedom of the numerical solution u i , i = 1, . . . , n, for the local linear elasticity model (LLEM), the extended domain method (EDM) without correction, and the variable horizon method (VHM).

Linear solution
The first example that we shall consider is defined when the body force density f b identically vanishes on Ω, i.e. f b (x) = 0, ∀x ∈ Ω, in which case the solution to Equations

Quadratic solution
The second example deals with the problem for which the body force density is given by: In that case, the exact solution to the local model equations (1)-(3) reads: We now solve the discrete problems with δ = 1/2, 1/4, 1/8, and 1/16, (i.e. h = 1/4, 1/8, 1/16, and 1/32, respectively) and report the maximum relative error in Table 1. We observe on the one hand that the relative errors for LLEM and VHM are identically zero independently of the grid size h, as expected. Indeed, since the derivatives of order three and higher of the exact solution u all vanish, the discrete methods are able in this case to exactly reproduce the solution.
On the other hand, we observe that the errors in the displacement when using EDM are nonzero and seem to decrease with a rate of convergence of order one with respect to h. In this case, the only sources of error arise from the fact that the second derivative is not properly estimated at x 1 = h and x n−1 = 1 − h, as explained by the result in Eq. (35). As a matter of fact, the error becomes zero if we multiply the left-hand side in the second and (n − 1)-th equations of Figure 4 by the correcting factor 8/7. Finally, we show in Figure 8 the error e i = u(x i ) − u i , i = 1, . . . , n, in order to illustrate how the sources of error committed at x 1 and x n−1 pollute the solution in the whole computational domain.

Cubic solution
We repeat the previous experiment with the following body force density:  for which the exact solution to Equations (1)-(3) is given by: The results are reported in Table 2. In the case of LLEM and VHM, we observe that the maximum relative errors decrease with a rate of convergence of order two with respect to h, as expected. Moreover, the errors for both methods are identical. This is explained by the fact that the source of error only arises from the approximation of the Neumann boundary condition. Indeed, the approximation being of order two, the leading term of the truncation error involves the third derivative of u, which is nonzero, while the leading term of the truncation error for the approximation of the second derivative involves the fourth derivative of u, which is in this case zero. Therefore, since the Neumann boundary condition is discretized in the same way for LLEM and VHM, the two methods provide the same approximate solution of the problem, hence the same relative error. As far as EDM is concerned, we note the same rate of convergence, approximately of order one, as in the previous problem with the quadratic solution. It is once again obvious from Figure 9 that the issue comes from the approximation of the second derivative at  slope of the displacement error at x n−1 . The issue is however less visible at x 1 as the solution is close to zero near the boundary.

Quartic solution
We repeat the previous experiment and consider the case of the quadratic body force density: which provides the following quartic solution to Equations (1)-(3):  Table 3: Maximum relative error E n for the quartic solution obtained with the local linear elasticity model (LLEM), the extended domain method (EDM), and the variable horizon method (VHM).
The maximum relative errors are shown in Table 3. We first note that the relative error for both LLEM and VHM decreases with a rate of convergence of order two similarly to the previous example. However, the errors in the solutions obtained with VHM are slightly greater than those in the solutions obtained by LLEM. Since the fourth derivative of u is no longer zero in Ω, the error is now partly due to the approximation of the second derivative. In the case of VHM, it is obtained as the average of two centered difference stencils of order O(h 2 ) and O((2h) 2 ), respectively (see e.g. Eq. (34)) at the grid points x i , i = 2, . . . , n − 2. It naturally follows that the solution is slightly less accurate than the solution obtained from LLEM, which only involves approximations of order O(h 2 ). However, the difference between the errors in the two solutions correspond to the modeling error of order δ 2 that arises from using the non-local model instead of the local model.
Finally, we show in Figure 10 the displacement error for the EDM solution. The conclusions that we drew in the cases of the quadratic and cubic solutions also hold, in addition to the fact that the behavior of the error seems even more unpredictable here. For comparison, we show in Figure 11 the displacement error in the VHM solution, which exhibits a monotonic and smoother behavior and decreases much faster as δ goes to zero.

Effect of correction factor in EDM
In the previous numerical examples, we considered   the extended domain method without any correction factor. In this example, we propose to study the effect of the analytical correction factor k(x)/κ, as given in (18) and (19) and shown in Figure 12 for δ = 0.125, and the numerical correction factor k(x 1 )/κ = k(x n−1 )/κ = 8/7 as explained by the result in (35). The correction methods with the analytical correction factor and the numerical factor shall be referred to here as EDM I and EDM II, respectively. We note that these correction factors are applied to the equations associated with node x 1 and node x n−1 only. We consider the problem with the quartic solution as described in the previous section and plot in Figures 13 and 14 the errors in the displacement when using EDM I and EDM II for various values of the horizon δ. We still observe some spurious effects in the solutions at x n−1 for both methods, as in EDM, but these tend to diminish as the horizon goes to zero. We compare in Figure 15 the absolute error in the displacement obtained with δ = 0.0625 for EDM without correction, EDM I, and EDM II. We observe that the maximum error in this case is reduced by a factor 5 for EDM I and a factor 10 for EDM II when compared to the EDM solution. Finally, the maximum relative errors are shown in Table 4. We observe that the rate of convergence for EDM without correction and EDM I is about one, as expected. Surprisingly though, the convergence rate for EDM II is approximately equal to two. We believe that, in this particular case, the solution benefits from favorable error cancellation as the relative error is even smaller than that obtained with LLEM and VHM, see Table 3 for δ = 0.125 and δ = 0.0625.
In conclusion, we have demonstrated that it is possible to improve EDM using the analytical correction factor or the numerical correction factor. The latter seems to yield better results, but its value is directly dependent on the numerical scheme and should be derived for each type of differentiation stencil. On the other hand, the use of the analytical correction factor, while improving on EDM without correction, still leads to errors that are larger than VHM. For that reason, we will consider only VHM in the following numerical examples.  Figure 13: Error in displacement using EDM I with correction factor k(x)/κ given by (18) and (19) for various values of the horizon δ.

m-convergence numerical study
The objective of this numerical example is to show that the approximate solution by the variable horizon method converges to the exact solution of the peridynamic model. We choose for the manufactured solution the quartic function of the previous example, for which the corresponding loading term in the peridynamic problem is given by: We consider here two values of the horizon, δ = 1/4 and δ = 1/8, and study the convergence in h, with h = δ/m and m = 2, 4, 8. The maximum relative errors E n and errors in the displacement field are shown in Table 5 and Figure 16 and in Table 6 and in Figure 17, for δ = 1/4 and δ = 1/8, respectively. We clearly observe that the solutions converge to zero with respect to the discretization parameter h at a rate of about two in both cases.

Manufactured solution with steep gradient
We consider in this example the manufactured solution to the classical linear elasticity problem: The corresponding load f b (x) was calculated using the computer algebra system Maxima and the boundary data g for the Neumann boundary condition is given by: We observe that the magnitude of the first derivative u at x = 1 increases as 1/ , as illustrated in Figure 18 in the cases where = 0.1 and = 0.01. Since the exact solution vanishes at x = 1, we evaluate in this example the maximum absolute error E n : rather than the maximum relative error as before. These errors are reported in Tables 7 and 8    with a rate of two for VHM with respect to the horizon δ.
We show the displacement obtained with VHM for various values of δ for = 0.1 and = 0.01 in Figures 19 and 20, respectively. It appears from these figures that the maximum error occurs mostly at the end point x = 1 and that a large number of degrees of freedom needs to be employed in order to recover the boundary layer. It implies that the source of discretization error at x = 1 severely pollutes the displacement field inside the domain.

Conclusion and future work
We have presented in this paper two methods, the so-called extended domain method and variable horizon method, to enforce boundary conditions within the non-local bond-based peridynamic model. The methods were developed based on the requirement that the solutions to the corresponding problems be compatible, as the horizon goes to zero, with the solution to the boundary-value problem derived from the classical linear elasticity model. The formulation of both methods was first derived at the continuous level, independently of the choice of the discretization scheme, and later discretized, as an example, using the finite difference method. The performance of both methods was assessed on a simple onedimensional linear elastic bar, fixed at one end and subjected to a traction at the other end, for various body force densities. The solutions of these manufactured problems were polynomial functions of degree one to four, which allowed us to study the advantages and disadvantages of each method.
The extended domain method (EDM) without correction provided the correct solution only in the case of the problem with a linear solution, for which the body force density is zero everywhere in the domain. In all other cases, we observed spurious artefacts in the displacement field near the boundaries. The reason for this behavior was explained by the fact that the extended domain method fails to correctly estimate the second derivative at the grid points closest to the boundaries, as it should in order to be compatible with the local model. We have seen that it was possible to improve the solutions obtained by EDM when using an analytical correction factor or a numerical correction factor. On the one hand, the analytical correction factor is derived at the continuous level and is valid for any type of discretization schemes. On the other hand, the numerical correction factor depends on the choice of the differentiation scheme. We observed on a numerical example that the numerical correction factor actually provided better results than the analytical correction factor. In view of this ambiguity, it is difficult to appreciate the superiority of one or the other correction approach.
The variable horizon method (VHM) delivered, for all the numerical examples, consistent solutions, if not the same, with those obtained from the local linear elasticity model. In that sense, VHM should be considered a better approach than EDM when applying boundary conditions in non-local models. It is also worth noting here that VHM can be viewed as a coupling method between a local model and a non-local model (see e.g. [32]), where, in that particular case, the coupling interface happens to be at the boundary of the domain.
VHM has shown its potential, through the preliminary one-dimensional numerical examples, for enforcing boundary conditions in the bond-based peridynamic model. In particular, we have theoretically and numerically shown that VHM is a method of order two in the horizon δ, with respect to the linear elasticity model, as well as in the discretization parameter h, with respect to the finite difference method. However, additional work should be carried out in order to confirm its performance in more general situations. For example, the method should be extended to twoand three-dimensional problems and assessed in the case of state-based peridynamic modeling. One could also consider other functions δ v (x) than the piecewise linear function studied here and analyze their effects on the resulting solutions. Finally, the method should be assessed when other discretization methods, e.g. the Finite Element method, are used for its implementation.