Implicit implementation of the nonlocal operator method: an open source code

In this paper, we present an open-source code for the first-order and higher-order nonlocal operator method (NOM) including a detailed description of the implementation. The NOM is based on so-called support, dual-support, nonlocal operators, and an operate energy functional ensuring stability. The nonlocal operator is a generalization of the conventional differential operators. Combined with the method of weighed residuals and variational principles, NOM establishes the residual and tangent stiffness matrix of operate energy functional through some simple matrix without the need of shape functions as in other classical computational methods such as FEM. NOM only requires the definition of the energy drastically simplifying its implementation. The implementation in this paper is focused on linear elastic solids for sake of conciseness through the NOM can handle more complex nonlinear problems. The NOM can be very flexible and efficient to solve partial differential equations (PDEs), it’s also quite easy for readers to use the NOM and extend it to solve other complicated physical phenomena described by one or a set of PDEs. Finally, we present some classical benchmark problems including the classical cantilever beam and plate-with-a-hole problem, and we also make an extension of this method to solve complicated problems including phase-field fracture modeling and gradient elasticity material.


Introduction
The nonlocal theory of elasticity [1][2][3][4][5] primarily considers distant action forces between objects. Different from the concept of local theory, an object can interact without physical interaction with a different object. The theory is of great significance to solve many physical problems, such as the law of universal gravitation. Some notable nonlocal numerical methods were proposed according to nonlocal interaction in nonlocal continuum field theories.
Peridynamics (PD) [6][7][8] is a formulation of continuum mechanics on the basis of the concept of nonlocal integration, PD avoids the singularity of the traditional local differential equations when solving discontinuous problems.
One key application of PD is fracture. It shares the same advantages as the Cracking Particles Method (CPM) presented in [9][10][11]. In contrast to many other discrete crack approaches as presented in [12][13][14][15][16][17][18][19][20], PD does not require the representation of the discrete crack surface and associated crack tracking algorithms. PD also has been successfully applied to rock fracture and soil damage analysis, such as impact fracture [21,22], composite material separation [23,24]and beam and plate structures [25]. However, to eliminate erroneous wave reflection and ghost force among particles, all traditional PD formulas must use the same horizon size. In many applications, to enhance calculating performance, its necessary to use different horizon sizes for the calculation of particles with non-uniform spatial distribution, such as adaptive encryption, multi-scale simulation, and multi-body analysis. In other words, in order to balance the calculation efficiency and calculation accuracy, we hope that PD can be based on the distribution characteristics of the particles, but if the size of the near field is used, it will result in the generation of false stress waves and the problem cannot be solved correctly. To address the aforementioned 1 3 issue, dual-horizon PD [26,27] was developed to improve computing efficiency and to allow for varying horizon sizes. The dual-horizon is the dual term of the horizon when variable horizons are used in the inhomogeneous discretization. It separates the horizon that exerts forces and counter forces between the particles, thereby solving the problem of false stress caused by the horizons. In addition, other nonlocal models mainly include nonlocal linear elasticity [28][29][30], dynamics of nonlocal fluid [31,32], electromagnetic nonlocal theory [33,34], nonlocal damage model [35][36][37] and nonlocal calculus [38][39][40][41].
In most of the current nonlocal numerical methods, the governing equation is generally solved by the explicit time integration method [25,42]. The explicit algorithm [43][44][45][46][47] is based on dynamic equations, without iteration and has good stability. However, the mass matrix is required to be a diagonal matrix in the explicit solution and the speed advantage can be exerted only when the unit-level calculation is as small as possible. Therefore, the reduced integration method is often used, which is easy to excite the hourglass mode and affect the stress and strain calculation accuracy. To ensure its conditional convergence, the integration step must be less than the minimum value of the free vibration period of all elements. The integration time step is generally 1/1000 to 1/100 of the implicit step, which is suitable for solving instant or short-time loading problems. For the quasi-static (quasi-static) problem, when the total load is fixed, the calculation takes a long time, which seriously affects the calculation efficiency. Compared with the explicit method, implicit method [48][49][50] has advantages for fast convergence and lower computational cost.
In recent years, several numerical approaches based on peridynamics operators have been proposed, see e.g. the contributions in [51][52][53]. A new computational method based on nonlocal operators is the NOM first proposed in [54] for electromagnetic problems. The approach has been later on extended to mechanical problems in [55,56]. The NOM can be considered as a generalization of non-ordinary state-based PD. It has been applied to numerous challenging problems in solid mechanics and can be a viable alternative to FEM or meshless methods. In order to acquire the differential operators, FEM and meshless methods are required to establish the shape functions as well as compute their derivatives, however, NOM can acquire the differential operators easily without the use of shape functions. The tangent stiffness is obtained naturally by simply defining an energy function thus drastically simplifying its implementation. In combination with the weighted residual method and variational principle, the residual and the tangent stiffness matrix can be established by NOM with ease. NOM is enhanced here also with operator energy functional to achieve the linear consistency of the field and avoid instabilities. Although the theoretical framework of the NOM has been proposed, other benefits of this numerical method have not been thoroughly discussed. Furthermore, the details of the derivation of the first-order and higher-order nonlocal differential operators, construction of the first-order and higher-order operator energy functional, the derivation of the residual and tangent stiffness matrix of operate energy functional and the detailed implementation procedure of first-order and higher-order NOM has not been shown before. The purpose of this paper is to describe in detail the method of implicitly implementing for first-order and higher-order NOM, which mainly including the derivation of the first-order and higher-order nonlocal differential operates, the detailed form of firstorder and higher-order tangent stiffness matrix for operate energy functional and elastic material constitutions in different conditions. The Mathematica code of first-order and higher-order NOM is presented and explained in detail, and it will be an effective tool for studying complicated physical problems.
The remainder of the paper is outlined as follows: In Section 2, we briefly reviewed the NOM and elaborated on the fundamental concept of support and dual-support. In Section 3, we derived the first and higher-order implicit nonlocal differential operators based on the Taylor series expansion. Then the elastic material constitutions are also derived. Hereafter, to remove the zero-energy mode, the first and higher-order operator energy functional by the nonlocal operator is constructed, combined with the method of weighed residuals and variational principles, residual and tangent stiffness matrix of operate energy functional are established. Finally, we present the detailed steps of the first-order and higher-order NOM implementation process. To demonstrate the capabilities of NOM, four numerical examples are presented in section 4. Finally, we conclude in section 5.
ij is the spatial vector established in S j . Figure 1b illustrates the concept of support and dual-support.
In calculus calculations, the NOM replaces the local operator with the fundamental nonlocal operators. By substituting the local differential operator with the corresponding nonlocal operator, the functional defined by the local differential operator can be utilized to establish the residual or tangent stiffness matrix.
The nonlocal gradient operator for a vector field and scalar field u at point i in support S i are defined as [55] where w( ij ) is the weight function for vector ij in support S i .
By nodal integration, the nonlocal gradient operator and its variation for a vector field in discrete form can be written as The operator energy functional for vector field at a point i is defined as (2) is a coefficient for the operator energy functional, is the normalization coefficient, i is a shape tensor and it can be defined as 3 Implementation

Derivation of the first-order implicit nonlocal differential operators
The displacement field at any point in the elastic body can be represented by three displacement components u, v, w along the rectangular coordinate axis, and it's three dimension vector form is The first-order Taylor series expansion at point (0,0) for a displacement field is given as ; ∶= (x, y, z) T denotes the initial bond vector, O( 2 ) represents higher-order terms, and for linear field O( 2 ) = 0. � − ⊗ can be obtained from Eq. 9 and can be shown as The gradient operate ∇ can be expressed as

3
The discrete form of Eq. 11 at point i can be shown as To simplify the equation more conveniently, letting The nonlocal gradient operator ∇ at point i can be rewritten as According to Eq. 13, the matrix form of nonlocal gradient operator ∇ at point i for vector field can be expressed as For the convenience of calculations, we transform the matrix form of nonlocal gradient operator ∇ i for vector field into vector form and it can be rewritten as Similarly, the nonlocal gradient operator ∇ u at point i for scalar field is written as The 2D form of the first-order nonlocal gradient operators is given in Appendix A.

Derivation of the higher-order implicit nonlocal differential operators
The higher-order NOM is based on higher-order Taylor series expansion of a multi-variable function. Consider a vector field at point j ( j ∈ S i ). For convenience, we shorthand ( j ) by j , which can be estimated via a Taylor expansion based on i in r dimensions with the highest order of derivatives as n: where where r represents the dimension and n denotes the highest order of derivatives. n r is a compilation of multi-indexes that have been flattened and it can be written as To eliminate the round-off error in Eq. 17, we here include the characteristic length scale l i of support S i at i , and the modified Taylor series expansion may be rewritten as (17) F o r a n y m u l t i -i n d e x (n 1 , … , n r ) ∈ n r , l i,n 1 …n r = l n 1 +⋯+n r i n 1 !…n r ! i,n 1 …n r , ∀(n 1 , … , n r ) ∈ n r . We let p l j , l i and i denotes the flattened polynomials, scaled partial derivatives, partial derivatives, respectively, according to multi-index notation n r and they can be shown as The l in the Eq. 23 allows the terms of the same characteristic scale for length thereby eliminating computational roundoff error. The current partial derivatives can be recovered by In which diag [ X 1 , … , X n ] signifies a diagonal matrix of X 1 , … , X n , whose diagonal entries begin from the upper left corner. Hence, the expansion of the Taylor series with i to the left of the Eq. 22 can be rewritten as where Therefore, the nonlocal operator ̃ i and its variation can be obtained as  The nonlocal operator gives all partial derivatives with the highest order up to n . In PDEs, the group of derivatives is a subgroup of the nonlocal operator. Each term in ̃ i corresponds to a row of i p l wi multiplied by Δ i . Equation 31 can be used to substitute the differential operators in PDEs to generate a strong form of algebraic equations. Meanwhile, we can solve the linear (nonlinear) weak formulations using the weighted residual methodology and the variational principle, hence the variation of i in Eq. 32 is required in these circumstances. Equation 31 can also be shown more succinctly as where i is the nonlocal operator coefficient matrix of point nodal values in support, the operator matrix generates all partial derivatives of maximal order smaller than | | + 1.

Elastic material constitution
For elastic material, the strain energy functional is a function of the deformation gradient F . According to the principles of traditional solid mechanics, the deformation gradient F in 3D form expressed as where 3×3 denotes the identity tensor. The first Piola-Kirchhoff stress can be derived from the directly derivative of strain energy (F) in the context of total Lagrangian formulation, and it can be derived as In addition, the 4th order elastic tensor 4 can be obtained using derivation of the first Piola-Kirchhoff stress To obtain the elasticity matrix, based on the Voigt notation, 4 can be written as a matrix form D 9×9 w h e r e (= 11 , 12 , ⋯ , 33 donates the flattened deformation gradient and flattened first Piola-Kirchhoff stress. In particular, for isotropic linear elastic material, regarding infinitesimal deformations of a continuous linear elastic material with a tiny displacement gradient relative to unity ( i.e. ∇ ≪ 1 ), any of the strain tensors utilized in finite strain theory (such as the Lagrangian strain tensor) may be geometrically linearized. The non-linear or second-order elements of the finite strain tensor are ignored in this linearization. Hence, we can obtain the Lagrangian strain tensor = 1 2 (F + F T ) − and stress tensor = ∶ .
The internal functional energy for 3D linear elastic material can be expressed as Likewise, the internal energy functional for plane stress and plane strain conditions can be shown as follows where For linear elastic material, the elastic matrix D in plane stress, plane strain and 3D conditions can be expressed as where , represent the Lamé constants, which are related to the Young's modulus E and Poisson's ratio :

Construction of the first-order and higher-order operator energy functional
As a particle-based method, when using node integration the first-order and higher-order NOM suffer from a zero-energy mode [57,58], which results in numerical instability. To eliminate the effect of the zero-energy mode, traditional PD and SPH introduce a penalty term to the force state [51]. Nevertheless, the approach described above is only applicable in the explicit time integration formulation. NOM employs operator energy functional for nonlocal gradients w h e r e ∇u i = [ It should be noted that the shape tensor i is involved in ∇ ⊗ i ∶ ∇ ⊗ i ⋅ i and the operator energy functional is valid in any dimension.
For convenience, we let To facilitate numerical implementation, we transform Eq. 47 into a detailed form and the point i displacement vector and first-order differential of displacement vector in 3D can be shown as � to achieve the linear field of the field and avoid numerical instabilities.
In first-order NOM, the operator energy functional at point i can be construct according to the first-order nonlocal operator ∇ ⊗ i , it can be expressed in discretization form as In this paper, the special variation ̄F , ̄2 F and F , 2 F are defined as Hence, the following relationships can be derived as According to Eqs. 53 and 55, Eq. 47 can be rewritten as The global tangent stiffness matrix of operate energy functional, internal residual and tangent stiffness matrix of physical energy functional in support S i can be obtained by Finally, the summation of 3D form for the first-order global tangent stiffness matrix and hourglass tangent stiffness matrix in support S i can be obtained The derivation of the first-order hourglass tangent stiffness matrix and the global tangent stiffness matrix in 2D form is given in Appendix B.
The operator energy functional for higher-order NOM at point i can be construct according to the higher-order nonlocal operator l i , it can be expressed as By submitting the Eq. 34 into the below Eq. 62 and it can be rewritten as The first and second derivative F hg i yield the higher-order residual and the tangent stiffness matrix of operator energy functional and can be expressed as The higher-order global tangent stiffness matrix of operator energy functional in support S i can be computed by Finally, the summation of the higher-order global tangent stiffness matrix and hourglass tangent stiffness matrix in support S i can be obtained

Numerical implementation with an open-source code
The numerical implementation of the first-order and higher-order NOM are summarized as the following steps, and a flow chart of the NOM implementation procedure is depicted in Fig. 2.

Step 1. Discretization of the solution domain.
Consider the solution domain to be a discrete domain consisting of discrete points of varying sizes and shapes (67) that are linked to one another. There are two methods to achieve the discretization of the solution domain. For the first method, a user-defined subroutine GridDomain is customized, which can discretize the solution domain evenly. This method is mainly used for the discretization of the rule solution domain. The second method can be achieved as follows: initially, the model is created using finite element software (such as ABAQUS), then the model is divided into grids of different sizes according to the characteristics of the model and export the model information into a inp format file. The model mesh and node information can be read through a user-defined subroutine ParseAbaqusFile, therefore the discrete solution domain with different densities can be achieved. The Mathematica code for discretization of the solution domain is shown below. In the user-defined subroutine GridDomain, where xmin, xmax are the minimum and maximum values of solution domain and dx is the spacing between points.
Step 2. Definition of the problem and solution domain. Determine the overall attributes and geometric boundary conditions of the solution domain based on the real case, and assign values to the corresponding parameters.
Step 3. Specify and search for the number of neighbors for each point in support and build index numbers for the specified neighbors.
In this step, a user-defined subroutine NeiList is customized by the built-in command Nearest of Mathematica to find the required number of neighbors, and the process of building index number for the specified neighbor points can be achieved by the user-defined subroutine NeiIIndex. The Mathematica code for finding the specified area points and specified number neighbor points in the solution domain is shown below, where coord is the coordinate of the points, numNei represents the number of specified neighbors and ndim denotes the dimension of coordinate.
Step 4. Establish the nonlocal operator coefficient matrix, hourglass tangent stiffness matrix, and the summation of global tangent stiffness matrix and hourglass tangent stiffness matrix in first-order and higher-order form.
1 First-order nonlocal operator method. As shown in below Mathematica code, three user-defined subroutines BHgmatrix, Dmatrix and Kmatrix are customized to calculate the nonlocal operator coefficient matrix B , hourglass tangent stiffness matrix hg K , elastic material constitution matrix D and the summation of hourglass tangent stiffness matrix and tangent stiffness matrix . In the user-defined subroutine BHgmatrix, initially the basic forms B i , hg K i and shape tensor i at point i are constructed, which are named as bmatrix, hgmatrix and kmatrix, respectively. Subsequently by Eq. 7, shape tensor can be calculated; by the equation (1, 1, 1) ⊗ (1, 1, 1) T and Eq. 56, the matrix j can be calculated and the hg K i matrix can be assembled. Similarly, by equation Eq. 15, the intermediate variables j can be obtained and the B i matrix is assembled. In the user-defined subroutine Dmatrix, by Eqs. 43, 44 and 45, the linear elastic solid elastic matrix D for plane stress, plane strain and 3D conditions can be obtained. Where type=1,2,3 represent the three constitutive models mentioned above. In the user-defined subroutine Kmatrix, by Eq. 67, the summation of global tangent stiffness matrix and hourglass tangent stiffness matrix can be established. In the below Mathematica code, where WeiF represents the weight function, voli represents the volume of each point, Es, mu represents the Young's Modulus and Poisson's ratio, respectively. 2 Higher-order nonlocal operator method. The main codes for higher-order NOM are presented in Appendix C, which includes constructing the global higherorder nonlocal operator coefficient matrix , global hourglass tangent stiffness matrix hg K and the summation of global tangent stiffness matrix and hourglass tangent stiffness matrix . Initially, a user-defined subroutine Multi-Index is customized to construct multi-index notation in

Step 5. Search boundary points in discrete domain and impose corresponding boundary conditions based on actual problems.
A user-defined subroutine FindPoints is customize for searching boundary points. When using user-defined subroutine FindPoints, it's need to call another customized subroutine LessThancoord at the same time. Subsequently, to employ Dirichlet and Neumann boundary conditions to the boundary points, two user-defined subroutines Dirichlet-BoundaryApply and NeumannBoundaryApply are customized. Where the penalty method is used in the user-defined subroutine DirichletBoundaryApply. The related Mathematica codes are shown as below, where Ksp represents the summation of global tangent stiffness matrix and hourglass tangent stiffness matrix , Rsp represents the global internal residual, pd represents the index of the points of application of specified displacement, pc represents the penalty coefficient and pf represents the index of the points of application of specified force. Step 6. Solve the global displacement in the discrete domain.
According to step 1 to step 5, the summation of global tangent stiffness matrix and hourglass tangent stiffness matrix and global residual vector R can be calculated. And the global displacement matrix can be calculated by solving the following linear algebra Eq. 68.
we employ the built-in command LinearSolve of Mathematica to solve the global displacement vector, then we use the user-defined subroutine UMatrix convert displacement vector into displacement matrix form. The related Mathematica codes are shown below.

Numerical examples
Four numerical examples in 2D or 3D are presented to validate the first-order and higher-order NOM in this section. The numerical results are compared with the analytical solution or that by FEM software to verify the feasibility.

A cantilever beam under shear load
In this section, a 2D cantilever beam loaded at the end with shear load is considered. To obtain the displacement of each point, we need to solve the Eq. 68, which can be shown in detail as Eq. 72. In this work, we use the penalty approach to apply Dirichlet boundary conditions. This method can be achieved through the following steps: When the displacement of point i u i = u i , we modify the i-th equation as follows, multiply it's diagonal element K ii by a penalty factor (In computations, is set to 10 10 .) and replace R i with k ii u i to obtain: The modified i-th equation can be expressed as: Since k ii ≫ k ij (i ≠ j) , the k ii u i term at the left end of the equation is much larger than the other terms, so it can be approximated k ii u i ≈ k ii u i . Then we can obtain u i ≈ u i . According to the Eq. 73, we can obtain the displacement of each point. It should be noted that this method is suitable for any given displacement of points, as the order of the formula and the displacement order of the points remain unchanged during the solution process. The number of target point's neighbors and the radius of support size in NOM can be very flexible and unlimited. However, a large number of neighbors is expensive in the k-nearest neighbor's search. Based on our numerical experience, eight points closest to the target point are selected to construct the target point's support domain. The difference between the numerical result and analytical solution is measured by the L2-norm, which is calculated by  Fig. 4 depicts the displacement for discretization at Δx = H∕50 . The displacement cloud diagram of the cantilever beam for discretization at Δx = H∕50 using first-order NOM numerical results and analytical results can be seen in Fig. 5. As shown in Table 1, Figs. 4 and 5, good agreements can be seen between the firstorder NOM numerical results and analytical results. It shows very close variations in L2-norm and the displacement errors in different directions.

Problem of an infinite plate with a hole in tension
Consider a 2D infinite plate with a circular hole, as shown in Fig. 6. The plate's length is L and there is a small circular hole with a radius of = a in the plate. Because the plate's thickness is substantially smaller than its length, it can be considered a plane stress problem. The force boundary condition of P = 1 MPa. Based on the symmetry of the structure and load, one-quarter of the model is used for analysis. When the plate's width is substantially more than the radius of the small hole, we can obtain the stresses analytical solution [61] in polar coordinate system according to classical elasticity theory, and it can be expressed as The highest hoop normal stress is achieved at the hole's edge, as illustrated in Fig. 6. When = a , = 2 ( 3 2 ) , attains the maximum value of three times the uniformly distributed stress P , ( ) max = 3P , and sharply approaches to P as it moves away to the edge. We characterize the stress concentration phenomena using the stress concentration factor K in this case. We solve stress concentration problems by the first-order NOM, plane stress condition is considered. The numerical results were compared with ABAQUS standard and analytical solutions. The parameters for the plate include E = 3 × 10 4 MPa, = 0.3 , and the radius of the hole is = 1 m. To facilitate the application of force boundary conditions, we converted the Eq. 76 from polar coordinate system to Cartesian coordinate system, it can be shown as To obtain a good discrete result, initially, we build the model and meshes in ABAQUS, then export it as an "model. inp" file, which includes the element and point coordinate information. We read the exported "model.inp" file in the environment of Mathematica. We can calculate the coordinates of each point in the model, the area of each element, and assign parameters such as area and force to the relevant points of the element according to the core principle. So as to realize the discreteness of the unit in the environment of Mathematica. It should be noted that only the nodes are used and no interaction between elements. The Mathematica software reads the "model.inp" model (L∕a = 5) and discretization of the model as shown in Fig. 7. In this study, we use the built-in command Nearest of Mathematica to find the required number of neighbor points, and eight neighbors are selected. We fix the diameter of the hole to 2m and change the L . Four cases of relative size of plate width and hole radius with L∕a=5,7,9,11 and four cases with total 525,2050,4575,8100 nodes are investigated. The plate is discretized according to ABAQUS mesh elements and CPS3 elements are adopted in ABAQUS [62] to calculate the reference results. At first we test the maximum value of stress in the x-direction ( xx ) max , according to the formulation ( xx ) max P , we can obtain the corresponding K and compared it to the analytic solutions. According to the first-order NOM, ABAQUS standard, and analytical solution, the statistical results of the stress concentration factor K, as shown in Table 2.
The displacement and stress distribution of the model(L∕a = 5) on polar coordinate system = 2 at the case of with total 4575 and 8100 nodes using first-order NOM are compared to the analytic solutions, as seen in Figs. 8 and 9. The displacement's L2 norm computed by Eq. 75 and stress field at four different discretization cases are shown in Table 3. These characteristics agree well with the theoretical analysis of this problem mentioned above, and it proved that the first-order NOM established in this paper has good applicability.

Crack propagation problem by NOM with phase-field modeling
Phase field fracture modeling [63][64][65][66][67][68] shows great advantages for modeling complex crack propagation, and its based on the Γ-convergence of the approximations of free discontinuity problems [69]. In this section, we use NOM to establish implicit phase field modeling and a phase field (x, t) ∈ [0, 1] is introduced to approximate the fracture surface, Γ . Where (x, t) = 1 denotes the complete damage of the material and (x, t) = 0 denotes the material is undamaged.
According to literature [70], the density of the crack surface per unit volume ( ) and the surface energy in the solid due to the formation of crack can be shown as   where G c represent critical energy release rate, l 0 is a length scale factor that governs the phase field's transition area, hence it can reflect the width of the crack, When l 0 → 0 , the Γ-limit recovers the sharp discontinuous interface of the crack.
To guarantee that the crack is only driven by tensile load in phase field simulation, the elastic energy must be decomposed into tensile and compressive components [71]. The strain tensor can be decomposed as follows by using Eigen-decomposition where the tensile and compressive parts of strain tensors are + and − , respectively; ⟨x⟩ ± ∶= (x ± �x�)∕2 ; the principal strain is a , and the principal strain direction is a .
The elastic energy density ( ) and Cauchy stress tensor for the tensile and compressive parts are represented by the decomposed strain tensor as following where 0 ( 0 > 0 and 0 ≪ 1 ) is a tiny positive factor that prevents the positive component of the elastic energy density from degrading completely, , are the Lamé parameters. I donates identity tensor.
The overall potential functional is calculated as the sum of the phase field approximations for the fracture energy, elastic energy, body energy, and potential energy caused by the surface load, it can be shown as where denotes the body force, u donates the displacement, and f donates the surface loads at the boundary, W int donates internal energy, W ext donates external work.
The Lagrange energy functional can be expressed as According to Hamilton's principle, the first variation of the Lagrange energy functional L should be 0, namely L = 0 , hence we can obtain the following phase-field governing equations and Neumann boundary conditions where n * is the outward-pointing normal vector for boundary Ω f . To maintain the phase field that increases monotonically [70], the local history field of strain H(x, t) is introduced and it can be shown as we use H(x, t) replace the + e in the Eq. 86, then we can obtain Since the implicit NOM is based on the governing equations weak form, and the weak forms of Eq. 88 are formulated as follows The continuous support S i is expressed in NOM by where j 1 , .., j k , .., j n denote the global indices of point i's neighbors.
The nonlocal operator for displacement ∇ ⊗ u i and phase field ∇ i in discrete form can be rewritten as where ≃ denotes discretization, ũ i (= i , j 1 , .., j k , .., j n ) is all the variations of the displacement in support S i , ̃i (= i , j 1 , .., j k , .., j n ) represents all variations of the phase field in support S i . ΔV j is the volume of neighbor point j, B u i and B i are the the coefficient matrices, as shown in Eqs. 15 and 16.
The research was based on a library of linear elastic material. According to the above expressions, Eq. 89 can be shown as We can obtain the point i residual of the overall systems are internal force, external force, respectively, which correspond to the displacement field.
Whereas F ,int i is the internal force corresponding to the phase field. The corresponding tangent stiffness matrices can be obtained based on the internal forces: where ijkl is a elasticity matrix, and it can be obtained using Eigen-decomposition algorithm for fourth-order isotropic tensor [70] The stress based on the Eigen-decomposition method can be expressed as where , are the Lamé constants, I 3×3 is the identity matrix, n 1 , n 2 , n 3 are Eigenvectors for principal strains 1 , 2 , 3 of (= ∑ 3 i=1 i n i ⊗ n i ) (Fig. 10). Then Eq. 110 can be shown as However, Eq. 89 suffers from the zero-energy mode, the penalty force from the nonlocal operator functional should be added, then the summation of global tangent stiffness matrix and hourglass tangent stiffness matrix in displacement and phase field as shown in following In this case, we consider a single notched square plate subjected to static tension loading. Figure 12 depicts the geometrical and boundary conditions of the plate. We consider the plate parameters as follows: E = 2.1 × 10 5 MPa , = 0.3 and the energy release rate g c = 2.7 × 10 −3 kN/mm . The plane strain condition is considered in this test. We deleted the points closest to the initial crack to produce the initial crack and each point has eight neighbors to construct their own support domain. The single edge notched square plate is discretized into 40401 points. The phase-field length scale is set to l 0 = 0.015 mm. Due to the maximal order of partial derivatives in Eq. 89 is 1, the first-order NOM is employed. The penalty method is used to enforce the Dirichlet boundary conditions for phase and displacement fields. In this case, to prevent the singularity during computation, 0 = 1 × 10 −6 is chosen. Due to the mutation property of crack propagation, in the staggered scheme, a small vertical displacement increment Δu = 1 × 10 −5 mm is applied to the upper boundary of the plate and the bottom boundary is fixed all the time. Figure 11 depicts the crack patterns in horizontal direction at various displacements for l 0 = 0.015 . Figure 12 depicts the reaction force-displacement curves. As expected in the literature [72], when adopted staggered scheme, the rate of the crack evolves seems to be delayed in comparison to a completely monolithic scheme. In comparison with the results by Miehe [70,72], the load curves obtained by the implicit nonlocal model align well with the reference result with the increase of the vertical displacement Δu . The crack tip still has a small bond force caused by the support and dual-support which leads to slight differences appear at post localization where FEM exhibits a sharper decay compared to the NOM.

3D gradient elasticity cantilever beam under shear load
To illustrate the feasibility of higher-order NOM, we expand it to handle the gradient elasticity beam issue. The isotropic elasticity gradient material's energy functional [73,74] can be represented as where ̄ represent the Cauchy-like stress tensor. represent the gradient material factor, = 1 2 (∇u + u∇) represent the infinitesimal strain tensor. b represent body force density, whereas ℙ and ℝ represent the traction force and double traction force act on on Ω . n * represents the outward-pointing normal vector for the boundary surface.
Based on the strain gradient elasticity theory [75], the variation of the strain energy and the external work can be written as where denotes the third rank double stress tensor.
In Eq. 103, the Cauchy-like quantities ̄ , work-conjugate to , and the double stress , work-conjugate to ∇ , are derived from the partial derivatives of the strain energy density with respect to and ∇ [74], respectively. Accordingly, we can obtain ̄ and as shown in follows where tr = ∇ ⋅ u , is the symmetric second rank unit tensor.
The Cauchy-like stress tensor ̄ and strain tensor are coupled typically by Hooke's law [76], but the double stress tensor is expressed on the basis of the Cauchy stress tensor, yielding the total stress tensor in the form According to the principle of virtual work, F int − F ext = 0 is valid for u in Ω and on surface Ω , which lead to the following governing equation and boundary conditions in terms of the Cauchy-like stress are written below To make the expression clearer, we simplify the Eq. 102 and make some modifications to following  According to the Eq. 108, we can establish the residual and the corresponding tangent stiffness matrix at point i of the overall system: i5 ΔS j A 3D gradient elasticity cantilever beam under uniform shear load is considered. The cantilever beam with dimensions of 25 × 10 × 3 m 3 . The material parameters are considered: E = 6 × 10 4 MPa, = 0.33 . The uniform shear force P y = −1.0 × 10 −3 MPa. The higher-order NOM based on numerical integration [77] is employed to investigate the effectiveness of the current method for gradient elasticity cantilever beam. The cantilever beam is discretized into 3744 points and hexahedral background meshes are generated for the numerical integration. We employ 40 Gauss neighbor points in the numerical test. For separate variables (u, v, w) in 3D cantilever beam, Eq. 108 include the differential operators ∇u , ∇v,∇w, ∇ 2 u, ∇ 2 v, ∇ 2 w, ∇ 3 u, ∇ 3 v, ∇ 3 w . The maximal order of partial derivatives in Eq. 108 is three, hence we select the third order of nonlocal operators in Eq. 31. In addition, various gradient coefficients 1 3 = 0, 1.5, 3.0, 4.5 are studied by current method and the displacement field for various gradient coefficients is given in Fig. 13. Figure 14 presents the displacement for points on the line (y = 10, z = 0) . The results obtained by higher-order NOM correspond well with those obtained by ABAQUS and with the gradient coefficient raising, the deform of the beam become more uniformly.

Conclusions
In this paper, we present the implementation procedure of first-order and higher-order NOM, and an open-source Mathematica code is presented and explained in detail. The performance of first-order NOM and higher-order NOM results are demonstrated compare with the corresponding analytical solutions or the results of the FEM commercial software. Concluding remarks can be stated as follows: Similar to the FEM and meshless method, NOM can establish the operator energy functional and tangent stiffness matrix by some matrix multiplications. However, unlike FEM and the meshless method, NOM can derive differential operators directly without employing shape functions. Hence, the complexity of the NOM is significantly reduced. The NOM requires only the definition of the energy, for a given energy functional, the nonlocal operators can be established automatically by the highest order of partial derivative and dimensions. Support, dual-support, nonlocal differential operators, and operator energy functional are the fundamental components of NOM. According to several numerical examples, which illustrate the method's high performance and capabilities. In conclusion, NOM is an easy-to-use, flexible, and efficient numerical method, further development of this method will incorporate anisotropic material and elasticplastic material.
Appendix A: The 2D form of the nonlocal gradient operators for vector field ∇ and scalar field ∇ u Appendix B: The derivation of the hourglass tangent stiffness matrix and the global tangent stiffness matrix in 2D form (112) u i x + k 12 where j = w( ij )ΔV j (1, 1) ⊗ (1, 1) T , i = k 11 k 12 k 12 k 22 .
Finally, we can obtain the summation of the first-order global tangent stiffness matrix and hourglass tangent stiffness matrix in support S i (121)