Accurate computation of partial volumes in 3D peridynamics

The peridynamic theory is a nonlocal formulation of continuum mechanics based on integro-differential equations, devised to deal with fracture in solid bodies. In particular, the forces acting on each material point are evaluated as the integral of the nonlocal interactions with all the neighboring points within a spherical region, called “neighborhood”. Peridynamic bodies are commonly discretized by means of a meshfree method into a uniform grid of cubic cells. The numerical integration of the nonlocal interactions over the neighborhood strongly affects the accuracy and the convergence behavior of the results. However, near the boundary of the neighborhood, some cells are only partially within the sphere. Therefore, the quadrature weights related to those cells are computed as the fraction of cell volume which actually lies inside the neighborhood. This leads to the complex computation of the volume of several cube–sphere intersections for different positions of the cells. We developed an innovative algorithm able to accurately compute the quadrature weights in 3D peridynamics for any value of the grid spacing (when considering fixed the radius of the neighborhood). Several examples have been presented to show the capabilities of the proposed algorithm. With respect to the most common algorithm to date, the new algorithm provides an evident improvement in the accuracy of the results and a smoother convergence behavior as the grid spacing decreases.


Introduction
The peridynamic theory provides a nonlocal reformulation of classical continuum mechanics: the internal forces are evaluated with integral equations, which are valid regardless of the presence of discontinuities in the displacement field. Hence, peridynamics can naturally model crack initiation, propagation and branching in solids. The first formulation of the peridynamic theory was the bond-based version [1], in which the Poisson's ratio is restricted to a fixed value. Subsequently, state-based peridynamics was developed [2], introducing the possibility of varying the Poisson's ratio. In the literature, there are many examples of applications [3,4], ranging from complex crack patterns, such as spontaneous branching [5], to multi-physics problems involving fracture [6,7].
Peridynamic points interact with each other up to a finite distance , called "horizon". The "neighborhood" of a point is the set of all the points interacting with that point. Therefore, the neighborhood has a circular shape in 2D problems and a spherical shape in 3D problems. The peridynamic equation of motion is based on the spatial integration over the neighborhood of the internal forces, which are generated by the interactions between neighboring points. In practice, the integration of the peridynamic equation of motion is carried out by means of numerical tools. The body can be discretized by a uniform or non-uniform grid (see for instance [8][9][10]). Various methods have been utilized to integrate numerically peridynamic equations: meshfree method with composite midpoint quadrature [11][12][13][14], Gauss-Hermite quadrature [12], finite element method [12,15,16], collocation method [17,18] and an adaptive integration method with error control [19]. Thanks to its simplicity of implementation and relatively low computational cost compared to other approaches, the meshfree method with a uniform grid is the most commonly used for peridynamic simulations. In 1 3 this method, the body is discretized in volume cells with a square shape in 2D problems and cubic in 3D, and the nodes lie at the center of the corresponding cells. The spatial integration over the neighborhood is transformed into a summation of integrals over cells and the midpoint quadrature rule is then applied in each cell, in which the nodes are employed as quadrature points.
However, near the boundary of the neighborhood, some cells are only partially within the neighborhood itself. Therefore, the quadrature weights related to those cells are computed as the fraction of cell volume which actually lies inside the neighborhood. The intersection area or volume of those cells with the neighborhood is also referred to as "partial area" in 2D and "partial volume" in 3D. The accuracy and convergence of the peridynamic results depend on the algorithm to compute the quadrature weights, i.e., the partial areas or volumes [13,14,20].
The first algorithm proposed in [11] considers the nodes within the neighborhood with their entire cell (even if a part of the volume is partially outside the neighborhood) and neglects the nodes outside the neighborhood (even if their cell is partially inside the neighborhood). This approach, under grid refinement, leads to an oscillatory convergence behavior in which the fluctuations seem rather random [13]. Many other algorithms to approximate the partial volumes have been proposed since then. An approximation based on the distance between neighboring nodes is proposed in [21], which improves the computation of the partial volumes of nodes within the neighborhood. A similar approximation is used in [22,23] to include the previously neglected partial volumes of nodes outside the neighborhood but with a part of the cell inside it. These algorithms reduce, but never eliminate, the seemingly random fluctuations of the convergence behavior. Subsequently, the algorithm to compute analytically the partial areas has been developed in [13]: the types of intersection between the neighborhood and the cells are rigorously categorized and then subdivided into domains of basic geometry (triangles, rectangles and circular segments), for which the analytical computation of the area is straightforward. Using this algorithm, the convergence behavior is smoothly oscillatory and the simulations yield results affected by a smaller error, on average, compared to the previously mentioned algorithms. In [13], it is also suggested to use the centroids of the partial areas as quadrature points, but this significantly increases the complexity of the computational model. The analytical computation of the partial volumes in 3D problems is much more complex and no algorithm is currently available for such purpose. The partial volumes can be computed numerically by two proposed algorithms, one based on the trapezoidal rule [19] and one based on a process of recursive subdivisions and sampling [14]. However, to reach the desired accuracy the computational cost may be very high. Other algorithms, specialized for non-uniform grids, are presented in [9,22,24,25].
The aim of this paper is to simplify the implementation of the algorithm to compute analytically the partial areas presented in [13] by skipping the step of subdivision of the intersection area in basic geometries, and to develop an algorithm for the analytical computation of the partial volumes. In order to achieve this, we solve directly the integrals which describe all the possible intersection areas or volumes. Actually, some integrals involved in the computation of the partial volumes are not explicitly solvable. Hence, we perform a Taylor series expansion of those functions and integrate the polynomials. The computation of the partial volumes converges to the analytical solution if the sum of infinite terms is not truncated. This, clearly, is not possible in a numerical model, but we will show that the algorithm is able to reach values of the error very close to machine precision with little computational effort. The numerical results obtained with the new algorithm show an evident improvement in the accuracy and in the convergence behavior when compared to the results obtained by the algorithm based on the approximation proposed by [22], which is arguably the most commonly used in engineering applications. We compared the numerical results by using the ordinary state-based version of the peridynamic theory, but the algorithm can be used with the bond-based version as well.
The paper is divided as follows. Section 2 reviews the basics of the state-based peridynamic theory and its discretized formulation. Section 3 presents the innovative algorithms for the computation of the partial areas and partial volumes. Section 4 contains several numerical examples that show the improvements provided by the proposed algorithm for the computation of the partial volumes with respect to the most commonly used algorithm. Section 5 draws the conclusions of the work.

Peridynamic theory
The peridynamic theory is a continuum theory based on nonlocal interactions between material points [1]. The derived numerical formulation is a very useful tool for simulating crack propagation in solid bodies. In the following, we present the fundamentals of the ordinary state-based peridynamics and the discretized formulae which could be implemented in a computational code.

Continuum model
The nonlocal interaction between two points, and ′ , in a peridynamic body B is described by a quantity named "bond": The relative displacement  vector is defined as where is the displacement field. Note that + is the relative position of points and ′ after the deformation occurred.
The state-based peridynamic equation of motion of a point within the body B is given by [2] where is the material density, ̈ is the acceleration field, is the force state, dV � is the differential volume of a point ′ within the neighborhood H and is the external force density field. The notation [ , t]⟨ ⟩ means that the state depends on the position of the point and on the time t, and operates on the bond . In an ordinary state-based peridynamic model, the force state is aligned with the corresponding bond for any deformation, as depicted in Fig. 1. For the purposes of the paper, it suffices to limit the study to quasi-static problems [26][27][28]. Hence, the peridynamic equilibrium equation is derived from Eq. (3) by dropping the dependence on time: The reference position scalar state x , representing the bond length, the extension scalar state e , describing the (1) ∶= � − , elongation (or contraction) of the bond in the deformed body, and the deformed direction vector state , the unit vector in the direction of , are, respectively, defined as The weighted volume m and the dilatation of a point are defined as where is a prescribed spherical influence function [29]. We adopt the Gaussian influence function Adopting the linear peridynamic solid model [2], the force state is computed as where K is the bulk modulus and is the shear modulus.
Since ⟨ ⟩ = − ⟨− ⟩ , the ordinary state-based peridynamic equilibrium equation becomes [2,30,31] Equation (12) relates the external forces to the displacement field, which might be a discontinuous function with respect to the spatial coordinates.

Discretized model
We adopt a meshfree method with a uniform grid spacing h to discretize the body domain (see, for instance, the discretization of a neighborhood in Fig. 2). Therefore, the cells surrounding each node are squares in 2D and cubes in 3D, respectively, with an area A = h 2 and a volume V = h 3 . Consider a node i as the node at which the peridynamic equilibrium equation should be computed and a node j with a portion of its cell inside the neighborhood H i of node i. The bond connecting node i to node j is described by Analogously, the relative displacement vector after the deformation of the body is defined as where i and j are the displacement vectors of nodes i and j, respectively.
Hence, the reference position scalar state and the influence function of bond ij can be computed as follows: Under the assumption of small displacements, the extension scalar state of bond ij is given as The non-local properties of node i, i.e., the weighted volume m i and the dilatation i , are determined numerically by transforming the integrals in Eqs. (8) and (9) into a summation of integrals over cells and applying a midpoint quadrature rule in each cell: where ij V represents the quadrature weight of the contribution of node j in the integral over the neighborhood of node i. The accurate computation of coefficients ij is the main result of the paper, which is presented in Sect. 3. In 2D problems under plane stress conditions, the volume of the cell is given as V = A t , where t is the constant thickness of the plate. Under the assumption of small displacements ( ⟨ ij ⟩ ≈ ij ∕‖ ij ‖ ), the peridynamic equilibrium equation in the discretized form is computed by using the quadrature scheme previously described as where m j and j are, respectively, the weighted volume and the dilatation of node j computed with Eqs. (18) and (19), and i is the external force density vector applied to node i.

Algorithms for the computation of the quadrature weights
The quadrature coefficient ij is the dimensionless factor defined as Similarly, in the discretized model of the neighborhood H i of a given node i, there are some nodes (solid circles) whose cell lies completely or partially inside the neighborhood and other nodes (empty circles) whose cell lies completely outside the neighborhood where Ṽ and Ã are, respectively, the partial volume and partial area of the cell, which should be computed from the intersection between the neighborhood H i and the cell itself. In particular, ij = 1 if the cell is completely inside H i , ij = 0 if the cell is completely outside H i and 0 < ij < 1 if the cell is partially inside H i . The set of nodes which constitutes H i ∶= j ∈ B ∶ ij > 0 depends on the algorithm used to compute ij .
We present one of the algorithms based on the approximation of the partial volume [22], used as a reference to compare our results. To compute analitycally ij in a convenient framework, we define a new reference system and exploit the cell-neighborhood symmetries. Then, we improve the algorithm for the analytical computation of partial areas by employing a simpler scheme and we use the same scheme to compute quasi-analytically the partial volumes. We utilize the expression "quasi-analytical" because the algorithm includes the truncation of the Taylor series expansions, but it is able to attain accurate results with a relatively small truncation order.

Approximated computation of partial areas or volumes
In the literature there are many algorithms that compute the partial areas or volumes as an approximation based on the distance between neighboring nodes [13]. The algorithm presented in [22] (see Algorithm 1) is arguably the most commonly used to date. This algorithm is based on the analytical computation of partial lengths in 1D problems, as shown in Fig. 3. If the distance between node i and the farthest side of the cell of node j is smaller than the horizon size, namely ‖ ij ‖ + h 2 < where h is the grid spacing, then ij = 1 (see Fig. 3a). If the distance between node i and the closest side of the cell of node j is greater than the horizon size, namely ‖ ij ‖ − h 2 > , then ij = 0 (see Fig. 3c). Otherwise, ij is computed as the difference between the horizon size and the distance of the closest side of the cell of node j from node i, divided by the grid spacing h (see Fig. 3b).
Algorithm 1 is applied without modifications to 2D and 3D problems. If the direction of ij lies along one of the axis, as for instance in Fig. 3b, the approximation is quite accurate. However, Fig. 4 shows other examples in which the computation of the quadrature weights with Algorithm 1 can be rather inaccurate. As a result, even if this algorithm is very simple, it leads to convergence issues [13,14]: for small variations of the grid spacing (considering fixed the horizon size) there are considerable variations of the computed mechanical properties, as the numerical results of Sects. 4.2 and 4.3 show.

Change of reference system
For simplicity sake, the concepts are hereinafter explained in the 2D case, but the generalization to a 3D case is straightforward. Since the focus is on the neighborhood of a node i, we adopt a new system of reference (x, y) with the origin at (x i , y i ) and the distances scaled by a factor 1/h, as shown  The coordinates of a node j in the new reference system are given as x j = (x j − x i )∕h and y j = (y j − y i )∕h , and the horizon size of the neighborhood becomes m = ∕h Fig. 5. Note that, since the quadrature coefficient ij is normalized with the area or volume of the cell (see Eq. 21), its value is not affected by the scaling of the distances. The coordinates of a node j in the new reference system are given as Since the grid is uniform, x j and y j are integer numbers.
As shown in Fig. 5b, the grid spacing in the new reference system is equal to 1, so that the area or the volume of each cell are A = 1 and V = 1 , respectively. Therefore, ij is simply computed as the area or volume of the intersection between the neighborhood of node i and the cell of node j: ij =Ã∕A =Ã or ij =Ṽ∕V =Ṽ . Furthermore, the only parameter which can change the values of the quadrature coefficients is the m-ratio, given as If m is unique for the whole peridynamic body (as it is often the case), then the values of ij can be computed only once and used for the neighborhoods of all the nodes.

Cell-neighborhood symmetries
The symmetries of the neighborhood with respect to the nodal grid, named "cell-neighborhood symmetries", can be exploited to reduce the number of cases to be considered. In [13] the symmetries with respect to the axes were used to compute the quadrature weights only in the first quadrant. On the other hand, we use four lines of symmetry for a 2D neighborhood (both the axes and the bisectors of the quadrants), as shown in Fig. 6. Therefore, the computation of the partial area can be carried out only for the nodes satisfying the following conditions: M ≥ y j ≥ x j ≥ 0 , where M ∶= ⌊m + 0.5⌋ ( ⌊⋅⌋ stands for the floor function and finds the greatest integer smaller than or equal to the input), and y j ≠ 0 . The latter condition is given by the fact that the central node does not interact with itself. On the other hand, the value M is used to provide an upper limit to the search for possible nodes inside the neighborhood. These nodes are enclosed by a red line in Fig. 6. Thanks to the cell-neighborhood symmetries, ij of the other nodes have the same values. For instance, the computation of the partial area for the node (x j = 2, y j = 3) is the same for nodes (3,2), Analogously, we exploit six planes of symmetry for a neighborhood in a 3D model (planes containing two axes or one axis and one bisector of the octants). Therefore, in this case, we consider only nodes that satisfy the following conditions: M ≥ z j ≥ y j ≥ x j ≥ 0 and z j ≠ 0 . The nodes considered in the proposed algorithm for the computation of the partial volumes for m = 3.2 are represented in Fig. 7.
Cell-neighborhood symmetries come into play also during the computation of some partial areas and volumes, as shown in Fig. 8. In 2D problems, the intersections between the neighborhood and cells of nodes with x j = 0 (see Fig. 8a) are symmetric with respect to the y-axis. Similarly, in 3D problems, the intersections between the neighborhood and cells of nodes with x j = 0 and y j ≠ 0 or x j = y j = 0 (see Fig. 8b, c) are symmetric with respect to the planes perpendicular, respectively, to the x-axis or both the x -and y -axis. These symmetries will be exploited in Appendix A and Sect. 3.5.

Computation of partial areas
The first algorithm proposing the analytical computation of partial areas can be found in [13]: the area of intersection between the neighborhood and the cells is subdivided into domains of basic geometry (triangles, rectangles and circular segments), for which the analytical computation of the area is straightforward. Reference [13] proposes 8 different cases of cell-neighborhood intersection. In Appendix A, we propose a novel approach to compute analytically the partial areas, which is based on the definition of the quadrature coefficient in an integral form. The integrals are then solved by distinguishing only 5 possible cases of cell-neighborhood intersections, for each of which an explicit analytical expression for the value of the quadrature coefficient is obtained. Algorithm 2 shows how to compute analytically the quadrature coefficients in 2D with the proposed approach. Function 1 is used to solve the only non-trivial integral derived from the computation of the quadrature coefficient. Please refer to Appendix A for the details of the analytical derivation, which is also very useful to better understand the extension of the formulae from 2D to 3D for the computation of the partial volumes shown in the next section. x y z Analytical computation of quadrature coefficients in 2D.
if y j = 0 and y j > M then 6: return β ij = 0 for exceeding the lower and upper limits of y j 7: end if x + y 2 2 ≤ m 2 then case-2 22: x e = m 2 − y 2 2 23:

Computation of partial volumes
We show hereinafter how to compute quasi-analytically the quadrature weights in 3D problems. The partial volumes are computed with the same approach explained in Appendix A for the partial areas. Therefore, the quadrature coefficients are computed as Possible cases of intersections between neighborhood and cell in 3D. Symmetric cell-neighborhood intersections are not shown here, but the unsymmetric portions of those intersections belong to one of the shown cases. Since the true origin of the reference system lies outside the images, it has been translated to one of the vertices of the cube for visualization clarity Figure 9 shows all the possible cases of intersection between the spherical neighborhood and a cubic cell. Note that in Fig. 9 only cells with no symmetries are illustrated, but the portions of the symmetric cell-neighborhood intersections that are used in the computation of the quadrature weights, i.e., the portions in the first octant of Fig. 8b, c, belong to one of those cases. The upper limit of the integral in x direction of Eq. (24) is the greatest x coordinate of the cell-neighborhood intersection, and it can be computed as Therefore, b x is equal to x 2 from case-2 to case-8, and to The computation of Eq. (26) at the beginning of the algorithm allows to compute case-8 and case-9 with the same formulae.
Similarly to Eq. (A1) in Appendix A, Eq. (24) can be solved for each case by splitting the integrals in correspondence of the intersections between the boundary of the neighborhood and the edges parallel to the x-axis and the faces parallel to the x-y plane (for details refer to Appendix B). There are 3 types of non-trivial integrals that derive from the previous step: where k 1 can be equal to a y , y 2 , z 1 or z 2 , k 2 to a y or y 2 , and k 3 to z 1 or z 2 . The parameters k 1 , k 2 and k 3 are defined in these ways to group the same types of integrals derived from Eq. (24). The explicit solution given in Eq. (27) is used in Function 2 to compute the integral in a general interval [a, b].
On the other hand, integrals in Eqs. (28) and (29) do not have an explicit solution. Therefore, we perform a Taylor series expansion centered at x 0 of the following functions: where c is a coefficient depending on the order n of the corresponding derivative and the indices p and q. For more details about the computation of the derivatives of a general order n (and the corresponding coefficients c(n, p, q)), please refer to Appendix C. The matrix containing all the coefficients c is obtained with Algorithm 3. If the order N of truncation of the Taylor series expansion tends to infinity ( N → ∞ ), then the solution of the integral is exact. Clearly, in a numerical algorithm N must be a finite number, but we will show that our approach is able to attain accurate results with little computational effort (i.e., with N relatively low).
Hence, we substitute Eqs. (30) and (31), respectively, in Eqs. (28) and (29) and solve the indefinite integrals: The integrals of Eq. (24) can be solved for each case by following the procedure described above (for details refer to Appendix B). Algorithm 4 illustrates how to compute the quadrature coefficients with this procedure. At the beginning of the algorithm, the coordinates of node j are considered only in their absolute value and sorted in ascending order to comply with the conditions imposed for symmetry: 0 ≤ x j ≤ y j ≤ z j . The quadrature coefficients ij computed for m = 3, 4, 6 are reported in Appendix D. 1: function Int1(a,b,k 1 ,m) : : : : :     Input: x j , y j , z j , m, N , c Output: β ij 1: x j = abs(x j ) absolute value of x j for symmetry (x j ≥ 0) 2: y j = abs(y j ) absolute value of y j for symmetry (y j ≥ 0) 3: z j = abs(z j ) absolute value of z j for symmetry (z j ≥ 0) 4: x j , y j , z j = sort(x j , y j , z j ) sort in ascending order for symmetry 5  x e 1 = m 2 − y 2 2 − z 2 2 28: x + a 2 y + z 2 2 ≤ m 2 and x 2 2 + y 2 2 + z 2 1 ≤ m 2 then case-4 32: x e 2 = m 2 − a 2 y − z 2 2 33: x + a 2 y + z 2 2 ≤ m 2 then case-5 35: x e 2 = m 2 − a 2 y − z 2 2 36: x + y 2 2 + z 2 1 ≤ m 2 then case-7 41: x e 3 = m 2 − y 2 2 − z 2 1 42: x + a 2 y + z 2 1 ≤ m 2 then case-8 or case-9 44:

Numerical results
To assess the accuracy of the computation of the partial volumes with respect to the order N of truncation of the Taylor series expansions, we compute the relative error on the total volume of the spherical neighborhood: The sum of all the quadrature coefficients ij and the volume V i = 1 of node i is the numerical value of the sphere volume, whereas (4∕3) m 3 is its analytical value. As shown in Fig. 10, the accuracy in the computation of the partial volumes is improved by increasing the value of N. Furthermore, the proposed algorithm reaches values of the errors very close to machine precision with m ≥ 3 and N = 20.
We show hereinafter several numerical results that confirm the improvements in the peridynamic integration in 3D problems. In particular, we compare the numerical results of Algorithm 1, arguably the most commonly used, with those of the new algorithm (Algorithm 4). We use N = 4 as the order of truncation for the Taylor series in Algorithm 4 since this value assures accurate results with low computational cost.

Geometrical quantities
To assess the performance of the proposed algorithm, we compute two geometrical quantities that reflect the accuracy in the computation of the partial volumes. The first one is the neighborhood volume, computed as: The analytical value of the neighborhood volume is equal to the volume of a sphere: V an = (4∕3) m 3 . The relative error on the neighborhood volume can be computed as in Eq. (34). The improvements obtained by the novel algorithm are evident in Fig. 11a, b. The second geometrical quantity for the comparison of the algorithms is the weighted volume m (Eq. (18)). The analytical computation of the weighted volume is carried out using spherical coordinated ( is the azimuthal angle and is the polar angle): where r = ‖ ‖ and erf (⋅) is the Gaussian error function. The relative error on the weighted volume can be computed as: (36) Table 1 Comparison of the cases in Fig. 4 of [13] with the cases in Fig. 15 derived from the new approach for the computation of the partial areas The abbreviation "symm." means that the cell-neighborhood symmetry should be used for the cases of the approach in [13] to obtain the cases in the new approach Approach in [13] New approach

Coefficients of the elasticity tensor
The mechanical properties of isotropic linearly elastic materials are described by the 4 th -order elasticity tensor . As in [13], we investigate the convergence behavior of the coefficients of the elasticity tensor. In state-based peridynamics, the elasticity tensor of a point in the bulk is given as [32] where is the modulus state which operates on two bonds, = � − and = �� − . The modulus state in a linear peridynamic solid model [2,32] is derived as where is the Dirac delta function defined as A component of the tensor is therefore computed as [32] (42) C an (43)  Figure 12a, c show the results of these computation, respectively, for the components C 1111 and C 1122 . The relative errors on these coefficients are computed as and are shown in Fig. 12b, d. It is evident that the proposed algorithm provides, on average, smaller errors. Furthermore, the oscillatory behavior as m increases is much smoother than that of Algorithm 1.

Manufactured problem
A "manufactured" problem, which consists in determining the force density distribution in a body under a prescribed displacement field, is solved analytically and numerically.
Following what was shown in [14] for a 2D problem, a body is subjected to the displacements where there are 18 independent coefficients for the quadratic terms. The force density derived from this displacement field is given as [14] where K = 555.56 MPa and = 416.67 MPa . Note that the force density vector an is constant for all the points in the bulk of the material, i.e., in all points that have a distance from the boundary of the body greater than or equal to 2 .
Since c u5 , c v6 and c w4 do not contribute to an , these coefficients are considered to be equal to 0. The values of the other coefficients are chosen randomly as follows: We consider only one node in the bulk of a body at which the force density vector is computed with Eq. (20). The relative errors on the components of are given as Figure 13 shows the results of the numerical computation. Algorithm 4 allows to obtain, on average, smaller relative errors and a smoother convergence behavior also in this case.

Conclusions
The peridynamic theory is a nonlocal reformulation of classical continuum mechanics based on integrals over the neighborhoods of nodes. Therefore, the numerical integration of peridynamic equations determines to a great extent the accuracy of the results. In particular, the quadrature weights, i.e., the partial areas in 2D (intersections between the circular neighborhood and the square cells) and the partial volumes in 3D (intersections between the spherical neighborhood and the cubic cells), should be computed accurately. We developed an innovative algorithm able to compute quasi-analytically the partial volumes. We use the expression "quasi-analytical" because a truncated Taylor series expansion of some functions, whose integrals are not explicitly solvable, is performed to carry out the integration. However, the new algorithm computes accurately the quadrature weights with very little computational effort, i.e., with a relatively low order of truncation of the Taylor series. The computational time required by the proposed algorithm is negligible compared to the time required to compute the bond forces and the peridynamic equilibrium equation. In particular, if m is constant in the whole domain, the same quadrature weights can be computed only once and used for every neighborhood in the body. A similar approach is also used to simplify the algorithm for the analytical computation of the partial areas, which was already developed in the literature.
Several examples have been presented to show the capabilities of the newly proposed algorithm. The numerical values of the geometrical quantities (volume of the neighborhood and weighted volume), of the coefficients of the elasticity tensor and of the manufactured problem obtained with the new algorithm are compared with those obtained with the most commonly used algorithm for the computation of the partial volumes. As the m-ratio increases, the proposed algorithm provides, on average, smaller errors and a smoother convergence behavior. For the sake of convenience, the quadrature coefficients for m = 3, 4, 6 are reported in Appendix D.

Appendix A: Analytical computation of partial areas
We propose a new algorithm for the computation of the quadrature coefficients in 2D (Algorithm 2) that solves the following integral: where x 1 = x j − 0.5 , x 2 = x j + 0.5 , y 1 = y j − 0.5 and y 2 = y j + 0.5 are the coordinates of the sides of the square cell and s is the number of symmetries of the cell with respect to the neighborhood ( s = 1 if x j = 0 and s = 0 if x j ≠ 0 , see Fig. 8a). The upper and lower limits of the outer integral, respectively, b x and a x , are scalar values that can be computed at the beginning of the algorithm. The value of the lower limit a x depends on whether the cell-neighborhood intersection is symmetric, as shown in Fig. 8a, or not. In the former case, the lower limit of the integration domain in x direction is set to a x = 0 and the value of the integral is multiplied by 2 (given that s = 1 ). In the latter case, the lower limit is the smallest x coordinate of the cell, i.e., a x = x 1 . These conditions can be written as Since only cells with y j ≥ x j ≥ 0 are considered for symmetry reasons (see Fig. 14 (a x , y 1 ) , the closest to node i, lies outside the neighborhood, i.e., a 2 x + y 2 1 > m 2 , as shown in Fig. 15e.
The comparison of these cases with the cases of the approach in [13] is summarized in Table 1.
Case-1 and case-5 are trivial since ij = 1 and ij = 0 , respectively. Furthermore, the quadrature coefficient in case-3 and case-4 can be computed with the same formulae if the upper limit of the integral in x direction is defined as The value of b x is indeed equal to x 2 in case-3, as shown in Fig. 15c, and to the x coordinate of the intersection between the boundary of the neighborhood and the lower side of the cell, as shown in Fig. 15d.
For later use, the following indefinite integral is analytically solved as This solution will be used to compute the definite integrals in different intervals. Now, we show how to compute the quadrature coefficient if the considered cell belongs to case-2. Equation (A1) can be rewritten as In case-2 (see Fig. 15b), there is an intersection between the boundary of the neighborhood and the upper side of the cell, which has an x coordinate equal to The integral in Eq. (A5) must be splitted because the integrand has different values in the intervals [a x , x e ] and [x e , b x ] . Then, the integral can be solved analytically by using Eq. (A4): Similarly, the quadrature coefficient for a cell belonging to case-3 or case-4 is analytically computed from Eq. (A5) as follows: The computation of the partial area of a cell is summarized in Algorithm 2. Function 1 is used to compute the integral of Eq. (A4) in a general interval between a and b. The proposed algorithm yields exactly the same results of the algorithm for the computation of the partial areas presented in [13], but its implementation is simpler because it is based on a smaller number of cases.

Appendix B: Quasi-analytical computation of partial volumes
We can rewrite Eq. (24) that describes the integral to compute the quadrature coefficients as follows: w h e r e x 1 = x j − 0.5 , x 2 = x j + 0.5 , y 1 = y j − 0.5 , y 2 = y j + 0.5 , z 1 = z j − 0.5 and z 2 = z j + 0.5 are the coordinates of the faces of the cubic cell and s is the number of symmetries of the cell with respect to the neighborhood (see Fig. 8b, c). a x , a y and b x are scalar values computed as shown in Eqs. (25) and (26).
Since the integrand and the upper limit of the inner integral in Eq. (B1) are discontinuous functions, the integral domain should be split into subdomains in which only continuous functions are integrated. The integral is split differently for each of the 10 cases represented in Fig. 9. The cases are differentiated depending on the intersections between the boundary of the neighborhood and the edges of the cell parallel to the x-axis and the faces parallel to the x-y plane. These intersections are indeed used to split the integral of Eq. (B1), as shown in next subsections. To begin with, we consider a cell completely within the neighborhood (case-1 in Fig. 9a) and decrease gradually the m-ratio to find other cases. We also need to check that, for the considered cells with z j ≥ y j ≥ x j ≥ 0 , other cases of intersections are not possible.

Case-1
In case-1, the cell is completely within the neighborhood, as shown in Fig. 9a. The condition for which this case is verified is that the farthest node with coordinates (x 2 , y 2 , z 2 ) is inside the neighborhood: In this case, the solution of Eq. (B1) is trivial:

Case-2
In case-2, the boundary of the neighborhood intersects the edge of the cell between the vertices (a x , y 2 , z 2 ) and (x 2 , y 2 , z 2 ) , as shown in Fig. 9b. Therefore, the condition in this case is To be sure that there are no other possible intersections with edges parallel to the x-axis and faces parallel to the x -y plane, we need to check that the vertices (x 2 , a y , z 2 ) and (x 2 , y 2 , z 1 ) lie within the neighborhood, so that the edge between (a x , a y , z 2 ) and (x 2 , a y , z 2 ) and the face on the plane z = z 1 are, respectively, within the neighborhood. The condition for the vertex (x 2 , a y , z 2 ) being within the neighborhood is given as We rewrite the condition in Eq. (B5) as By comparing this inequality with the condition on the lefthand side in Eq. (B4), we obtain: If the cell has two symmetries as shown in Fig. 8c ( s = 2 , a x = a y = 0 ), Eq. (B7) becomes y 2 2 ≥ x 2 2 which is always verified. If the cell has one symmetry as shown in Fig. 8b  ( s = 1 , a x = 0 , x 2 = 0.5 , y j ≥ 1 ), Eq. (B7) becomes 2y j ≥ 0.5 2 which is always verified. If the cell has no symmetries ( s = 0 ), Eq. (B7) becomes 2y j ≥ 2x j which is always verified.
Similarly, one can show that also the vertex (x 2 , y 2 , z 1 ) lies within the neighborhood. In this case (see Fig. 9b), the outer integral of Eq. (B1) should be split at the intersection of the boundary of the neighborhood with the edge between the vertices (a x , y 2 , z 2 ) and (x 2 , y 2 , z 2 ) , namely at Furthermore, the inner integral should be split in correspondence of the arc of circle with equation which is given by the intesection of the boundary of the neighborhood with the plane z = z 2 (where the face of the cell lies). Therefore, Eq. (B1) becomes: Note that the upper limit of the inner integral is y 2 because the face of the cell on the plane z = z 1 is completely within the neighborhood. Since the solution of the inner integral of the last term in Eq. (B10) is similar to the integral solved in Eq. (27), the quadrature coefficients in case-2 can be computed as (B9) The integrals in Eq. (B11) are of the types solved in Eqs. (27), (32) and (33). To shorten the formulae, we define the following functions (similar to Functions 2, 3 and 4 used in Algorithm 4): where f (x) and g(x) are functions defined in Appendix C and their corresponding n-th derivatives f (n) (x) and g (n) (x) are computed as in Eqs. (C2) and (C7). The parameter k 1 can be equal to a y , y 2 , z 1 or z 2 , k 2 to a y or y 2 , and k 3 to z 1 or z 2 depending on the cases. Using the previous definitions, Eq. (B11) is solved as follows:

Case-3
In case-3, the boundary of the neighborhood intersects the face of the cell on the plane z = z 2 , but does not intersect any edge parallel to the x-axis, as shown in Fig. 9c. Therefore, the condition in this case is The condition on the right-hand side implies that there are no intersections of the boundary of the neighborhood with the edge of the cell between the vertices (a x , y 2 , z 2 ) and (x 2 , y 2 , z 2 ) . Similarly, the condition on the left-hand side implies that there are no intersections with the edge between the vertices (a x , a y , z 2 ) and (x 2 , a y , z 2 ) . With the same approach illustrated at the beginning of Sect. B.2, one can also show that the vertex (x 2 , y 2 , z 1 ) lies within the neighborhood for the considered cells ( z j ≥ y j ≥ x j ≥ 0 ), so that the face of the cell on the plane z = z 1 is completely inside the neighborhood.

Case-4
In case-4, the boundary of the neighborhood intersects the edge of the cell between the vertices (a x , a y , z 2 ) and (x 2 , a y , z 2 ) and the face on the plane z = z 1 lies completely within the neighborhood, as shown in Fig. 9d. Therefore, the conditions in this case are The intersection of the boundary of the neighborhood with the edge between the vertices (a x , y 2 , z 2 ) and (x 2 , y 2 , z 2 ) is given as Given the conditions in previous case, there are no intersections of the boundary of the neighborhood with the edge of the cell between the vertices (a x , y 2 , z 2 ) and (x 2 , y 2 , z 2 ).

Case-5
In case-5, the boundary of the neighborhood intersects both the edge of the cell between the vertices (a x , a y , z 2 ) and (x 2 , a y , z 2 ) and the face on the plane z = z 1 , as shown in Fig. 9e. Therefore, the conditions in this case are With the same approach illustrated at the beginning of Sect. B.2, one can also show that the vertices (a x , y 2 , z 1 ) and (x 2 , a y , z 1 ) lie within the neighborhood for the considered cells ( z j ≥ y j ≥ x j ≥ 0 ), so that there is an intersection with the edge between (a x , y 2 , z 1 ) and (x 2 , y 2 , z 1 ) and the edge between (a x , a y , z 1 ) and (x 2 , a y , z 1 ) lies completely within the neighborhood. The intersection of the boundary of the neighborhood with the edge between the vertices (a x , y 2 , z 2 ) and (x 2 , y 2 , z 2 ) is given in Eq. (B20), whereas the intersection with the edge between (a x , y 2 , z 1 ) and (x 2 , y 2 , z 1 ) is given as Note that x e 2 ≤ x e 3 since: If the cell has two symmetries ( s = 2 , a y = 0 ), Eq. (B25) becomes 2z j ≥ 0.5 2 which is always verified ( z j ≥ 1 ). If the cell has one or no symmetries, Eq. (B25) becomes 2z j ≥ 2y j which is always verified.