A peridynamic-based machine learning model for one-dimensional and two-dimensional structures

With the rapid growth of available data and computing resources, using data-driven models is a potential approach in many scientific disciplines and engineering. However, for complex physical phenomena that have limited data, the data-driven models are lacking robustness and fail to provide good predictions. Theory-guided data science is the recent technology that can take advantage of both physics-driven and data-driven models. This study presents a novel peridynamics-based machine learning model for one- and two-dimensional structures. The linear relationships between the displacement of a material point and displacements of its family members and applied forces are obtained for the machine learning model by using linear regression. The numerical procedure for coupling the peridynamic model and the machine learning model is also provided. The numerical procedure for coupling the peridynamic model and the machine learning model is also provided. The accuracy of the coupled model is verified by considering various examples of a one-dimensional bar and two-dimensional plate. To further demonstrate the capabilities of the coupled model, damage prediction for a plate with a preexisting crack, a two-dimensional representation of a three-point bending test and a plate subjected to dynamic load are simulated.


Multiple linear regression
In this section, the basic concept of multiple linear regression is presented. According to Montgomery et al. [57] and Alpaydin [58], in the multiple linear regression, the numeric output y is assumed to be written as a function of several input variables, x 1 , x 2 , . . . , x N and noise as y = w 0 + w 1 x 1 + w 2 x 2 + · · · + w N x N + ε (1) where x 1 , x 2 , . . . , x N represent the input variables that can be called the regressors. The parameters w 0 , w 1 , . . . , w N represent regression coefficients that need to be determined. The parameter ε represents the noise of the model. The coefficients given in Eq. (1) can be determined based on the least square criterion [57,58]. The sum of the squares of the differences between the predicted values and the correct outputs is minimized. The minimization of the squared errors can be expressed as follows; Let's assume that we have training data as x 1(i) , x 2(i) , . . . , x N (i) , y (i) , i = 1, . . . , M. Therefore, Eq. (1) can be rewritten for every single data as The sum of squared errors can be calculated as In order to minimize the squared errors given in Eq. (3), the coefficients w j , j = 1, 2, . . . , N must satisfy the following conditions ∂ S ∂w j = 0, j = 0, 1, 2, . . . , N By substituting Eq. (3) into Eq. (4), the following relations are obtained x 1(i) x N (i) . . .
Therefore, the vector of coefficients, w, can be obtained as Fig. 1 PD material points and interaction of material points at x (k) and x ( j) [5] 3 Peridynamic theory Peridynamics (PD) is a nonlocal theory representing material behavior by using integro-differential equations. In peridynamics, a material point k located at coordinates x (k) has interactions with its surround points within a distance δ as shown in Fig. 1. The circle, H x (k) surrounding material point k with the radius of δ is called the horizon of this material point. The radius δ is called horizon size and material points inside the horizon H x (k) , except material point k, are family members of the material point k. The bond-based peridynamic equation of motion introduced by Silling [1,3] is an integro-differential equation in time and space as which can be rewritten in the discretized form as where N represents the number of family members of the material point k and j represents a family member of material point k. Each material point is identified by its coordinates, x (k) , volume V (k) and mass density ρ (k) . The displacement and acceleration are denoted by u (k) andü (k) , respectively. The term V ( j) represents the volume of the material point j, and f (k)( j) represents the vector of interaction forces. The terms b (x, t) or b (k) represent the body forces applied to the material point k.

PD equation of motion for 1D structures
The PD equation of motion given in Eq. (8b) can be rewritten for a linear elastic bar as [1,4,24,59] where b x(k) represents the axial body force applied on the material point k, ξ represents the distance between material points k and j. The term C ax represents the PD bond constant. For a static state withü (k) = 0, Eq. (9a) can be rewritten as

PD equations of motion for 2D structures
According to Silling [1], Silling and Lehoucq [4], Diyaroglu [59], Nguyen and Oterkus [31], the PD equation of motion given in Eq. (8b) can be rewritten in the linearized form for 2D structures as whereb x(k) andb y(k) represent applied force per unit area [31], x represents uniform mesh size in the PD discretized model. The term s (k)( j) represents linearized bond stretch, h represents the thickness of the plate, C represents bond constant, and x (k) , y (k) and x ( j) , y ( j) represent coordinates of material points k and j, respectively. In Eqs. (11a-b), the parameter μ (k)( j) represents the interaction state between material points k and j which can be defined as [3] μ (k)( j) = 1 interaction exists 0 otherwise The state of interaction, intact or broken, can be defined based on the critical stretch criteria as [1,5] where s c represents the critical stretch. This critical stretch can be evaluated for 2D bond-based PD model as [5] where G c represent the critical energy release rate of the material, E represents the elastic modulus and ν represents Poisson's ratio. If a static condition withü (k) =v (k) = 0 is considered, Eqs. (11a-b) can be rewritten as

PD-based machine learning model
In this section, the PD-based machine learning models for 1D and 2D structures are presented. The displacement of one material point is expressed as a linear function of displacements of its family members in the horizon size of δ = 3.015 x and the external forces applied to it. The training data are generated by using modal analysis in ANSYS. The PD-based machine learning model is obtained from the data set by using linear regression.

PD-based machine learning model for one-dimensional structures
As presented in Eq. (10), within the linear elasticity, the relation between applied body force and displacements of material points is linear. Moreover, the relationship between the displacement of a material point and its family members can also be linear. Therefore, the displacement value for a material point can be expressed as a linear function of displacements of its family members and the external force applied to it as where m i , i = 1, . . . , 7 represent coefficients that need to be determined for the ML model. The term F x(k) represents the axial force applied on the material point k, x represents mesh size in the discretized model, A represents the cross-sectional area of the 1D bar, and E represents material elastic modulus. In Eq. Note that with the mesh size, elastic modulus and cross-sectional area presented in the last term of Eq. (16a), the coefficient m 7 is expected to be independent of these geometrical and material properties. As a result, the expected ML model can be applied for any cross-sectional area, mesh size and elastic modulus.

Generating the training and testing data set
In this section, the training and testing data for Eq. (16a) are obtained from the modal analysis. The data set includes seven input variables , F x(k) x/AE) and one output variable (u (k) ). By using Eq. (9c), the last input variable in Eq. (16a) can be represented as By utilizing the PD relation given in Eq. (10), the input variable given in Eq. (17) can be calculated as In order to generate the data set for 1D structures, a bar with a length of L = 6 m a cross-sectional area of A = 0.01 m 2 and the elastic modulus of E = 2 × 10 11 N/m 2 is chosen. The bar is discretized with 6 elements (7 material points) with a uniform mesh size of x = 1 m and the LINK180 element is used by allowing only axial deformations. The bar is considered for 10 basic boundary conditions as shown in Fig. 3. By considering all possible vibration modes of the bar for 10 boundary conditions as shown in Fig. 3, 50 data sets are obtained. The data sets are arranged as shown in Table 1. Specifically, for each vibration mode of the bar, the displacements of all nodes are obtained. The displacement of node k is added into the output variable column in Table 1. The displacements of the remaining nodes ( u (k+3) ) are added into the first 6 columns in Table 1. The input variable in 7 th column is calculated by using Eq. (18).

Finding the coefficients for the ML model
The data are split into two parts for training and testing purposes. The first 45 data sets are used for training and the remainders are used for testing the accuracy of the ML model. By using linear regression, the coefficients for the relation given in Eq. (16a) are obtained as In order to evaluate the accuracy of the regression model, the obtained regression model is used to predict output for the test data. By using Eq. (3), the mean-squared error between the predicted values and the original output values in the testing data is calculated as 9.09 × 10 −20 .

Extending the PD-based ML model for dynamic problems
The machine learning model provided in Eq. (16a) can be used for static problems. The model can also be further developed by adding the inertia term for dynamic problems. The relationship provided in Eq. (16a) can be rewritten as Next, by adding the inertia term, Eq. (20b) can be extended for dynamic problems as where the linear regression coefficients are presented in Eq. (19).

PD-based machine learning model for two-dimensional structures
In this section, first, the machine learning model for 2D is presented. Next, training and testing data set are obtained from modal analyses by using ANSYS. The data set is then used for linear regression analysis to obtain the regression coefficients for the machine learning model. The displacement of the material point k can be determined based on the displacements of its family members and the applied forces. In the discretized model, the material point k with a horizon size of δ = 3.015 x has 28 family members that can be numbered in the order as shown in Fig. 4.
Similar to the 1D structure, the displacement of a material point k can be assumed as a linear function of the displacements of its family members and the applied forces as where a j , b j , m j , n j ,c, d, p, q are coefficients that need to be determined for the ML model. The terms F x(k) , F y(k) represent external forces in x, y directions that are applied on the material point k.
Note that the vectors a j , b j , m j , n j store the coefficients a j , b j , m j , n j which are associated with 1st,. . .,28th family members of the material point k. Similarly, vectors u j and v j store displacement components of 1 st , . . . , 28 th family members of the material point k as shown in Fig. 4.

Generating the training and testing data set
In this section, the training and testing data for Eq. (22) are obtained from the modal analysis. The data set includes 58 input variables and two output variables (u (k) , v (k) ).
By using Eqs. (11c, d, g), the last two terms on the right-hand side of Eqs. (22c-d) can be represented as where C represent the bond constant given in Eq. (11g). By using the PD relationships given in Eq. (15), the input variables given in Eq. (23) can be calculated as In order to generate data set for 2D structures, a 6 × 6 m 2 square plate is chosen as shown in Fig. 4. The plate has a thickness of 0.1 m, the elastic modulus of E = 2 × 10 11 N/m 2 and Poisson's ratio of ν = 1/3. The plate is discretized with a mesh size of x = 1 m. The material point k located at (x = 3 m, y = 3 m), shown in blue, and its 28 family members, shown in red, are considered for the data set. Note that the dimensions of the plate, as well as the elastic modulus, are chosen arbitrarily for obtaining the data set. The ML model is expected to be applicable for any geometry and linear elastic material.
The data set is obtained from modal analyses for the plate by using the PLANE182 element in ANSYS. The plate is considered in 16 basic boundary conditions as shown in Fig. 5. In each boundary condition, the 20 possible vibration modes of the plate are considered. The displacements of the material point k and its 28 family members are obtained and added into the data set. Therefore, the data set includes 320 deformation states and it is arranged as shown in Table 2. Specifically, in each vibration mode of the plate, the displacements of all nodes are obtained. The displacements of node k (Fig. 4), (u (k) , v (k) ), are added into the two last columns in Table 2 for the output variables. The displacements of the remaining nodes (u 1 , u 2 , . . . are added for the first 56 columns in Table 2. The input variables in 57th and 58th columns are calculated by using Eq. (24).

Finding the coefficients for the ML model
In order to obtain the coefficients for the ML model, the obtained data set is split into two parts in which the 310 deformation states are used for training purpose and the remainders are used for testing the accuracy of    Similar to the 1D model, by using the testing data, the mean-squared error between the predicted values and the original values from testing data is 6.42 × 10 −23 .

Extending the PD-based ML model for dynamic problems
As given in Eq. (25), the coefficients d (k) and p (k) are very small which can be ignored as Therefore, Eq. (22) can be rewritten as Therefore, by reintroducing the inertia forces on the left-hand side of Eqs. (27e-f), the machine learning model can be extended for dynamic problems as

Numerical implementation
As given in Sect. 4, the PD-based ML models for 1D and 2D structures are obtained for material points with full interactions with their family members. Specifically, the 1D ML model is applicable for material points with 6 interactions and the 2D ML model is applicable for material points with 28 interactions. However, for material points that have some missing family members or broken interactions such as material points near boundary surfaces or near crack surfaces, the developed ML models can produce significant errors. Moreover, generating training data for all of these special cases can be very time-consuming. Therefore, a hybrid approach that couples the ML model with the PD model is used. The behaviors of the material points with full interactions are predicted by using the ML model. Meanwhile, all other material points are predicted using the PD model. Similar to the conventional PD solution, the deformation behavior of structures can be obtained by using a meshless scheme. The domain is divided into a uniform mesh. In this study, for static and quasi-static loading conditions, the adaptive dynamic relaxation (ADR) method [60,61] is used. For dynamic problems, the explicit time integration scheme [5] is used. The numerical procedure for the 1D structure is shown in Fig. 6 for static problems and Fig. 7 for dynamic problems. The numerical procedure for a 2D structure is shown in Fig. 8 for static problems and Fig. 9 for dynamic problems.
As shown in Fig. 6 for 1D static problems, at each ADR iteration, first, the boundary and loading conditions are applied. Next, the PD and ML regions for the discretized model are defined. Specifically, the PD region is defined by material points with less than 6 interactions (near boundary regions). Meanwhile, the ML region is defined by material points with 6 interactions (regions that are far from boundaries). After defining the regions, by using the displacement values of the PD family members, the PD forces for each PD material points are calculated. Next, by solving the PD equations of motion given in Eq. Similar to static problems, for 1D dynamic problems for each time step, first the boundary and loading conditions are applied as shown in Fig. 7. Next, the PD and ML regions for the discretized model are defined by using the same method used for the 1D static problems. After defining the regions, the force densities for all material points, including PD and ML regions, are calculated. If a material point belongs to the PD regions, the PD force densities for that material point are obtained by using the displacement values of its PD family members according to Eq. (9) as Otherwise, the PD force densities for that material point are approximated by using the displacement values of its ML family members according to Eq. (21b) as or  Numerical procedure for coupled approach in 2D dynamic problems Similar to 1D dynamic case, Fig. 9 presents the algorithm for 2D dynamic problems. The force densities for all material points, including PD and ML regions, are calculated according to Eq. (11) or Eq. (28) as After the force densities for all nodes are obtained, the displacements, velocities and accelerations of material points are obtained by solving the 2D equation of motion as

Numerical results
In this section, first, the PD-based machine learning models are verified by considering various examples of 1D and 2D structures. As presented in Sect. 5, a hybrid approach for coupling the machine learning models and bond-based PD models is used. The results obtained by the coupled approach are compared to PD and finite element analysis (FEA) solutions. The FEA solutions are conducted by using ANSYS commercial software with the LINK180 element for the 1D bar and PLANE182 element for the 2D plate. To further verify the capabilities of the coupled approach, damage predictions on a plate with preexisting crack subjected to tension, on a 2D representation of a three-point bending test and on plate subjected to dynamic loading are performed.

Verification for 1D model
In this section, a 1D structure subjected to different loading conditions is investigated. The displacements of four material points on the left end and four material points on the right end of the bar are obtained by solving the PD equation of motion. On the other hand, the displacements of the remaining material points are obtained by using the machine learning model.

Bar subjected to simple axial loading
A bar with a cross-sectional area A = 0.2 × 0.2 m 2 and a length of L = 1 m is investigated as shown in Fig. 10. The bar has an elastic modulus of E = 69 × 10 9 N/m 2 . The bar is fixed on the left end and it is subjected to two different loading conditions which are the tensile load of F x = 5 × 10 7 N and the compressive load ofF x = − 5 × 10 7 N. The bar is discretized with uniform 100 integration points. To implement the boundary condition, three fictitious points [62,63] are added on the left end of the bar as shown in Fig. 10b. Displacements of three fictitious points as well as the displacement of the material point located at x = 0 are set equal to zero. Therefore, as shown in Fig. 10b, the first four material points on the left end of the bar are subjected to the zero displacement condition. In Fig. 10b, the red points represent the material points in the real region, and black points represent the material points in the fictitious region. In the FEA model, the bar is discretized with 100 elements by using the LINK180 element. Figure 11 shows the displacement variations along the bar for two loading conditions. As can be seen from the figures, the results captured by the coupled machine learning and PD models match very well with the results from FEA.

Bar subjected multiple axial loading
For further verification, another bar with different geometrical and material properties is investigated as shown in Fig. 12. The bar has a cross-sectional area of A = 0.02 m 2 and length of L = 2 m and elastic modulus of E = 3.8 ×10 9 N/m 2 . The bar is subjected to forces F 3 = 5 ×10 6 N at the right end, F 1 = 0.4F 3 at x 1 = 0.6 m and F 2 = 0.5F 3 at x 2 = 1.4 m.  The bar is discretized with uniform 200 integration points. Figure 13 shows the displacement variation along the bar. As can be seen from the figure, the results captured by the coupled model match very well with the FEA results.

Bar vibration
In this example, the vibration of a bar with a length of L = 1 m, a cross-sectional area of A = 0.05 m 2 is investigated. The bar has an elastic modulus of E = 2 × 10 11 N/m 2 , and it is discretized with a mesh size of x = 0.005 m. Initially, the bar is subjected to a displacement gradient of ∂u/∂ x = 0.1. Later, the bar is left to freely vibrate; meanwhile, the left end is fixed [5]. Figure 14 shows the displacement variation of a material point located at x = L/2. As can be seen from the figure, the results captured by the coupled model match very well with the results in both PD and FEA. Therefore, the accuracy of the coupled ML and PD models is verified.

Verification for 2D model
In this section, the coupled 2D model is verified by considering various examples of 2D structures with different geometrical and material properties.

Plate subjected to tension
A square plate with dimensions of L = W = 1 m and thickness of h = 0.01 m is investigated as shown in Fig. 15a. The plate has an elastic modulus of E = 69 × 10 9 N/m 2 and Poisson's ratio of ν = 1/3, and it is subjected to tensional force per unit length of f x = 2 × 10 8 N/m.
The plate is discretized uniformly with a mesh size of x = L/100. Similar to 1D problems, the plate is discretized with two regions as shown in Fig. 15b. The PD regions, shown in red, include three layers of material points on four edges of the plate. The displacements of these points are obtained by solving the PD equations of motion. Since the effects of near-surface boundaries are significant in the PD models, the surface correction [5] is adopted for the PD regions. The remaining material points, shown in blue, belong to the ML region. The displacements of these material points are obtained by using the machine learning model. The numerical procedure for the coupled approach is described in Sect. 5. Figures 16 and 17 present the variations of displacement components of the plate. Figure 18 shows variations of displacement components along two centerlines x = L/2 and y = W/2. As can be seen from the figures, the results captured by the coupled model match very well with FEA results. Figure 19 shows the comparison of the computational time for the bond-based PD model and the coupled approach. The relationships between run-time per time step and the total number of material points in the discretized models are considered. In the PD solution, the adaptive dynamic relaxation (ADR) method is used [60,61]. In the coupled solution, the ADR method is also used for the PD region. As can be seen from the figure, the solution by using coupled ML and PD models requires less computational time for each time step than the conventional PD solution. Therefore, by using the coupled ML and PD models, the computational cost can be reduced.

Plate with square cutout subjected to tension
In this section, a square plate with dimensions of L = W = 0.05 m and thickness of h = 0.0005 m is investigated as shown in Fig. 20a. The plate has 0.01×0.01 m 2 a square cutout in the middle, and it is subjected    to loading condition of u 0 = ±0.005 m at two ends. The plate has an elastic modulus of E = 192 × 10 9 N/m 2 and Poisson's ratio of ν = 1/3. The plate is discretized uniformly with a mesh size of x = L/100. As shown in Fig. 20b, the PD regions are shown in red and the ML regions are shown in blue. Similar to the previous example, the surface correction [5] is adopted for the PD regions. Figures 21 and 22 show the displacement fields on the plate. Figure 23 shows variations of displacement components along two centerlines x = L/2 and y = W/2. As can be seen from the figures, the results captured by using coupled ML and PD models and the FEA have a good agreement which shows the accuracy of the coupled model for 2D.  After verifying the accuracy of the hybrid approach by coupling of ML and PD models, in this section, damage predictions for 2D plates are presented. In order to properly capture the behavior of the structures with progressive damages, the PD regions and ML regions are updated adaptively. At each time step, material points with 28 intact interactions are updated and the behaviors of these material points are obtained by using the 2D ML model. On the other hand, the behaviors of material points with less than 28 intact interactions, which are either near boundary surfaces or near crack surfaces, are obtained by using PD solution. Similar to the previous examples, the surface correction [5] is adopted for the PD regions near boundary surfaces.

Plate with preexisting crack subjected to tension
In this example, a plate with dimensions of L × W = 0.5 × 0.806 m 2 and thickness of h = 0.01 m is investigated as shown in Fig. 24 [64]. The plate has a preexisting crack with a length of 2a = 0.1 min the middle. The material properties are represented by the elastic modulus of E = 2.16 × 10 11 N/m 2 , Poisson's ratio of ν = 1/3, fracture toughness K c = 70 × 10 6 MPa √ m [65], and the critical energy release rate of G c = 2.1233 × 10 4 J/m 2 . The critical stretch can be calculated by using the relation given in Eq. (26) as s c = 0.0021.
The plate is subjected to incremental displacements on top and bottom by | v| = 10 −8 m per each time step as shown in Fig. 24. In the numerical model, the plate is discretized with a mesh size of x = L/150 and quasi-static solution described in Fig. 8 is used. As shown in Fig. 24b, in order to apply loading conditions, three layers of material points, shown in red, are added on top and bottom of the plate and incremental displacements are applied to these material points as where n v (top) and n v (bot) represent displacements at the current time step of material points on the top and bottom boundaries, respectively. The terms n−1 v (top) and n−1 v (bot) represent displacements at the previous time step of material points on the top and bottom boundaries, respectively. The time step size is chosen as dt = 1s. Figure 25 shows the damage evolution on the plate. As shown in Fig. 25a, the crack starts propagating when the applied displacements equal to 3.5 × 10 −4 m. As the applied displacements increase, the crack propagates horizontally as expected and it nearly reaches two sides of the plate when applied displacements are 5.3 × 10 −4 m. This observation has good agreement with the experimental results captured by .Simonsen, Törnqvist [64]. Figure 26 shows the machine learning and PD regions at different load steps. As can be seen from the figure, the PD and ML regions are adaptively updated as damage progresses. All material points with less than 28 intact interactions are determined at each load step and they are defined as PD region. The remaining material points, shown in blue, are defined as the ML region. As shown in Fig. 26, there are a small number of material points, shown in red, that need to use PD solution. As shown in Fig. 26d, when the applied displacement equals to 5.3 × 10 −4 m, the PD region includes 1602 material points; meanwhile, the total number of material points is 37750. Therefore, the computational cost for the simulation can be reduced by using this coupled approach.

Three-point bending test
In this example, a 2D representation of a three-point bending test for a concrete beam conducted by Jenq [66] is investigated as shown in Fig. 27. The dimensions of the plate are defined by S = 0.3048 mand H = 0.0702 m. The plate has an initial crack located at X = 0.25S with a crack length of a = 0.5H . The material properties are represented by the elastic modulus of E = 30 × 10 9 N/m 2 , Poisson's ratio of ν = 1/4, and the critical energy release rate G c = 20.7368 J/m 2 [66]. The critical stretch can be calculated by using the relation given in Eq. 19) as s c = 4.4376 × 10 −4 .
The plate is discretized with uniform 301 × 58 material points. The quasi-static loading is applied by increasing the displacement by w = −10 −9 at x = 0.6S, y = H for each load step. A zero vertical displacement v = 0 is applied at x = 0.1S, y = 0 and x = 1.1S, y = 0. The numerical procedure given in Fig. 8 is used for this quasi-static problem. Figure 28 shows the damage evolution on the plate. As expected, the crack propagates towards the middle position of the specimen. As shown in Fig. 28b, the angle between the crack path and vertical axis is approximate 35 0 which shows good agreement with the experimental result [66]. Figure 29 shows the adaptive PD and ML regions at different load steps. Similar to the previous example, the PD and ML regions are automatically updated based on the progressive damages on the structure. As shown in Fig. 29c, when the applied displacement

Kalthoff experiment
In this section, the experiment presented by Kalthoff and Winkler [67], Kalthoff [65,68] for a prenotched plate subjected to dynamic load is simulated by using the coupled ML and PD models. Since the problem is symmetric, only the upper half plate with dimensions of L = W = 0.1 m and thickness of h = 0.009 m is modeled as shown in Fig. 30. The plate is made of steel with the elastic modulus of E = 2 × 10 11 N/m 2 , Poisson's ratio of ν = 1/3. The fracture toughness of steel is K I c = 70 × 10 6 Nm −3/2 [65], and the critical energy release rate is G c = 2.2714×10 4 J/m 2 . The critical stretch can be calculated by using the relation given with v 0 = 16.5 m/s, t 0 = 1 µs. The plate is discretized into 150 × 150 material points, and the solution results are obtained by using a dynamic explicit time integration scheme with a time step of 0.01 µs and total simulation time of 80 µs. Similar to the previous examples, the solution results are obtained by solving Eq. (11) for the PD region and by using Eq. (28) for the ML region. Figure 31 presents the damage evolution on the plate. As can be seen from the figure, under dynamic loading conditions, the crack propagates up 67.3 0 orientation with respect to the horizontal axis. After 80µs, the crack propagates nearly to the top edge of the plate as shown in Fig. 31d. As can be seen from the figure, the crack paths captured by coupled ML and PD models match very well with the experimental results in [65,67,68]. Figure 32 shows the PD and ML regions at different time steps. As shown in Fig. 32d for the coupling model at 80µs, 12.68% of the total number of material points (2909 per 22950 material points) belong to the PD region.

Conclusion
In this study, the PD-based machine learning models for 1D and 2D structures are presented and verified by conducting various problems. The main conclusions arising from the present study are listed below: 1. The PD-based machine learning models are obtained for linear elastic material for the first time in the PD literature by using linear regression. The linear relationships between displacements of a material point with displacements of its family members and its applied forces are presented in the analytical form which can be used straightforward. 2. The machine learning models for both 1D and 2D structures can be coupled with the traditional PD solutions.
The numerical procedures for coupling the PD models and the machine learning models are also provided for both static and dynamic problems. The results captured by the coupling models show good agreements with both PD and FEA results. 3. The capability of the coupled model is also further verified by considering various fracture problems. The results captured by the coupled model have very good agreement with the experimental results. 4. The PD and ML regions can be updated adaptively at each time step. Therefore, the coupled ML and PD model can use the advantage of the PD model in terms of capturing complex fracture problems, and the advantages of the machine learning model in terms of saving computational cost. 5. The hybrid approach proposed in this study can be further applied in the PD literature for other types of structures.