PeriFast/Dynamics: a MATLAB code for explicit fast convolution-based peridynamic analysis of deformation and fracture

We present PeriFast/Dynamics, a compact and user-friendly MATLAB code for fast peridynamic (PD) simulations for deformation and fracture. PeriFast/Dynamics uses the fast convolution-based method (FCBM) for spatial discretization and an explicit time marching scheme to solve large-scale dynamic fracture problems. Different from existing PD solvers, PeriFast/Dynamics does not require neighbor search and storage, due to the use of the Fast-Fourier Transform and its inverse. Run-times and memory allocation are independent of the number of neighbors inside the PD horizon, leading to faster computations and lower storage requirements. The governing equations and discretization method are briefly reviewed, the code structure explained, and individual modules described in detail. A 3D demonstrative example on dynamic brittle fracture is solved using three different constitutive models (a bond-based, an ordinary state-based, and a correspondence model). Differences between the results are explained. Users are provided with a step-by-step description of the problem setup and execution of the code. PeriFast/Dynamics is a branch of the PeriFast suite of codes, and is available for download at the GitHub link provided in reference [1].


Introduction
Computational modeling of damage and fracture has been one of the most challenging areas in computational mechanics.Classical theories with the governing equations expressed in terms of partial differential equations (PDEs) are not fully capable of describing fracture since cracks are, in fact, evolving discontinuities in the continuum, and spatial derivatives at discontinuities in the displacement field are not defined.Peridynamic formulations for mechanics [2] offer alternative nonlocal approaches in which spatial derivatives are replaced with volume integrals of the primary unknowns over a certain finite region around each point, hence, allowing discontinuities (in the unknown field) to emerge and evolve in a mathematical consistent way since integration is not affected by discontinuities.PD makes seamless modeling of fracture and damage possible.In PD, cracks can naturally emerge, propagate, branch, and coalesce without the need of external, ad-hoc rules and conditions (e.g., see [3][4][5]).Significant interest on modeling fracture with PD has been observed [6][7][8].
The most common, straightforward and functional discretization for PD equations is the so-called meshfree method.In this, one approximates the integral over the nonlocal region (the PD horizon region) with a Riemann-type sum, normally using the one-point Gaussian integration, or a slight modification of that to account for nodal volumes (usually cubes) that are only partially covered by the PD horizon region [9].Note that the commercially available computer-aided engineering (CAE) software is mostly based on the finite element method (FEM) and classical PDEs.Consequently, they are inherently different from meshfree PD in terms of data structures for geometry (elements and quadrature nodes in FEM, versus nodes in meshfree PD), and in terms of solvers used, since they are based on different numerical approximation methods.There have been several attempts to manipulate commercial FEM packages to perform PD analyses (e.g., see [10,11]).Some commercial codes, e.g., LS-Dyna, have added PD capabilities as separate modules in their platform.In LS-Dyna, for example, the Discontinuous Galerkin method is used to approximate solutions to PD models ( [12,13]).U.S. National labs like Sandia and Oak Ridge National Laboratories, and research groups in academia and research labs in industry developed in-house codes for PD.Peridigm [14] is one of the few open-source PD software available from Sandia.The MOOSE-based PD code for implicit thermomechanical analysis by Idaho National Laboratory [15] is another example.
Because of its versatility in solving problems in fracture and damage, the meshfree method with direct summation for the quadrature is adopted by most existing PD in-house codes.In this approach, at every node, a loop is performed over all nodes in its "family" (neighboring nodes positioned within a finite size distance from the current node).If N is the total number of nodes and M is the number of nodes in the family of an arbitrary node, the nested loops result in solvers with the computational complexity of, at best, O(NM) .In 3D PD simulations with coarsest grids, M is at least in order of hundreds, which make PD simulations costly when compared with, for example, FEM solvers for corresponding local models.Using FEM solvers for PD is, obviously, an option but the complexity would be the same; in addition, FEM solvers are not practical for solving problems with discontinuities.That is where the advantage of the meshfree method comes in.These observations show the need for faster solvers for PD models, especially for problems involving discontinuities, like fracture and damage.
Various attempts have been made to reduce the cost of PD simulations.One popular approach is the local-nonlocal coupling where only areas around cracks are modeled by PeriFast/Dynamics analysis aims to solve PD equations of dynamic deformation and damage subjected to initial conditions (IC) and volume constraints (VC), a.k.a.nonlocal boundary conditions.Consider a 3D peridynamic body ( B ), with constrained volumes Γ 1 , Γ 2 , and Γ 3 on which the displacement field components u 1 , u 2 , and u 3 are respectively prescribed.The constrained volumes usually coincide with one another, but they do not have to.Figure 1 shows a generic 2D PD body with constrained volumes.
Let x(t) = x 1 (t), x 2 (t), x 3 (t) be the position vector of a material point at time t , with i = 1, 2, 3 corresponding to the three Cartesian coordinate directions in 3D.The PD initial- value volume-constrained (IVVC) problem for dynamics is [31] where u i is the displacement in the i-direction, v i (velocity) is the time-derivative of u i , g i is a given volume constraints on Γ i , and b i is the body/external force density in the i-direction.L i denotes the internal force density in i-direction and is defined as where H x is the finite size neighborhood of x where the nonlocal interactions pertaining to x occur.H x is known as the family or the horizon region of point x and is usually a sphere in 3D, centered at x with the radius referred to as the horizon size.x ′ denotes the position vector for family nodes in H x .f i x, x ′ , t is the dual force density: the net force between the material volume at x and the material volume at x ′ , and is determined by a PD consti- tutive model.PD material models, which define the expression for f i x, x ′ , t in Eq. ( 2), can be of two types: bond-based (BB) and state-based (SB).In BB-PD, the dual force density for each pair of nodes depends on the displacement of those nodes only, whereas in the more general SB-PD, the dual force density for each pair of nodes can depend on the deformation of the entire families of x and x ′ .In SB-PD, PD states are introduced as general nonlinear mappings, generalizations of tensors, which are linear mappings, in the (1) classical continuum mechanics theory [32].The constitutive relationships define the PD "force-state" as a function of the PD "deformation-state" and other quantities.f i in Eq. ( 2) is defined based on the force-states at x and x ′ [32].The relationship between PD force and deformation states can either be directly constructed/obtained in the nonlocal setting (the "native PD approach"), or it can be derived by a conversion (or "translation") method from a classical (local) constitutive model.The latter is known as the PD correspondence approach, which usually leads to non-ordinary state-based (NOSB) PD models.In ordinary state-based (OSB) PD models, the force vector between x and x ′ is collinear to the bond vector connecting the two points, while NOSB-PD models this does not necessarily happen [32].Correspondence models are convenient since they can use existing constitutive local models, but can suffer from numerical instabilities (zero energy modes, see [33][34][35]), and tend to have a higher computational cost than corresponding OSB ones.The constitutive model formulas used in this work are given in Appendix.
The function in Eq. ( 2) is a history-dependent bond-level damage function with the following binary definition normally used for brittle-type damage: Note that PD bonds refer to pairs of family points.A broken bond means that the interaction between the two family points that the bond connects no longer exists.In PeriFast/ Dynamics, we use the energy-based damage model proposed in [23], which is consistent with the FCBM discretization.In this model, once the strain energy density ( W ) at a point reaches a critical strain energy density ( W c ), that point loses all of its bonds irreversibly, i.e., it is completely detached from the body.The definition for in this approach can be expressed as where The definition of W(x, t) depends on the constitutive model.For the material models implemented in PeriFast, W is provided in Appendix.The threshold W c is calibrated to the critical fracture energy G 0 ∶ The details of calibration can be found in [23].Note that the calibrated formula shown in reference [23] does not include the 2 in the denominator.There, the 2 was incorporated into formula for W , leading to equivalent results as here.We prefer the current formula for clarity.
Most engineering measurements are taken on surfaces of the domain, leading to mathematical descriptions in terms of (classical) Dirichlet, Neumann, or mixed boundary conditions.In order to approximate a (classical) Dirichlet boundary condition in the PD nonlocal settings described by Eq. ( 1), one can impose displacements on a -thick volumetric layer at the boundary: this is known as the "naïve approach" [36].For more accurate enforcement of local boundary conditions in PD models, please see, e.g., [36][37][38][39][40].In the current version of PeriFast/Dynamics, we use the naïve approach.The mirror-based fictitious nodes method (FNM) [36] is also compatible with FCBM and has been implemented in the PeriFast/Corrosion branch [41].
Traction boundary conditions (Neumann type) are usually implemented as body force densities applied on a -thick layer at the corresponding boundary.Other options can be used, for example, one can specify a certain profile for g i in Eq. (1), that approximates, for example, the desired Dirichlet and Neumann boundary conditions, see [22,36].The body force approach is implemented here in PeriFast/Dynamics.
In order to be able to use the FCBM-PD, a constitutive model needs to be setup in convolutional form.For the PeriFast/Dynamics code, the linearized BB, linearized native OSB-PD model, and PD correspondence model are implemented based on the formulations presented in [23,42], where the convolutional forms for each of these constitutive models, including brittle fracture, have been derived.
While for linearized PD models and PD correspondence models of the form shown in [32], convolutional structures are easy to obtain (see [42][43][44]), a case-by-case investigation is needed for general nonlinear models to find a convolutional form to which FCBM can be applied.One example for a nonlinear bond-based model is provided in [23], while references [43,44] show the procedure for obtaining the convolution form in the case of elastoplasticity and ductile failure.

Review of the Fast Convolution-Based Discretization Method (FCBM)
PeriFast/Dynamics uses the fast convolution-based method (FCBM) to solve the PD-IVVC problem in Eq. (1).In FCBM, the convolution theorem and efficient FFT algorithms are employed to evaluate the mid-point quadrature at significantly lower costs compared to the direct summation that is traditionally used.Details of the method are given in [23] and briefly summarized below.Identification and looping over neighbors of a given node are no longer needed in FCBM, making the method independent of the neighbor numbers.The initial family search is eliminated, and memory allocation is significantly reduced, since neighbor information does not need to be stored.We aim to approximate the integral over the horizon region in Eq. (1) using mid-point integration (one-point Gaussian quadrature) but evaluated using the Fast Fourier Transform (FFT) and its inverse, instead of the regular direct summation through a nested loop over the horizon region.For FFT to be applicable in computing the convolution sums, the problem needs to be extended by periodicity to the entire space.This is done by first embedding the PD domain in a rectangular box (with a buffer of at least between the surface of the domain and the edge of the box), which is then extended by periodicity to the entire space.
Figure 2 shows the box (delineated by the dashed line) with the actual domain contained in it, extended by periodicity as depicted in Fig. 1.Note that the box edges should be at least one horizon ( ) away from the boundary of the body B .This will ensure that there will be no wrap-around effect in the circular convolution discussed below.
After extension of the body to , the following characteristic functions, are defined for distinguishing various subdomains (partitioning in a way): B is defined for eliminating any interaction between the PD body and the rest of the box and Ω i is for applying the BCs.
Using B and Ω i , the PD IVVC problem in Eq. ( 1) is modified as follows: where w i (x, t) is known from the given data: Changing the domain of integration from H x in Eq. ( 2) to in Eq. (10) does not alter the integral because f i is zero outside of the horizon region.
The solution to Eq. ( 10) on Ω i is the same as the solution to Eq. (1).Equation (10), however, is defined over a periodic domain ( ) which allows for utilizing FFT for fast evaluation of the circular convolutions arising from discretization of PD integrals.
PeriFast/Dynamics uses uniform grid spacing for spatial discretization at this stage.The discrete coordinates are defined: L 1 , L 2 and L 3 are the dimensions of the box in 3D, and N 1 , N 2 , and N 3 are number of nodes in each coordinate direction.Note that FCBM might be compatible with nonuniform discretizations if the nonuniform FFT is employed, but this is an area for future research.Using mid-point quadrature for the integral in Eq. ( 9).One gets where . Note that to compute PD integrals more accurately, one can use the partial-volume correction algorithms [45,46].These algorithms can be easily incorporated to the FCBM framework by introducing a volume correction function to Eq. ( 12).The correction functions can be defined similar to the one defined in [45].This is not done here because we will tend to use relatively large m-values (m is the ratio of horizon size to grid spacing), reducing the error in that way.An analysis of the influence of partial-volume algorithms on FBCM results is planned in the future.
The key step in FCBM is to express the summation in the equation above in terms of linear combinations of convolutions, in the following general form: where N c is a positive integer that denotes the number of convolutions, and for each l = 1, … , N c : a l is a function of point x, b l is function of x ′ , and c l is a function of (x − x ′ ) .Here c l functions are referred to as the kernel functions.Note that different constitutive models lead to different a l , b l , and c l functions that need to be defined in the code.Convo- lutional forms for the linearized bond-based, linearized native state-based, and PD correspondence models used in this work are provided in Appendix.Generally, a convolutional structure is natural for the integral operator in linear PD formulations [23,47].For nonlinear PD models, one needs to either linearize them or investigate on a case-by-case basis to see if such a structure can be found.In our previous publication [23] (see also Eq. ( 13) in the present manuscript) we showed how to obtain the convolutional structure for a large class of nonlinear PD problems.For problems that do not fall directly into this general setting, like the PD model with critical bond-strain damage criterion, we had to introduce a modified damage criterion (based on critical nodal strain energy density, instead of critical bond-strain) which allowed us to recast the formulation into that general setting and easily obtain the convolutional structure needed.Several examples for constructing convolutionbased discretizations for nonlinear PD problems have been shown for nonlinear diffusion [22], and nonlinear elasticity (bond-based) with brittle fracture [23].Notably, PD correspondence models of the form presented in [38] also fall into the general setting mentioned above and, therefore, it is easy to derive their convolutional structure (see [42][43][44]).
Using the discrete convolution theorem, Eq. ( 13) can be computed as where and −1 denote the FFT and inverse FFT operations, and c s l is the shifted kernel with respect to the box coordinates.c s l is the periodic version of c l function over , where the origin of c l is shifted to coincide with the corners of .This is necessary for the circular convolution operation to represent the PD convolution integrals.Figure 3 shows the original and the shifted version of a generic 2D radial kernel.In Section 4.6 the operation of generating the shifted kernel from the given kernel function in PeriFast is described.
By comparing Eqs. ( 13) and ( 15), we can see that the summation over the neighbors of x nmp no longer appears in the fast convolution computation, and therefore FCBM is independ- ent of the number of neighbors of a given node.As a consequence, there is no need to search, identify, and store neighbor information, leading to important CPU and storage savings.
The displacement and velocity fields at each time step Δt are updated explicitly via the velocity-Verlet algorithm (see [4] for details): Remark: in addition to the internal force density, all other PD integrals, if used (e.g., PD strain energy density), need to be expressed in the form of Eq. ( 13), in order to be computed using the FCBM.
Remark: in order to impose periodic BC in FCBM, one takes B (x) = Ω i (x) = 1 for all x .This implies that the body becomes a torus/periodic box: B ≡ Ω i ≡ .In this case, the Fourier basis functions employed in the FFT operations naturally capture the "wrap-around" effect expected in a periodic setting.This is in contrast with other discretization methods such as FEM or other meshfree methods where the periodic/wrap-around condition needs to be explicitly enforced on the boundary nodes.In the case of Fourier-based discretizations, like Fourier spectral methods and FCBM, periodic BCs are naturally captured and there is no need to explicitly enforce any type of conditions or constrains.The characteristic function (15)  introduced in our FCBM approach allows the extension of this Fourier-based method to bounded domains with non-periodic BCs.
Remark: convergence studies for FCBM have shown a quadratic rate for diffusion problems in 1D and 2D [21,22] and a super linear rate in a 3D elasticity example [23].For fracture problems, an m-convergence study in 2D was reported in [23], and a -convergence study can be found in [43,44].Convergence studies for certain variations of FCBM can also be found in [24][25][26].

PeriFast/Dynamics Code Description
In this section, we describe the data structures used in the discretization in PeriFast/Dynamics, discuss the overall structure of the code and provide details of each of its modules (m-files).

Data Structure for PD Nodes
PeriFast/Dynamics stores the PD nodal positions and nodal values for different quantities in a consistent way with MATLAB's multi-dimensional FFT operations.
Let = [x min x max ] × [y min y max ] × [z min z max ] denote a periodic box in 3D, with the uniform discretization given below: In PeriFast/Dynamics the x , y , and z-coordinates of all nodes are stored in three distinct 3D arrays of size where X jik , Y jik , and Z jik respectively denote the x, y, and z-coordinates of node x ijk .Note that the index in y-direction precedes the x-direction index, due to the way MATLAB's meshgrid function generates 3D arrays.While in traditional solvers all nodal data are usually vectorized regardless of the spatial dimension, in FCBM it is necessary to work with multi-dimensional arrays, because of using multi-dimensional FFT operations.
In PeriFast/Dynamics, functions of space and time are defined as functions of , , , and t , and return outputs in the form of 3D N 2 × N 1 × N 3 arrays containing their nodal values.
For example, let C(x) be a 3 × 3 tensor-valued function defined in 3D.The discrete version of this function in PeriFast/Dynamics is For each p and q (each component of the tensor C ), pq is a 3D N 2 × N 1 × N 3 array returned by a function of , , and .See PeriFast's nodes_and_sets.mmodule for exam- ples of such definitions.( 17) Naturally, performing the visualization during the analysis slows down the run time.For speed tests, or if solving a larger problem, it is recommended to turn off visualization_ during_analysis (in inputs.m).One can postprocess the recorded output data once the simulation is completed.The option to generate a Tecplot output (tecplot_output in inputs.m)file may also affect the speed of the analysis.
For a given problem, the user needs to specify input data in inputs.m and geometrical data in nodes_and_sets.m.Currently, the geometry data (the characteristic function, the boundary regions, the box domain coordinates) is setup manually, on a case by case basis.The users are invited to contribute functions to the code that would automate this step, for example, to directly import various CAD systems representations of the geometry data.
In this version, three material models have been implemented in PeriFast/Dynamics: (1) linearized bond-based isotropic elastic; (2) linearized state-based isotropic elastic; and (3) PD correspondence model for a hyperelastic material.We model brittle damage in all of the three cases (see Section 2 and Appendix for the damage models).New material models can be added to PeriFast/Dynamics by defining additional material types in pre_constitutive.m and constitutive.m.The user can also easily specify additional variables to output in dump_output.m(e.g., internal variables in history dependent material models) and customize visualization.m,open_Matlab_video.m,create_matlab_video.m and close_ Matlab_video.m as desired.
In the following, we take a closer look at each m-file.

Description of main.m
Box 1 shows the structure of main.m as the executable main file of the program.The file consists of three stages: reading input information, initialization, and the time computation loop.After the time loop a outputs are saved in a file named: results.mat.

Description of inputs.m
In the inputs.mfile, user-prescribed data are assigned to variables and passed onto the main program.The user needs to directly insert the input data in this file.
The terms in the parentheses denote the MATLAB variable names used in the code.Props is a 1D array, while Fb, IC_u and IC_v are structure arrays-type variables, each containing three functions corresponding to each vector component.The function of the body force's x-component for example is Fb (1).func.The variables used for traction and displacement boundary conditions are also of struct type.For BCs, however, each coordinate direction has a distinct variable associated with it, containing the number of the prescribed BCs in that direction and the corresponding functions.For example, if one needs to enforce two traction BCs in y direction, one sets trac_y.No = 2, and then define the two functions trac_y (1).func and trac_y (2).func.
The desired number of data dumps and frames for visualization is selected by the user through variables number_of_data_dump and number_of_visualization_frames.
Variable tecplot_output can be either 1 or 0. Choosing value 1 leads to the selected outputs being saved as a Tecplot file during the analysis (in the current version of the code we choose damage index as an output to be saved in a Tecplot file as an example; users can select any desired outputs).Using value 0 cancels the Tecplot output.Variable visualization_during_analysisis either 0 or 1, with 1 requesting Matlab visualization during the analysis phase, and 0 leaving out run-time visualization.Note that plots/animations can be obtained by postprocessing the output saved in Results.mat file.The variable visualization_during_analysis can be set to 0 so that the results can be plotted/animated by running postprocess.mafter saving Results.mat.This is, in fact, the recommended option when solving larger problems since plotting during analysis slows down the solver.In this version of the code the default option sets the visualization_during_analysisis variable to 0.
The user can choose the desired outputs among dumped outputs (dump_output.m)to be plotted/animated by defining a vector of integers, outputs_var_for_visualization in inputs.m.A number is assigned for each output and used for defining outputs_var_for_ visualization in inputs.m.In the current version of PeriFast/Dynamic code, u 1 ,u 2 ,u 3 ,u_ mag,v 1 ,v 2 ,v 3 , v_mag, W, d (i.e., displacement vector components, displacement magnitude, velocity vector components, velocity magnitude, strain energy density, damage index) and lambda are dumped as outputs and the assigned number for them is 1 through 11, respectively.For example, if users want to visualize u 1 , u 2 , and d among these dumped outputs, they need to set outputs_var_for_visualization = [1,2,10].Note that in order to add any other outputs for visualization which is not defined in the current version of dump.output.m,users first need to add it to dump.output.m and then modify visualization.m,open_Matlab_video.m,create_Matlab_video.m and close_Matlab_video.m to visualize that as well.

Description of nodes_and_sets.m
nodes_and_sets.m contains nodal coordinates and the geometrical information of the problem.Before describing its details, we first point out how the domain extension required by Define the following variables and return as the funcƟon output: material properƟes (props) simulaƟon Ɵme (t_max) Ɵme step (dt) GPU run switch (run_in_gpu) number of dump data (number_of_data_dump) visulalizaƟon frames frequency (number_of_visualizaƟon_frames) Tecplot switch (tecplot_output) visualizaƟon during analysis switch (visualizaƟon_during_analysis) select desired outputs for being ploƩed/animated (outputs_var_for_visulazaƟon) body force density (Fb) iniƟal condiƟons: -displacements as funcƟons of space (IC_u) -velociƟes as funcƟons of space (IC_v) tracƟon boundary condiƟons (trac_x, trac_y, trac_z): -number of tracƟons -each tracƟon as a funcƟon of space and Ɵme displacement boundary condiƟons (dispBC_x, dispBC_y, dispBC_z): -number of displacement BCs -each displacement BC as a funcƟon of space and Ɵme Box 2 Structure of inputs.mFCBM (see Fig. 2) is implemented in PeriFast/Dynamics.Given a PD body ( B ) defined by the original PD-IVVC problem, it is first assumed that the body is enclosed in a rectangular box, as tight as possible to the body.This enclosing box is shown in Fig. 4 for the 2D case with dash-line.Note that this is different from the box that is repeated by periodicity.Assuming a coordinate origin, we define the coordinates of the enclosing box vertices.To construct the periodic box , the enclosing box is extended along each direction/axis with an extension at least as large as the horizon size to avoid the "wrap-around" effect in the circular convolution.Figure 4 shows a PD body, the enclosing box (the dash-line), and the extended periodic box in 2D.l e denotes the extension length (which should be selected larger than ).Note that the best choice of the enclosing box (and the coordinate system directions) is the one that leads to the least extra space between the body and the box.Considering a fixed nodal spacing, less gap results in less excess degrees of freedom in FCBM.If the body itself is a rectangular box, then the enclosing box would be the body itself.
Box 3 shows the structure of nodes_and_sets.m.
In nodes_and_sets.m,the user first defines the horizon size ( ), the enclosing box dimensions, the extension length ( l e in Fig. 4), and the number of nodes in each direc- tion.The program then extends the enclosed box to find , and then create nodes according to Eqs. ( 14) and (15).Next, the various characteristic functions/node sets are defined by the user to describe different subdomains corresponding to the original body, traction forces, volume constraints, and pre-damage.At the end, node sets representing displacement BCs in the same directions are merged to form three distinct node sets Γ 1 , Γ 2 , Γ 3 .Then Ω 1 , Ω 2 , Ω 3 are obtained by Eq. ( 5).The horizon, box info, nodal coordinates, and the characteristic functions are passed onto main.m to be used in the analysis.chit_x,chit_y, chit_z, and chiG_x, chiG_y, chiG_z are all struct type variables and include the number of BCs in their specific direction, as well as the node sets for each of those.For example, if there are two traction BCs given in the y direction, one needs to set chit_y.No = 2, and define chit_y (1).set and chit_y (2).set,where each of these sets are 3D N 2 × N 1 × N 3 arrays with value 1 for nodes in the node set and zero elsewhere.

Description of pre_constitutive.m
This m-file contains the time-invariant functions needed for evaluation of the PD constitutive terms such as the internal force, strain energy, etc., available in the form of form of Eqs.(10) and (11).For most well-known material models, kernel functions ( c l is Eq. ( 10)) are invariant in time and should be defined in this module.Note that this module returns the FFT of the kernel functions in their shifted forms ( c s l ) described in previous section (see Fig. 3).Box 4 gives the structure of pre_constitutive.m.
Here is how the "shift operation" shown in Fig. 4 is carried out in PeriFast/Dynamics: for obtaining c s l , first, c l is translated such that its origin coincides with the center of the box: c l − x c , − y c , − z c .Then the fftshift MATLAB function is used on the trans- lated c l .The fftshift command breaks down the array from mid-planes of the box and swap the partitions, resulting in the desired shifted form: c s l .More information on fftshift is provided in the MATLAB documentation.
The coded PD correspondence model for the hyperelastic material (material ID = 2) uses St. Venant-Kirchhoff classical model for finite deformation elasticity.The implemented correspondence model includes the stability term introduced in [33] to suppresses zero energy modes and stabilizing the PD correspondence solutions.

Description of constitutive.m
This module takes the displacement field, history-dependent variables such as the old damage parameter, material properties (defined in inputs.m),discretization info (defined in nodes_and_sets.m),and the invariant parameters the constitutive response (from pre_ constitutive.m)as inputs, and returns the internal force density, strain energy density, and updated history-dependent variables (e.g.damage) as outputs.Box 5 presents the structure of this module.Note that user-defined material models are allowed in PeriFast/Dynamics and can be introduced by defining appropriate functions in pre_constitutive.m and constitutive.m,with additional material IDs, in the If-statements.
While PeriFast/Dynamics can adopt different user-defined damage models along with the user-defined constitutive laws, in the current version, for the three included constitutive models, we used the same energy-based pointwise damage model introduced in [23].In this damage model, the parameter that store damage information is a binary variable denoted by lambda which is 0 for a damaged node and 1 otherwise.A damaged node is a node for which its strain energy density exceeds a threshold calibrated to the critical fracture energy of the material.The damage index (here tracked by the variable named damage) varies between 0 and 1 and it is computed from lambda using the following relation [23]: In Eq. ( 20), is the influence function and = x ′ − x denotes the bond vector.The influence function | | is used in this work.Lambda, damage, and any other history-dependent quantities, are defined in a structuretype variable named history_var.
Remark: if one intends to study stress waves only (i.e., deformation without damage/ fracture), one can either comment out the commands corresponding to updating damage, or just prescribe a very large fracture energy value in inputs.m.

Description of update_VC.m
This module takes the displacement BCs as functions of space and time (from inputs.m),and also their corresponding node sets (from nodes_and_sets.m),and returns the nodal values for functions w i ( i = 1, 2, 3 ) in Eq. (15) as outputs.Box 6 shows the structure of update_VC.m.

Description of update_tractions.m
In PeriFast/Dynamics, traction BCs are enforced as body forces applied uniformly on a -thick layer of the body at the boundary (distributed uniformly through the thickness of the layer).The body force nodal value is obtained by dividing the traction force at a point by .
The structure of update_tractions.m is very similar to update_VC.m.

Description of initial_gpu_array.m
To accelerate computations using GPUs, one needs to convert variables involved in the convolution operations to MATLAB's "gpuarray" type using the file initial_gpu_array.m.Then, calls to MATLAB's FFT and inverse FFT functions will automatically use the GPU for these operations.Note that the Parallel Computing Toolbox needs to be installed to enable GPU computing in MATLAB.

Description of dump_output.m
This module gets the snapshot number (ks), displacements and velocities in x, y, and z directions, strain energy density, damage index, and lambda as inputs.These variables along with other post-processed quantities such as displacement magnitude, are stored in a single structure-type MATLAB variable named Output.If the visualization switch is on (if visualization_during_analysis = = 1) this variable is passed onto the visualization module for creating MATLAB plots during the analysis.The frequency of visualization of outputs is dependent on the number _of _visualization_frames defined in inputs.m.Also, if the Tecplot switch is ON in inputs.mfile, the desired output is saved as a Tecplot file (.plt).Results stored in Output can be used for any desired post-processing operation.dump_output.mmodule can be easily modified by the user to store other user-defined outputs.

Description of visualization.m and postprocess.m
This module takes the outputs from dump_output.m, the snapshot number, nodal coordinates, and the body node set, and uses them to visualize the results.This module, too, can be easily Module inputs: funcƟons of displacement BCs, their node sets, nodal coordinates, Ɵme Assemble (wx) using displacement funcƟons in x (dispBC_x) and their node sets (chiG_x) Assemble using displacement funcƟons in y (dispBC_y) and their node sets (chiG_y) Assemble using displacement funcƟons in z (dispBC_z) and their node sets (chiG_z) Return assembled volume constraints (wx, wy, wz) Box 6 Structure of update_VC.mModule inputs: funcƟons of traction BCs, their node sets, horizon, nodal coordinates, Ɵme Assemble (btx) using tracƟon funcƟons in x (trac_x), their node sets (chit_x), and horizon (delta) Assemble using tracƟon funcƟons in y (trac_y), their node sets (chit_y), and horizon (delta) Assemble using tracƟon funcƟons in z (trac_z), their node sets (chit_z), and horizon (delta) Return assembled tracƟon body force (btx, bty, btz) Box 7 Structure of update_tractions.m modified by the to plot the desired figures and/or record animations (user can select the desired output for visualization in inputs.m),and to export files in user-defined formats for further processing in external software.In order to record Matlab videos from the snapshots, create_Matlab_video.m is used.There is an option in input.m to select whether the user desires to visualize the results during the analysis or after.The default is to perform the visualization after the analysis by running postprocess.mand using the data saved in Results.matfile.

Description of open_Matlab_video.m, create_Matlab_video.m and close_
Matlab_video.m These modules are used for creating Matlab videos from the outputs.For every desired output to be animated, first, a video file needs to be opened using open_Matlab_video.m.Next, by calling create_Matlab_video.m, the sequence of frames from the desired output is written to the video file.Finally, the video file needs to be closed by using close_Matlab_video.m.In the current version of PeriFast/Dynamics, a video file for damage evolution is created.Users can easily add any other desired output for creating a video by modifying outputs_var_for_visualization in inputs.m.For example, for the nodal velocity vector components and the strain energy density, one can define outputs_var_for_visualization = [5][6][7]9], where v 1 , v 2 , v 3 , and W are assigned the indices 5, 6, 7, and 9 in this version of code.

Example of Running PeriFast/Dynamics: 3D Dynamic Analysis of Brittle Fracture in a Glass Plate
In this section, we show how a particular problem on dynamic fracture in glass is setup and run with PeriFast/Dynamics.The physical problem is an example of dynamic brittle fracture in which crack branching takes place, when the applied loading is sufficiently high.For the crack to grow straight, one needs to lower the applied stress, see below.These type of problems, until the advent of PD, have been especially difficult to correctly simulate [12,49].

Problem Setup
We consider a thin single-edge glass plate of size 0.1 × 0.04 × 0.002 m 3 with a pre-crack, subjected to sudden uniaxial tensile stress of 0 = 4 MPa on its top and bottom edges (see Fig. 5).These types of boundary conditions are not easily replicated in experiments, with crack surface ramped-up loadings being a more realizable scenario [50].However, these To execute the code, we run the main.m file.

Visualization of Results
Figure 6 shows the damage index 3D MATLAB profiles obtained by the bond-based and the state-models (native and correspondence).Evolution of velocity fields, as well as strain energy density and damage index during fracture are provided in Videos 1, 2, and 3 for these PD models, respectively.

Explanations of Differences Between Models
The results shown in Fig. 6a for the bond-based model are similar to those obtained with a 2D plane stress simulation in [49].This is a good verification of the PeriFast/Dynamics' implementation.The slight differences between damage patterns (branching near the edge) from the three FCBM-based models stem from the small actual difference between the PD constitutive models.
Although the force density in the state-based and the bond-based models in Eqs.(A-1 and A-5) are different in general, for the linearized versions in Eqs.(A-1 and A-6, if the Poisson ratio is chosen as ¼ in the state-based model, the first term in Eq. (A-6) vanishes and the bond-based formula is recovered, for points in the bulk.These models, however, even for the one-quarter Poisson ratio value, are slightly different near surfaces.The root cause for this difference is in the different PD elastic micro-moduli computed in these two models.In the bond-based formulation (see [47]) the micro-modulus is computed based on a calibration for points in the bulk, and, assuming no surface correction is used for points near boundaries, has the value equal to 12E 4 in 3D.In the state-based formulation, the bondlevel elasticity constant, 30  m , depends on the weighted volume at a node, denoted by m .The weighted volume in our model is obtained numerically by approximating the following integral over the horizon (see [47]): We can easily show the equivalency of the elastic constants in the native state-based model ( 30m ) to the bond based micromoduli at the continuum level for points in the bulk by computing m for nodes in the bulk (over a full spherical neighborhood) and using the following influence function The domain of integration in computing m , i.e., neighborhood H x , varies, however, for nodes near surfaces, including original domain boundaries and growing crack surfaces, compared to the nodes in the bulk, leading to automatically modified bond-level elastic properties near the surfaces for the native state-based models.In other words, for points near the boundary, the function m , according to Eq. ( 21) in the state-based model, changes value, while in bond-based models, unless PD surface correction algorithms (e.g., see [37]) are enforced, the bulk parameters are used everywhere.We a state-based model, to compare with the results from the bond-based shown in Fig. 6a, by setting = 0.25 and m = 4 at all points in the domain (independent on whether they are near a boundary or not).We obtained results identical to the bond-based model.While the bond-based and native state-based models differ mostly near surfaces as described above, the correspondence model is intrinsically different from the other two, making use a "translation" between PD concepts (force and displacements maps) and classical continuum mechanics quantities (stresses and strains tensors) and employing, a local constitutive model for defining the stress-strain relationship.
As the horizon goes to zero, one expects the bond-based and native state-based models approach identical solutions since their near-the-surface differences vanish.The correspondence model, in the limit of -convergence, and for well-behaved problems, also converges to the classical solution of the corresponding problem.For problems with damage/fracture, this statement needs further investigation, which is outside the scope of the current work.
The 3D PD dynamic brittle fracture analyses, using a single processor, with over 2 × 10 6 nodes and over 660 time steps took about 1.15, 1.67, and 2.87 h to complete, with the bondbased, native state-based, and the correspondence models, respectively.When employing GPU-based calculations, the computational time is around 5 min, 6 min, and 11 min, for the three different constitutive models, respectively.Computations were performed on a Dell-Precision T7910 workstation PC, Intel(R) Xeon(R) CPU E5-2643 W v4 @3.40 GHz logical processors, and 128 GB of installed memory and NVIDIA Quadro M4000 GPU with 8 GB memory.

Summary and Possible Extensions of PeriFast/Dynamics
We introduced a compact Matlab-based code, PeriFast/Dynamics, which is an implementation of the Fast Convolution-Based Method (FCBM) for dynamic deformations and fracture problems in 3D.The current version of the code uses explicit time integration and offers three different options in terms of peridynamic (PD) material models: the linearized bond-based and ordinary state-based models for isotropic elastic materials, and the PD correspondence model for isotropic hyperelastic materials.Each of these comes with a model for brittle damage based on nodal strain energy density.The code is modularized with the explicit purpose to make it user-friendly and easier to adapt, modify, and extend to other problems.As long as the PD formulation for a particular problem can be setup to exhibit a convolutional structure, one can simply update/modify the MATLAB files defining the constitutive model for that particular problem.For example, elasto-plastic and ductile failure problems can easily be implemented with the structure of our code.The code could also be extended to include a pre-processor step that reads CAD-generated sample geometries and boundary conditions and automatically determines the characteristic functions that identify the domain and boundary regions in the computational box.
Because of the FCBM used to discretize the PD formulations, PeriFast/Dynamics' simulation run-times and memory requirements are independent of the number of neighbors of a node.Previous studies showed that the FCBM leads to speedups of tens to thousands compared against the traditional meshfree method, depending on the number of neighbors used.
We have briefly reviewed the PD governing equations for dynamic brittle fracture and the FCBM discretization, followed by describing the data structures used in the code.The general structure of PeriFast/Dynamics and detailed descriptions of each of the m-files contained in the code have been given.A demonstrative example of dynamic brittle fracture in glass in 3D, solved using three different constitutive models, has been provided, with step-by-step descriptions for input data and choices of outputs.

Possible Extensions
Note that the current version uses damage models with a single parameter, which can be calibrated to the critical fracture energy (material fracture toughness).models work well in problems with pre-cracks, but when applied to problems with no pre-cracks, a higher and higher effective strength is found if one uses smaller and smaller horizon sizes (for a discussion of how to select a "proper" horizon size please see [51,52]).For quasibrittle fracture problems in bodies without pre-cracks we recommend using (and implementing), for example, the two-parameter bond-failure model (see [53]).Such an extension is immediate by defining lambda in constitutive.mas a non-binary variable with a gradual transition from 1 to 0, capturing a softening behavior at the microscale.
To implement ductile failure models, one can use, for example, the new PD correspondence model introduced and verified in [43,44].The PeriFast version presented here uses an explicit time integration scheme (velocity Verlet) and solves dynamic problems.Implicit solvers using iterative methods such as the nonlinear conjugate gradient method have been used with FCBM before (see [23]) and can be easily added to the code to perform static and quasi-static analyses.
PeriFast/Corrosion is one branch of the PeriFast suite of Matlab-based codes that implement the FCBM for PD models.The PeriFast/Corrosion branch solves corrosion damage problems (pitting corrosion, including with formation of lacy covers) and is described in [41].By coupling the/Corrosion and/Dynamics code branches of PeriFast, one can solve, for example, stress-corrosion cracking problems like those in [5].Because the code is fast and memory requirements are relatively low, one can solve such problems for samples at engineering-relevant scales.
Another possible extension of the code presented here is to model thermomechanical fracture and damage.Using the diffusion-type solver structure implemented in the/Corrosion branch of PeriFast, one can easily write a similar solver for transient thermal transport and couple it with the mechanics code/Dynamics to simulate thermomechanical fracture.
While not immediate, other interesting extensions may be possible: (1) fracture in heterogeneous materials (these could use, for example, the masking functions used in [41] to generate a polycrystalline microstructure); (2) impact and fragmentation (contact detection algorithms would be required for such models).

Linearized bond-based elastic material model
This model is basically the linearized version of the micro-elastic solid (see [47]).The internal force density for this material is where Is the bond vector, is the relative displacement and

Fig. 1
Fig. 1 Schematic of a 2D peridynamic body ( B ), consisting of the domains Ω 1 and Ω 2 , where displacement components u 1 and u 2 are unknown, respectively, and the constrained volumes ( Γ 1 and Γ 2 ) where u 1 and u 2 are independently prescribed.(Figure adopted from [23])

Fig. 2
Fig. 2 Extension of a generic peridynamic body to a periodic box in 2D.(Figure adopted from [23])

Fig. 3
Fig. 3 On the left: a 2D generic kernel function in its original form centered at zero ( c l ).On the right: the periodic shifted version ( c s l ) used in the fast convolution on = [x min x max ] × [y min y max ] .The colored disk denotes the non-zero part of the kernel function

Fig. 4 A
Fig. 4 A 2D generic PD body ( B ), the enclosing box (shown by the dash-line), and one possible extension to the periodic box Module inputs: material properƟes, horizon, nodal coordinates, periodic box dimensions Read material ID from the properƟes Calculate the coordinates of the center of the box (x_c, y_c, z_c) If material ID = 0 (Linearized bond-based elasƟc material) -Compute the PD elasƟcity constants from material properƟes -Define kernel funcƟons (see Eq. -Perform shiŌ operaƟon on to obtain -Compute FFT of funcƟons: -Return as the module output If material ID = 1 (Linearized state-based elasƟc material) -[same procedure as BB, but with different funcƟons] If material ID = 2 (PD correspondence model for hyperelasƟc material) -[same procedure as BB, but with different funcƟons] Box 4 Structure of pre_constitutive.m Module inputs: material properƟes, displacements, history-dependent variables, horizon, body node set, nodal coordinates, nodal volume Read material ID from the properƟes If material ID = 0 (Linearized bond-based elasƟc material) -Compute the frequently repeated terms, and store for the following computaƟons -Compute the internal force density: (L1, L2, L3) (from Eq. (14)) -Compute strain energy density (W) -Update the damage parameter (lambda) given the old lambda and W -Compute the damage index (damage) -Return L1, L2, L3, W, damage, lambda If material ID = 1 (Linearized state-based elasƟc material) -[same procedure as BB, but with different funcƟons] If material ID = 2 (PD correspondence model for hyperelasƟc material) -[same procedure as BB, but with different funcƟons] Box 5 Structure of constitutive.m

Fig. 5
Fig. 5 Problem description for the 3D numerical example of dynamic brittle fracture.The thickness of the sample along the z-direction is exaggerated for visibility

4.2 The Overall Code Structure The
current version of PeriFast/Dynamics consists of 14 MATLAB m-files: main.m,inputs.m,nodes_and_sets.m,pre_constitutive.m,constitutive.m,update_tractions.m,update_VC.m,initial_gpu_arrays.m,dump_output.m,visualization.m,open_Matlab_video.m,create_matlab_video.m close_Matlab_video.m and postprocess.m.main.m is the script that executes the program.inputs.mcontains certain input data including material properties, simulation time, time steps, initial and boundary conditions, visualization parameters.nodes_and_sets.mcontains the PD horizon and discrete geometrical data including nodal coordinates and discrete characteristic functions that define various subdomains: the original body, constrained volumes, pre-damaged regions, and subregions where tractions are applied as a body force.Pre_constitutive.m, and constitutive.mcontainthematerial model information (available in the form of Eq. (13)).Functions that are independent of the field variables and time, i.e., are not changing during the simulations, are defined in pre_constitutive.m.The kernel functions are usually of this type and are defined in this module.The precomputed functions in pre_constitutive.m as well as the displacement field and other inputs are passed onto the module constitutive.m,where the internal force density, strain energy density, and damage are computed.constitutive.m is the module that is called in each time step to update the material response.Files update_tractions.
m, and update_VC.mare modules called when traction and displacement boundary conditions need to be updated, respectively.initial_gpu_arrays.m,converts variables involved in the convolution operations to MATLAB's "gpuarray" type to use GPU-based computations.Dump_output.mscript is called every several time steps (frequency defined by the user in inputs.m) to record output data into a Matlab variable and into a Tecplot 360 [48] file (user can determine in inputs.m if Tecplot file is desired).If visualization is requested by the user (in inputs.m),visualization.m is called as well to plot results in Matlab at every snapshot during the analysis (number of data dump and visualization frames can be set by the user in inputs.m).