SYSTEMS, METHODS, AND MEDIA FOR MORE EFFICIENT PERIDYNAMIC MODELING OF BOUNDED DOMAINS

Abstract

In accordance with some embodiments, systems, methods, and media for efficient peridynamic modeling with bounded domains are provided. In some embodiments, a method comprises: receiving a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, wherein the plurality of parameters includes a shape of a domain Ω, a horizon size δ, and one or more boundary constraints; drawing a periodic box that includes the domain Ω and a boundary layer Γ extending from the domain Ω; imposing the one or more boundary constraints on the boundary layer Γ; performing a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and causing a result of the simulation to be presented.

Description

BACKGROUND

Peridynamics (sometimes referred to herein as PD) is a nonlocal extension of classical continuum mechanics that has been used in modelling of fracture and damage. In contrast with classical (i.e., local) models described using spatial derivatives, governing equations in PD can replace spatial derivatives with integral operators. For example, in lieu of partial differential equations used in classical models, PD can use integro-differential equations. In PD, nonlocal interactions can be considered between points within a certain distance, replacing the usual spatial derivatives with integral operators. This can eliminate the usual requirements of continuity and smoothness for problem's solution function (e.g., displacement field), allowing the emergence and evolution of discontinuities, such as cracks and damage, to be modeled using PD as natural parts of the solution to the governing equations. However, numerical solutions to such nonlocal formulations are relatively expensive, due in part to the computation of volume integrals (or area integrals in 2D) performed at each node.

Nonlocal diffusion-type equations are frequently used to describe various phenomena such as swarm of organisms, flocking of birds, pedestrian traffic, delayed reaction-diffusion in biology, material dissolution and corrosion, etc. A peridynamic diffusion equation has been used a nonlocal formulation for modeling diffusion in domains with evolving discontinuities, such as heat diffusion in a cracking domain or mass-transfer in corrosion damage propagation.

Meshfree techniques based on one-point Gaussian quadrature are perhaps the most widely used numerical method for discretizing PD equations due to its benefits in modeling damage and fracture. Continuous and discontinuous Galerkin finite element methods (FEM) are other approaches used for discretizing PD equations. Continuous FEM models for PD formulations are rarely used in simulating fracture, while the discontinuous Galerkin-based discretizations of PD models have shown some difficulties in correctly predicting dynamic brittle fracture. When these discretization techniques are used, nonlocality increases the cost for a PD model, relative to the cost of numerically approximating local models. For a certain size of the nonlocal interaction region (sometimes referred to as the “horizon” size), the computational complexity of the meshfree and FEM approaches is O(N²) with N being the total number of nodes employed to discretize the domain. Thus, for largescale problems in 3D, the cost becomes prohibitive, even on massively parallel computers.

Various attempts have been made to reduce the cost of peridynamic simulations. For example, coupling the local theory with PD is one approach that uses a local model for parts of the domain, and uses the PD model only at locations near cracks/damage. This approach does not work well for problems in which damage/cracks are (or become) widely distributed throughout the domain, such as in problems like impact fragmentation, etc.

Accordingly, new systems, methods, and media for efficient peridynamic modeling with bounded domains are desirable.

SUMMARY

In accordance with some embodiments of the disclosed subject matter, systems, methods, and media for efficient peridynamic modeling with bounded domains are provided.

In accordance with some embodiments of the disclosed subject matter, a method is provided, comprising receiving a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, wherein the plurality of parameters includes a shape of a domain Ω, a horizon size δ, and one or more boundary constraints; drawing a periodic box that includes the domain Ω and a boundary layer Γ extending from the domain Ω; imposing the one or more boundary constraints on the boundary layer Γ; performing a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and causing a result of the simulation to be presented.

In some embodiments, the non-local model comprises a peridynamic model.

In some embodiments, the plurality of parameters includes one or more initial conditions.

In some embodiments, performing the simulation comprises: discretizing the non-local model using one or more Fourier transforms, wherein integral operators associated with the model are expressed in terms of convolutions; and solving the discretized non-local model on T.

In some embodiments, the non-local model comprises a model of diffusion.

In some embodiments, the model of diffusion is non-linear.

In some embodiments, the non-local model comprises a model of deformation and fracture.

In some embodiments, the model of deformation and fracture is non-linear.

In some embodiments, the non-local model comprises a model of dissolution damage.

In some embodiments, the non-local model comprises a model of corrosion damage.

In some embodiments, the domain Ω comprises a three-dimensional (3D) shape.

In some embodiments, the time step size satisfies a constraint.

In some embodiments, the time step size Δ t satisfies Δt≤1/vβ, where v is a scalar diffusivity constant associated with the non-local model, and β corresponds to an integral of a kernel function associated with the non-local model over a region H_0 of the domain Ω, H_0 having a size based on the horizon size δ and the dimensionality of the domain Ω.

In some embodiments, causing the result of the simulation to be presented comprises: causing a visual representation of at least a portion of the non-local model at a particular time step to be presented.

In some embodiments, causing the result of the simulation to be presented comprises: causing a visual state of at least a portion of the non-local model at a particular load step to be presented.

In some embodiments, receiving the plurality of parameters associated with the peridynamic model to be simulated, comprises: receiving the plurality of parameters from a remote computing device over a wide area network (WAN).

In accordance with some embodiments of the disclosed subject matter, a system for peridynamic modeling is provided, the system comprising: at least one processor that is configured to: receive a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, wherein the plurality of parameters includes a shape of a domain Ω, a horizon size δ, and one or more boundary constraints; draw a periodic box T that includes the domain Ω and a boundary layer Γ extending from the domain Ω; impose the one or more boundary constraints on the boundary layer Γ; perform a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and cause a result of the simulation to be presented.

In accordance with some embodiments of the disclosed subject matter, a non-transitory computer readable medium containing computer executable instructions that, when executed by a processor, cause the processor to perform a method for peridynamic modeling is provided, the method comprising: receiving a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, including a shape of a domain Ω, a horizon size δ, and one or more boundary constraints; drawing a periodic box T that includes the domain Ω and a boundary layer Γ extending from the domain Ω; imposing the one or more boundary constraints on the boundary layer Γ; performing a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and causing a result of the simulation to be presented.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

Various objects, features, and advantages of the disclosed subject matter can be more fully appreciated with reference to the following detailed description of the disclosed subject matter when considered in connection with the following drawings, in which like reference numerals identify like elements.

FIG. 1 shows an example of a system for efficient peridynamic modeling with bounded domains in accordance with some embodiments of the disclosed subject matter.

FIG. 2 shows an example of hardware that can be used to implement a computing device and a server, shown in FIG. 1 in accordance with some embodiments of the disclosed subject matter.

FIG. 3 shows an example of a peridynamic body and a horizon of a generic point in the body.

FIG. 4 shows an example of a peridynamic domain and an associated boundary layer that can be used in connection with efficient peridynamic modeling of transient diffusion in accordance with some embodiments of the disclosed subject matter.

FIG. 5 shows an example of the peridynamic domain of FIG. 4 surrounded by a periodic box that can be used to extend the bounded-domain peridynamic model into a periodic domain in accordance with some embodiments of the disclosed subject matter.

FIG. 6 shows an example of a process for efficient peridynamic modeling with bounded domains in accordance with some embodiments of the disclosed subject matter.

FIG. 7 shows examples comparing numerical solutions to a 1D nonlocal diffusion problem with two Dirichlet (non-periodic) boundary conditions determined using mechanisms described herein and the exact solution for the problem at various times.

FIG. 8 shows example convergence studies in terms of spatial resolution for three time step sizes based on results determined using mechanisms described herein and results determined using another technique.

FIG. 9 shows an example convergence study in terms of spatial resolution based on results determined using mechanisms described herein for a 1D diffusion example with non-periodic boundary conditions.

FIG. 10 shows examples of computational complexity of a meshfree technique and techniques described herein for solving a 1D diffusion example with non-periodic boundary conditions.

FIG. 11 shows examples of numerical solutions of a 2D nonlinear volume-constrained peridynamic transient diffusion problem determined using mechanisms described herein for efficient peridynamic modeling with bounded domains.

FIG. 12 shows examples of convergence studies in terms of spatial and temporal resolutions for the nonlinear 2D volume-constrained peridynamic transient diffusion problem determined using mechanisms described herein for efficient peridynamic modeling with bounded domains.

FIG. 13 shows example schematics of a cube with an insulated cut-out, and an extended domain with nonlocal boundary layers for imposing boundary conditions using mechanisms described herein for efficient peridynamic modeling with bounded domains.

FIG. 14 shows examples of 3D peridynamic diffusion simulation results for a cube with an insulated cutout determined using mechanisms described herein for efficient peridynamic modeling with bounded domains.

FIG. 15 shows examples comparing the 3D solution in a cross-section determined using mechanisms described herein for efficient peridynamic modeling with bounded domains, and results generated using conventional techniques.

FIG. 16 shows examples of efficiency gains relative to m-factor values for the 3D solution determined using mechanisms described herein for efficient peridynamic modeling with bounded domains over the meshfree technique, and the speedup versus problem size (N) using mechanisms described with parallel fast Fourier transform solvers.

DETAILED DESCRIPTION

In accordance with various embodiments, mechanisms (which can, for example, include systems, methods, and media) for efficient peridynamic modeling with bounded domains are provided.

Reducing the computational cost for PD models can be achieved with Fourier-based techniques. For example, the PD operator in the PD diffusion equation can be expressed in terms of convolution integrals. Fourier-based techniques can be used to transform the convolution integral operator into a multiplication in the Fourier space. Such techniques can have a complexity of O(N log₂ N), with the bulk of the computations used to compute the Fourier transform and its inverse (e.g., via fast Fourier transform (FFT)). However, Fourier-based techniques are normally only useable in periodic domains, limiting the use of such techniques to a small subset of problems to which PD can be applied.

For example, Fourier-based techniques can be used to solve the nonlocal Allen-Cahn equation, fractional-in-space reaction diffusion, and PD diffusion and wave equations with linear operators, over periodic domains. A certain class of PD models discretized with the FEM or meshfree approaches, are amenable to fast solutions using a Fourier-based approach.

As another example, Fourier-based techniques can be applied to certain non-periodic domains and arbitrary boundary conditions using a boundary adapted spectral method (BASM), which can also be referred to as a fast convolution-based method (FCBM) for 1D linear peridynamic diffusion problems, using a volume penalization (VP) technique. However, while such an approach can produce significant efficiency gains compared with meshfree discretization for PD diffusion in 1D, it was limited to 1D and linear problems.

In some embodiments, mechanisms described herein can utilize an embedded constraint(s) that can facilitate use of Fourier-based techniques for diffusion problems in multiple dimensions (e.g., in 2D and 3D problems), and in nonlinear PD diffusion problems. Note that mechanisms described herein can also be applied to problems with linear PD diffusion. Additionally, in some embodiments, mechanisms described herein can utilize an embedded constrain(s) that can facilitate enforcement of boundary conditions and use of Fourier-based techniques for at least elastic (linear or nonlinear) deformations, fracture and damage, corrosion, and dissolution.

In some embodiments, mechanisms described herein can be used to generate constitutive models and material models with convolutional structures that can facilitate efficient peridynamic modeling with bounded domains for a variety of problems. Additionally, mechanisms described herein can configure nonlinear problems, such as non-linear deformation, corrosion damage, and fracture (e.g., dynamic fracture in a brittle material), as a convolution structure amenable to modeling using techniques described herein. Note that mechanisms described herein can also be applied to classical local models. For example, mechanisms described herein can be used with the PD correspondence model that places the classical local model into a peridynamic framework.

FIG. 1 shows an example 100 of a system for efficient peridynamic modeling with bounded domains in accordance with some embodiments of the disclosed subject matter. As shown in FIG. 1, a computing device 110 can receive model data from a model data source 102. In some embodiments, computing device 110 can execute at least a portion of a peridynamic modeling system 104 to simulate a problem (e.g., a diffusion problem, a fracture problem, a deformation problem, a corrosion problem, an erosion problem, an electrostatics problem, etc.) using peridynamic modeling techniques.

In some embodiments, computing device 110 can execute at least a portion of peridynamic modeling system 104 to efficiently (e.g., more quickly, using less computational resources, etc.) simulating a problem with bounded or periodic domains and linear or non-linear behavior using mechanisms described herein.

Additionally or alternatively, in some embodiments, computing device 110 can communicate data received from model data source 102 to a server 120 over a communication network 108, which can execute at least a portion of peridynamic modeling system 104. In such embodiments, server 120 can return information to computing device 110 (and/or any other suitable computing device) indicative of an output of a peridynamic modeling process. In some embodiments, peridynamic modeling system 104 can execute one or more portions of process 600 described below in connection with FIG. 6.

In some embodiments, computing device 110 and/or server 120 can be any suitable computing device or combination of devices, such as a desktop computer, a laptop computer, a smartphone, a tablet computer, a wearable computer, a server computer, a virtual machine being executed by a physical computing device, etc.

In some embodiments, model data source 102 can be any suitable source of model data (e.g., data related to one or more materials that can be used to perform peridynamic modeling, data related to one or more bodies that be used to perform peridynamic modeling) and/or other data that can be used to simulate a problem. In a more particular example, model data source 102 can include memory storing modeling data (e.g., local memory of computing device 110, local memory of server 120, cloud storage, portable memory connected to computing device 110, portable memory connected to server 120, etc.). In another more particular example, model data source 102 can include an application configured to generate model data (e.g., computer-aided design (CAD) software) being executed by computing device 110, server 120, and/or any other suitable computing device.

In some embodiments, model data source 102 can be local to computing device 110. For example, model data source 102 can be incorporated with computing device 110 (e.g., computing device 110 can include memory that stores model data, and/or can execute a program that generates model data). As another example, model data source 102 can be connected to computing device 110 by a cable, a direct wireless link, etc. Additionally or alternatively, in some embodiments, model data source 102 can be located locally and/or remotely from computing device 110, and can communicate model data to computing device 110 (and/or server 120) via a communication network (e.g., communication network 108).

In some embodiments, communication network 108 can be any suitable communication network or combination of communication networks. For example, communication network 108 can include a Wi-Fi network (which can include one or more wireless routers, one or more switches, etc.), a peer-to-peer network (e.g., a Bluetooth network), a cellular network (e.g., a 3G network, a 4G network, a 5G network, etc., complying with any suitable standard, such as CDMA, GSM, LTE, LTE Advanced, NR, etc.), a wired network, etc. In some embodiments, communication network 108 can be a local area network (LAN), a wide area network (WAN), a public network (e.g., the Internet), a private or semi-private network (e.g., a corporate or university intranet), any other suitable type of network, or any suitable combination of networks. Communications links shown in FIG. 1 can each be any suitable communications link or combination of communications links, such as wired links, fiber optic links, Wi-Fi links, Bluetooth links, cellular links, etc.

FIG. 2 shows an example 200 of hardware that can be used to implement computing device 110 and/or server 120 in accordance with some embodiments of the disclosed subject matter. As shown in FIG. 2, in some embodiments, computing device 110 can include a processor 202, a display 204, one or more inputs 206, one or more communication systems 208, and/or memory 210. In some embodiments, processor 202 can be any suitable hardware processor or combination of processors, such as a central processing unit (CPU), a graphics processing unit (GPU), an application specific integrated circuit (ASIC), a field-programmable gate array (FPGA), etc. In some embodiments, display 204 can include any suitable display devices, such as a computer monitor, a touchscreen, a television, etc. In some embodiments, inputs 206 can include any suitable input devices and/or sensors that can be used to receive user input, such as a keyboard, a mouse, a touchscreen, a microphone, etc.

In some embodiments, communications systems 208 can include any suitable hardware, firmware, and/or software for communicating information over communication network 108 and/or any other suitable communication networks. For example, communications systems 208 can include one or more transceivers, one or more communication chips and/or chip sets, etc. In a more particular example, communications systems 208 can include hardware, firmware and/or software that can be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, etc.

In some embodiments, memory 210 can include any suitable storage device or devices that can be used to store instructions, values, etc., that can be used, for example, by processor 202 to perform a computer vision task, to present content using display 204, to communicate with server 120 via communications system(s) 208, etc. Memory 210 can include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 210 can include random access memory (RAM), read-only memory (ROM), electronically-erasable programmable read-only memory (EEPROM), one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, etc. In some embodiments, memory 210 can have encoded thereon a computer program for controlling operation of computing device 110. For example, in such embodiments, processor 202 can execute at least a portion of the computer program to transmit model data to server 120, receive model data from server 120, efficiently simulate a physical process (e.g., diffusion, fracture, etc.), generate results related to the simulation, receive results related to the simulation from server 120, and/or present results related to the simulation. As another example, processor 202 can execute at least a portion of the computer program to implement peridynamic modeling system 104. As yet another example, processor 202 can execute at least a portion of process 600 described below in connection with FIG. 6.

In some embodiments, server 120 can include a processor 212, a display 214, one or more inputs 216, one or more communications systems 218, and/or memory 220. In some embodiments, processor 212 can be any suitable hardware processor or combination of processors, such as a CPU, a GPU, an ASIC, an FPGA, etc. In some embodiments, display 214 can include any suitable display devices, such as a computer monitor, a touchscreen, a television, etc. In some embodiments, inputs 216 can include any suitable input devices and/or sensors that can be used to receive user input, such as a keyboard, a mouse, a touchscreen, a microphone, etc.

In some embodiments, communications systems 218 can include any suitable hardware, firmware, and/or software for communicating information over communication network 108 and/or any other suitable communication networks. For example, communications systems 218 can include one or more transceivers, one or more communication chips and/or chip sets, etc. In a more particular example, communications systems 218 can include hardware, firmware and/or software that can be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, etc.

In some embodiments, memory 220 can include any suitable storage device or devices that can be used to store instructions, values, etc., that can be used, for example, by processor 212 to present content using display 214, to communicate with one or more computing devices 110, etc. Memory 220 can include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 220 can include RAM, ROM, EEPROM, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, etc. In some embodiments, memory 220 can have encoded thereon a server program for controlling operation of server 120. For example, in such embodiments, processor 212 can receive model data from computing device 110, transmit model data to computing device 110, efficiently simulate a physical process (e.g., diffusion, fracture, etc.), generate results related to the simulation, transmit results related to the simulation to computing device 110, and/or cause results related to the simulation to be presented (e.g., by computing device 110). As another example, processor 212 can execute at least a portion of the computer program to implement peridynamic modeling system 104. As yet another example, processor 212 can execute at least a portion of process 600 described below in connection with FIG. 6.

FIG. 3 shows an example of a peridynamic body and a horizon of a generic point in the body. In peridynamics, the peridynamic body can be specified as a domain Ω, and a spatial point within the domain Ω can be specified by a position vector x∈³ (e.g., when domain Ω is specified in 3D). The spatial point x can interact with other points within a finite size region (sometimes referred to herein as a “neighborhood”) centered at X. Such a region can be any suitable shape, such as a disk in 2D, or a sphere in 3D. This neighborhood can referred to as a horizon region of x, which can be specified as _x, and the radius of the neighborhood can be referred to as the horizon size, or simply, the horizon, and can be specified as δ. Points located inside _x can be denoted by x′ and can be referred to as the family of x. FIG. 3 schematically shows a peridynamic body with a generic point x, its family and its horizon. Objects that carry pairwise nonlocal interactions between x and x′ can be referred to as bonds. A bond vector for each family pair can be specified as=x′−x. In general, plain (i.e., non-bold) italic letters denote scalars, boldface italic letters denote vectors, and boldface non-italic letters denote tensors. Peridynamic states, described below, are denoted by underlined letters.

FIG. 4 shows an example of a peridynamic domain and an associated boundary layer that can be used in connection with efficient peridynamic modeling of transient diffusion in accordance with some embodiments of the disclosed subject matter.

The general nonlinear anisotropic inhomogeneous case of the classical diffusion equation can be expressed as:

$\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}=\nabla ·\left[C\left(x,t\right)·\nabla u\left(x,t\right)\right]+r\left(x,t\right)& \left(1\right)\end{array}$

where u(x,t) is a function describing the diffusing quantity (e.g., concentration in mass transport, temperature in heat transfer, etc.) at point X and time t, r is a source term, ∇ is the gradient operator, the symbol (⋅) denotes inner (dot) product, and C is the second-order diffusivity tensor, which may depend on position and/or time. Nonlinearity in EQ. (1) can arise if C depends on u (e.g., concentration-dependent diffusivity, or temperature-dependent conductivity). If C is independent of u the equation is linear. In this case, for an isotropic and homogenous diffusion defined by ∇ as the scalar diffusivity constant, EQ. (1) can become:

$\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}=\nu {\nabla }^{2}u\left(x,t\right)+r\left(x,t\right),& \left(2\right)\end{array}$

where ∇² denotes the classical Laplacian operator.

A general PD diffusion equation can be expressed as a nonlocal version of EQ. (1):

$\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}={ℒ}_{\gamma }u\left(x,t\right)+r\left(x,t\right),& \left(3\right)\end{array}$

where the classical diffusion operator ∇·(C·∇) in EQ. (1) can be replaced with the PD diffusion operator _γ:

$\begin{array}{cc}{ℒ}_{\gamma }u\left(x,t\right)=\gamma \left(x,x^{\prime },t\right)\left[u\left(x^{\prime },t\right)-u\left(x,t\right)\right]{\mathrm{dx}}^{\prime }.& \left(4\right)\end{array}$

In EQ. (4), the kernel γ can be assumed to be nonnegative and symmetric in the first two variables: γ(x,x′,t)=γ(x′,x,t). Kernel γ can define the nonlocal interaction between X and X′. PD diffusion can be nonlinear if γ depends on u(x,t) and/or u(x′,t), and otherwise can be linear.

For isotropic homogeneous diffusion with constant diffusivity, the peridynamic diffusion equation can be expressed as:

$\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}={\mathrm{\nu ℒ}}_{\mu }u\left(x,t\right)+r\left(x,t\right),& \left(5\right)\end{array}$

where _μ is the PD Laplacian operator, a nonlocal version of the Laplacian in EQ. (2). The PD Laplacian can be expressed as:

$\begin{array}{cc}{ℒ}_{\mu }u\left(x,t\right)=\mu \left(x^{\prime }-x\right)\left[u\left(x^{\prime },t\right)-u\left(x,t\right)\right]{\mathrm{dx}}^{\prime },& \left(6\right)\end{array}$

where a kernel function μ is nonnegative and symmetric (e.g., μ(x)=μ(−x)).

In order to take advantage of increased efficiency facilitated by Fourier-based techniques for solving PD diffusion problems, the PD operator can be expressed in terms of convolutions. For example, the PD Laplacian given in EQ. (6) can be decomposed into a convolution integral term and a linear term:

$\begin{array}{cc}{ℒ}_{\mu }u\left(x,t\right)=\mu \left(x-x^{\prime }\right)u\left(x^{\prime },t\right){\mathrm{dx}}^{\prime }-\left(\mu \left(x^{\prime }-x\right){\mathrm{dx}}^{\prime }\right)u\left(x,t\right),& \left(7\right)\end{array}$

Additionally, since μ(x-x′) is zero outside of _x, EQ. (7) can be rewritten to:

_μ u(x,t)=[μ*u](x,t)−βu(x,t)  (8)

where (*) denotes the convolution operator,

$\begin{array}{cc}\left[\mu *u\right]\left(x,t\right)=\mu \left(x-x^{\prime }\right)u\left(x^{\prime },t\right){\mathrm{dx}}^{\prime },& \left(9\right)\end{array}$ $\mathrm{and}:$ $\begin{array}{cc}\beta =\mu \left(x\right)\mathrm{dx}.& \left(10\right)\end{array}$

Note that for examples other than this linear case, a similar decomposition can be performed for certain classes of nonlinear problems for which such convolutional expressions exist. For example, for kernels of the form:

γ(x,x′,t)=σ(x,x′,u,u′,t)ω(x′−x,t)  (11)

with Ω(X′x,t)=Ω(x−x′,t),σ(x,x′,u,u′,t)=σ(x′,x,u′,u,t) and σ being a linear combination of functions (or products of such functions) that depend on (x,u,t) or (x′,u′,t). The decomposition into convolutions for a specific example of such PD nonlinear diffusion operators is described below in connection with EQ. (12).

If the kernel γ in EQ. (4) is defined as:

$\begin{array}{cc}\gamma \left(x,x^{\prime },t\right)=\left[\frac{f\left(x,u,t\right)h\left(x^{\prime },u^{\prime },t\right)+f\left(x^{\prime },u^{\prime },t\right)h\left(x,u,t\right)}{2}\right]\omega \left(x^{\prime }-x,t\right)& \left(12\right)\end{array}$

where u′ denotes u(x′,t) for notation simplicity. Similarly, by denoting f′ and h′ for f(x′,u′,t) and h(x′,u′,t) respectively, the nonlinear diffusion PD operator becomes:

$\begin{array}{cc}{ℒ}_{\gamma }u\left(x,t\right)=\left\{\left[\frac{f\left(x,u,t\right)h\left(x^{\prime },u^{\prime },t\right)+f\left(x^{\prime },u^{\prime },t\right)h\left(x,u,t\right)}{2}\right]\omega \left(x^{\prime }-x,t\right)\right\}\left[u\left(x^{\prime },t\right)-u\left(x,t\right)\right]{\mathrm{dx}}^{\prime }=\frac{1}{2}\left({\mathrm{fh}}^{\prime }+f^{\prime }h\right)\omega \left(x-x^{\prime },t\right)\left(u^{\prime }-u\right){\mathrm{dx}}^{\prime }=\frac{1}{2}\left({\mathrm{fh}}^{\prime }{u}^{\prime }-{\mathrm{fh}}^{\prime }u+f^{\prime }{\mathrm{hu}}^{\prime }-f^{\prime }\mathrm{hu}\right)\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }=\frac{1}{2}\left[{\mathrm{fh}}^{\prime }{u}^{\prime }\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }-{\mathrm{fh}}^{\prime }u\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }+f^{\prime }{\mathrm{hu}}^{\prime }\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }-f^{\prime }\mathrm{hu}\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }\right]=\frac{1}{2}\left[f{h}^{\prime }{u}^{\prime }\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }-\mathrm{fu}{h}^{\prime }\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }+h{f}^{\prime }{u}^{\prime }\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }-\mathrm{hu}{f}^{\prime }\omega \left(x-x^{\prime },t\right){\mathrm{dx}}^{\prime }\right]=\frac{1}{2}\left\{f\left[\omega *\left(\mathrm{hu}\right)\right]-\mathrm{fu}\left(\omega *h\right)+h\left[\omega *\left(\mathrm{fu}\right)\right]-\mathrm{hu}\left(\omega *f\right)\right\}& \left(13\right)\end{array}$

Note that EQ. (8), the linear case, can be recovered from EQ. (13) when f=h=1, and Ω=μ.

In some embodiments, mechanisms described herein can be generalized and can be applied to any integral operator with a convolutional structure. For example, as described in Appendix D, which is hereby incorporated by reference herein, mechanisms described herein can be applied to PD equations of motion (e.g., to model deformation, damage, fracture, etc.). Mechanisms described herein can be used with models that are elastic or inelastic, linear or non-linear, and dynamic or static. As another example, as described in Appendix E, which is hereby incorporated by reference herein, mechanisms described herein can be applied to dissolution problems and corrosion damage problems.

In connection with EQS. (14) to (34), the convolution form given in EQ. (8) for the linear case and in EQ. (13) for the non-linear case are used to describe a Fourier-based technique for peridynamic models of diffusion in 3D periodic domains. As described in Appendix D, similar techniques can be used to apply a Fourier-based technique for peridynamic models of deformation and/or fracture (e.g., based on peridynamic equations of motion).

As described above, EQ. (5) expresses the peridynamic diffusion equation for isotropic homogeneous PD diffusion. In order to apply Fourier-based techniques, mechanisms described herein can assume periodicity of the 3D domain in the Cartesian directions. If u(x,t) is a complex-valued scalar function defined over the periodic domain =[0, L₁]×[0, L₂]×[0, L₃], 0 can be identified with L₁, L₂, and L₃ in all three directions to create the assumed domain periodicity (techniques that can be used to extend this to arbitrary domains are described below in connection with FIG. 5). Then u(x,t) can be expressed using the following infinite Fourier series in space:

$\begin{array}{cc}u\left(x,t\right)=\sum _{k_3=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_1=-\infty }^{+\infty }\hat{u}\left(k,t\right)e^{2\mathrm{\pi \zeta }\left(\frac{k_1{x}_1}{L_1}+\frac{k_2{x}_2}{L_2}+\frac{k_3{x}_3}{L_3}\right)}& \left(14\right)\end{array}$

where k={k₁,k₂,k₃} is the integer vector of Fourier modes, ζ=√{square root over (−1)}, and:

$\begin{array}{cc}\hat{u}\left(k,t\right)=\frac{1}{L_3{L}_2{L}_1}{\int }_0^{L_3}{\int }_0^{L_2}{\int }_0^{L_1}u\left(x,t\right)e^{-2\mathrm{\pi \zeta }\left(\frac{k_1{x}_1}{L_1}+\frac{k_2{x}_2}{L_2}+\frac{k_3{x}_3}{L_3}\right)}{\mathrm{dx}}_1{\mathrm{dx}}_2{\mathrm{dx}}_3& \left(15\right)\end{array}$

are the Fourier coefficients of U for different values of k. EQ. (15) can also be referred to as the Fourier transform of u while EQ. (14) can be referred to as the inverse Fourier transform of û. The kernel function μ (see EQ. (6)) can also be expressed using a similar Fourier series if it is assumed to be a periodic complex-valued scalar function defined over :

$\begin{array}{cc}\mu \left(x\right)=\sum _{k_3=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_1=-\infty }^{+\infty }\hat{\mu }\left(k\right)e^{2\mathrm{\pi \zeta }\left(\frac{k_1{x}_1}{L_1}+\frac{k_2{x}_2}{L_2}+\frac{k_3{x}_3}{L_3}\right)},\mathrm{and}& \left(16\right)\end{array}$ $\begin{array}{cc}\hat{\mu }\left(k\right)=\frac{1}{L_3{L}_2{L}_1}{\int }_0^{L_3}{\int }_0^{L_2}{\int }_0^{L_1}\mu \left(x\right)e^{-2\mathrm{\pi \zeta }\left(\frac{k_1{x}_1}{L_1}+\frac{k_2{x}_2}{L_2}+\frac{k_3{x}_3}{L_3}\right)}{\mathrm{dx}}_1{\mathrm{dx}}_2{\mathrm{dx}}_3& \left(17\right)\end{array}$

In this case, the PD Laplacian operator of u becomes:

_μ u(x,t)=[μu](x,t)−βu(x,t)  (18)

Where () denotes the 3D circular (sometimes referred to as cyclic or periodic) convolution on :

[μu](x,t)=u(x−x′)u(x′,t)dx′  (19)

In the Fourier space, the circular convolution is transformed into multiplication of the convoluted functions:

{[μu](x,t)}=L₃L₂L₁(μ)(u)=L₃L₂L₁{circumflex over (μ)}(k)û(k,t)  (20)

where denotes the Fourier transform operation.

Using the Discrete Fourier Transform (DFT) can facilitate simulation of a problem using mechanisms described herein. Accordingly, u and μ can be approximated using truncated (finite) Fourier series. For example, if N₁, N₂, and N₃ are three integer numbers up to which the Fourier series are truncated in the three directions, then the truncated Fourier series can be represented as:

$\begin{array}{cc}{u}^N\left(x,t\right)=\sum _{k_3=-N_3/2}^{+N_3/2}\sum _{k_2=-N_2/2}^{+N_2/2}\sum _{k_1=-N_1/2}^{+N_1/2}\hat{u}\left(k,t\right)e^{-2\mathrm{\pi \zeta }\left(\frac{k_1{x}_1}{L_1}+\frac{k_2{x}_2}{L_2}+\frac{k_3{x}_3}{L_3}\right)}& \left(21\right)\end{array}$ $\begin{array}{cc}{u}^N\left(x\right)=\sum _{k_3=-N_3/2}^{+N_3/2}\sum _{k_2=-N_2/2}^{+N_2/2}\sum _{k_1=-N_1/2}^{+N_1/2}\hat{\mu }\left(k\right)e^{-2\mathrm{\pi \zeta }\left(\frac{k_1{x}_1}{L_1}+\frac{k_2{x}_2}{L_2}+\frac{k_3{x}_3}{L_3}\right)}& \left(22\right)\end{array}$

A uniform discretization of the spatial domain can be selected, with the same number of nodes as the Fourier modes in the truncations of EQS. (21) and (22) in each direction:

$\begin{array}{cc}{x}_{\mathrm{ijl}}=\left(i\Delta x_1,j\Delta x_2,l\Delta x_3\right),\mathrm{where}\Delta x_1=\frac{L_1}{N_1};\Delta x_2=\frac{L_2}{N_2};\Delta x_3=\frac{L_3}{N_3}& \left(23\right)\end{array}$ $\mathrm{and}=\left\{0,\dots ,N_1-1\right\};j=\left\{0,\dots ,N_2-1\right\};l=\left\{0,\dots ,N_3-1\right\}.$

Thus, the DFT for u and μ can be represented as:

$\begin{array}{cc}\tilde{u}\left(k,t\right)=\sum _{l=0}^{N_3-1}\sum _{j=0}^{N_2-1}\sum _{i=0}^{N_1-1}{u}^N\left(x_{\mathrm{ijl}},t\right)e^{-2\mathrm{\pi \zeta }\left(\frac{k_{1}i}{N_1}+\frac{k_{2}j}{N_2}+\frac{k_{3}l}{N_3}\right)}& \left(24\right)\end{array}$ $\begin{array}{cc}\tilde{u}\left(k\right)=\sum _{l=0}^{N_3-1}\sum _{j=0}^{N_2-1}\sum _{i=0}^{N_1-1}{u}^N\left(x_{\mathrm{ijl}}\right)e^{-2\mathrm{\pi \zeta }\left(\frac{k_{1}i}{N_1}+\frac{k_{2}j}{N_2}+\frac{k_{3}l}{N_3}\right)},& \left(25\right)\end{array}$

and the inverse discrete Fourier transform (iDFT) for ũ and {tilde over (μ)} can be represented as:

$\begin{array}{cc}{u}^N\left(x_{\mathrm{ijl}},t\right)=\frac{1}{N_3{N}_2{N}_1}\sum _{k_3=-\frac{N_3}{2}}^{\frac{N_3}{2}-1}\sum _{k_2=-\frac{N_2}{2}}^{\frac{N_2}{2}-1}\sum _{k_1=-\frac{N_1}{2}}^{\frac{N_1}{2}-1}\tilde{u}\left(k,t\right)e^{2\mathrm{\pi \zeta }\left(\frac{k_{1}i}{N_1}+\frac{k_{2}j}{N_2}+\frac{k_{3}l}{N_3}\right)}& \left(26\right)\end{array}$ $\begin{array}{cc}{u}^N\left(x_{\mathrm{ijl}}\right)=\frac{1}{N_3{N}_2{N}_1}\sum _{k_3=-N_3/2}^{+N_3/2}\sum _{k_2=-N_2/2}^{+N_2/2}\sum _{k_1=-N_1/2}^{+N_1/2}\tilde{u}\left(k\right)e^{2\mathrm{\pi \zeta }\left(\frac{k_{1}i}{N_1}+\frac{k_{2}j}{N_2}+\frac{k_{3}l}{N_3}\right)}& \left(27\right)\end{array}$

The PD Laplacian of u^N on the discretized domain can be expressed via one-point Gaussian quadrature as:

$\begin{array}{cc}\begin{array}{c}{ℒ}_{\mu }{u}^N\left(x_{\mathrm{ijl}},t\right)=\sum _{m=0}^{N_3-1}\sum _{q=0}^{N_2-1}\sum _{p=0}^{N_1-1}{\mu }^N\left(x_{\mathrm{ijl}}-x_{\mathrm{pqm}}\right)[u^N\left(x_{\mathrm{pqm}},t\right)\\ -{\mu }^N\left(x_{\mathrm{ijl}},t\right)]\Delta x_1\Delta x_2\Delta x_3\\ =\sum _{m=0}^{N_3-1}\sum _{q=0}^{N_2-1}\sum _{p=0}^{N_1-1}{\mu }^N(x_{\mathrm{ijl}}\\ -x_{\mathrm{pqm}})u^N\left(x_{\mathrm{pqm}},t\right)\Delta x_1\Delta x_2\Delta x_3\\ -\{\sum _{m=0}^{N_3-1}\sum _{q=0}^{N_2-1}\sum _{p=0}^{N_1-1}{\mu }^N(x_{\mathrm{ijl}}\\ -x_{\mathrm{pqm}})\Delta x_1\Delta x_2\Delta x_3\}{\mu }^N\left(x_{\mathrm{ijl}},t\right)\\ =\left[u^N{*}_{𝕋}{u}^N\right]\left(x_{\mathrm{ijl}},t\right)-{\beta }^N{u}^N\left(x_{\mathrm{ijl}},t\right)\end{array}& \left(28\right)\end{array}$

Using the definitions for the truncated Fourier series representations and DFT operations described above, the circular convolution of μ and u can be represented as:

$\begin{array}{cc}\begin{array}{c}\left[u^N{u}^N\right]\left(x_{\mathrm{ijl}},t\right)=\mathrm{iDFT}\left\{\mathrm{DFT}\left\{\left[u^N{u}^N\right]\left(x_{\mathrm{ijl}},t\right)\right\}\right\}\\ =\mathrm{iDFT}\{\Delta x_1\Delta x_2\Delta x_3\mathrm{DFT}[\sum _{m=0}^{N_3-1}\sum _{q=0}^{N_2-1}\sum _{p=0}^{N_1-1}{\mu }^N(x_{\mathrm{ijl}}\\ -x_{\mathrm{pqm}})u^N\left(x_{\mathrm{pqm}},t\right)]\}\\ =\mathrm{iDFT}\left[\Delta x_1\Delta x_2\Delta x_3\tilde{\mu }\left(k\right)\tilde{u}\left(k,t\right)\right]\end{array}& \left(29\right)\end{array}$

Using Fast Fourier Transform (FFT) algorithms to implement the DFT and iDFT, the cost for the DFT operation and its inverse is O(N log₂ N), where N=N₁N₂N₃ is the total number of nodes. Using a simplify notation u^N (X_ijl,t)=u_ijl,t^N and μ^N(X_ijl)=μ_ijl^N, the linear nonlocal operator in discrete form can be represented as:

$\begin{array}{cc}\begin{array}{c}{ℒ}_{\mu }{u}^N{❘}_{\mathrm{ijl},t}={\mathrm{FFT}}^{-1}\left[\Delta x_1\Delta x_2\Delta x_3\mathrm{FFT}\left({\mu }_{\mathrm{ijl}}^N\right)\mathrm{FFT}\left(u_{\mathrm{ijl},t}^N\right)\right]\\ -{\beta }^N{u}_{\mathrm{ijl},t}^N,\end{array}& \left(30\right)\end{array}$

where FFT and FFT⁻¹ refer to the FFT and inverse FFT operations, respectively, and:

$\begin{array}{cc}{\beta }^N=\sum _{l=0}^{N_3-1}\sum _{j=0}^{N_2-1}\sum _{i=0}^{N_1-1}{u}^N\left(x_{\mathrm{ijl}}\right)\Delta x_1\Delta x_2\Delta x_3,& \left(31\right)\end{array}$

For singular kernels

$\underset{x\to 0}{(\lim }\mu \left(x\right)=\infty ,$

e.g. the kernel used in the example described in connection with FIGS. 13-16), μ^Nu^N and β^N become singular which make

${ℒ}_{\mu }{u}^N{❘}_{\mathrm{ijl},t}$

in EQ. (28) appear to be undefined, but in fact it is defined. The canceling pairs of terms leads to canceling singularities in μ^Nu^N and β^N, and leads to a value of zero for

${ℒ}_{\mu }{u}^N{❘}_{\mathrm{ijl},t}$

at the singularity point (ijl=pqm in EQ. (28)). In the actual implementation described below, the value of was set to zero (e.g.,

${\mu }^N\left(x_{\mathrm{ijl}}\right){❘}_{x_{\mathrm{ijl}}=0}=0).$

For integrable kernels, β can be calculated analytically at the continuum level formulation in EQ. (18). When using the discretized formulation, however, the β integral can be computed using the same discretization as the convolution term μ^Nu^N. Otherwise, discretization inconsistency between μ^Nu^N and β^N can undermine the accuracy of the computed

${ℒ}_{\mu }{u}^N.$

Using r_ijl,t in place of r(x_ijl,t), the spatially-discretized PD diffusion equation with the Fourier-based technique can be represented as:

$\begin{array}{cc}\begin{array}{c}\frac{{\mathrm{du}}_{\mathrm{ijl},t}^N}{\mathrm{dt}}={\mathrm{vFFT}}^{-1}\left[\Delta x_3\Delta x_2\Delta x_1\mathrm{FFT}\left({\mu }_{\mathrm{ijl}}^N\right)\mathrm{FFT}\left(u_{\mathrm{ijl},t}^N\right)\right]\\ -v{\beta }^N{u}_{\mathrm{ijl},t}^N+r_{\mathrm{ijl},t}\end{array}& \left(32\right)\end{array}$

Using, for example, the forward Euler method in time for approximating the solution of the ordinary differential equation (ODE) in EQ. (32), yields:

$\begin{array}{cc}\begin{array}{c}{u}_{\mathrm{ijl},t+\Delta t}^N=u_{\mathrm{ijl},t}^N\\ +\Delta t\{{\mathrm{vFFT}}^{-1}\left[\Delta x_3\Delta x_2\Delta x_1{\tilde{\mu }}_{k_1{k}_2{k}_3}\mathrm{FFT}\left(u_{\mathrm{ijl},t}^N\right)\right]\\ -v{\beta }^N{u}_{\mathrm{ijl},t}^N+r_{\mathrm{ijl},t}\}\end{array}& \left(33\right)\end{array}$

Note that {tilde over (μ)}_k1k2k3=FFT(μ_ijl^N) can be precomputed just once, before the time integration loop, since it does not depend on time in this case.

The procedure described above for the isotropic homogeneous linear PD diffusion given in EQ. (5) can be extended to nonlinear cases in which the nonlinear PD operator is decomposable into convolutional terms. For example, in the case of the nonlinear diffusion with the time dependent properties given by EQS. (4) and (13), Fourier-based techniques can leads to the following discretization:

$\begin{array}{cc}\begin{array}{c}{u}_{\mathrm{ijl},t+\Delta t}^N=u_{\mathrm{ijl},t}^N\\ +\Delta t\{\frac{1}{2}\Delta x_3\Delta x_2\Delta x_1\{f_{\mathrm{ijl},t}^N{\mathrm{FFT}}^{-1}\left[\mathrm{FFT}\left({\omega }_{\mathrm{ijl},t}^N\right)\mathrm{FFT}\left(h_{\mathrm{ijl},t}^N{u}_{\mathrm{ijl},t}^N\right)\right]\\ -f_{\mathrm{ijl},t}^N{u}_{\mathrm{ijl},t}^N{\mathrm{FFT}}^{-1}\left[\mathrm{FFT}\left({\omega }_{\mathrm{ijl},t}^N\right)\mathrm{FFT}\left(h_{\mathrm{ijl},t}^N\right)\right]\\ +h_{\mathrm{ijl},t}^N{\mathrm{FFT}}^{-1}\left[\mathrm{FFT}\left({\omega }_{\mathrm{ijl},t}^N\right)\mathrm{FFT}\left(f_{\mathrm{ijl},t}^N{u}_{\mathrm{ijl},t}^N\right)\right]\\ -h_{\mathrm{ijl},t}^N{u}_{\mathrm{ijl},t}^N{\mathrm{FFT}}^{-1}\left[\mathrm{FFT}\left({\omega }_{\mathrm{ijl},t}^N\right)\mathrm{FFT}\left(f_{\mathrm{ijl},t}^N\right)\right]\}+r_{\mathrm{ijl},t}\}\end{array}& \left(34\right)\end{array}$

Note that the only difference between the linear and nonlinear case is the discretized diffusion operator, because the nonlinear case utilizes more FFT and FFT⁻¹ calculations. In the example described above, the nonlinear case costs four times more than the linear case. However, the complexity for the nonlinear case remains O(N log₂ N).

In general, boundary conditions in PD nonlocal approaches are defined as volume constraints, defined over a δ thick volume layer on the domain exterior. For example, FIG. 3 shows a PD domain (Ω) and a boundary layer (Γ) extending a distance (δ) from a boundary of the domain Ω. A boundary between the domain Ω and the boundary layer Γ can be defined as ∂Ω. Note that although problems are generally described herein in up to three dimensions (e.g., in 1D, 2D, or 3D), these are examples, and mechanisms described herein can be used in connection with problems in any suitable number of dimensions.

According to nonlocal vector calculus, volume-constrained peridynamic problems can be defined analogous to boundary value problems with partial differential equations (PDEs) in the local theory. For example, a volume constrained peridynamic transient diffusion problem can be generally expressed as:

$\begin{array}{cc}\{\begin{array}{ccc}\frac{\partial u\left(x,t\right)}{\partial t}={ℒ}_{\gamma }u\left(x,t\right)+r\left(x,t\right)& x\in \Omega ,& t>0\\ u\left(x,0\right)=u_0\left(\mathrm{initial}\mathrm{conditions}\right)& x\in \Omega & \\ G\left(u\right)=0\left(\mathrm{volume}\mathrm{constraints}\right)& x\in \Gamma ,& t\ge 0\end{array}& \left(35\right)\end{array}$

where G is a function that prescribes the constraints of u on Γ. A solution for a problem with such volume constraints is described below in connection with FIGS. 11 and 12.

In many engineering problems, measurements are taken at the surfaces of a body, not through a layer near the surface. A natural representation of such measurements can be implemented via local boundary conditions. Techniques for imposing local boundary conditions on peridynamic bodies have been proposed. One such technique is a fictitious layer/nodes technique (e.g., as described below in connection with FIG. 22) where a certain volume constraint is considered on Γ, leading to desired boundary conditions on ∂Ω. This technique is exact in 1D and reasonably accurate in 2D and 3D (except at regions with high gradients, corners and sharp curvatures). This is used in examples described below in connection with FIGS. 7-10 and 13-16 to enforce local boundary conditions in two PD problems, by determining the corresponding volume constraint values on the boundary layer Γ. Details for imposing Dirichlet and Neumann BCs to peridynamic models are described in Appendix A, which is hereby incorporated by reference herein in its entirety.

FIG. 5 shows an example of the peridynamic domain of FIG. 4 surrounded by a periodic box that can be used to extend the bounded-domain peridynamic model into a periodic domain in accordance with some embodiments of the disclosed subject matter. FIG. 5 illustrates a technique that can be used to extend the Fourier-based techniques described above in connection with FIG. 4, to bounded domains with arbitrary boundary conditions.

As described above, while domain periodicity is a requirement for Fourier-based techniques, modifications have been introduced for local theories to solve problems with non-periodic boundary conditions. For example, as described above, a volume penalization technique can be used to solve problems with non-periodic volume constraints in 1D peridynamic diffusion using a Fourier-based technique. In some embodiments, mechanisms described herein can be used to implement simpler, more general, and more efficient techniques to impose non-periodic boundary conditions in the solution of PD models.

Using a generic irregular domain (Ω) and its boundary layer (Γ) as an example, to utilize mechanisms described herein, for solving the following volume-constrained problem:

$\begin{array}{cc}\{\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}={ℒ}_{\gamma }u\left(x,t\right)+r\left(x,t\right)& x\in \Omega ,t>0\\ u\left(x,0\right)=u_0& x\in \Omega \\ u\left(x,t\right)=ℊ\left(x,t\right)& x\in \Gamma ,t\ge 0\end{array}& \left(36\right)\end{array}$

a periodic box () can be placed around the domain and boundary layer (e.g., Ω∪Γ), as shown in FIG. 5. The region(s) that is bounded by the periodic box, but not part of either the domain or the boundary layer (e.g., \(Ω∪Γ)) can be denoted as Λ and is sometimes referred to herein as the gap region.

The volume-constrained problem in EQ. (36) can be replaced with the following periodic problem:

$\begin{array}{cc}\{\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}=\chi \left(x,t\right)\left[{ℒ}_{\gamma }u\left(x,t\right)+r\left(x,t\right)\right]+\left[1-\chi \left(x,t\right)\right]\frac{\partial w\left(x,t\right)}{\partial t}& x\in ,t>0\\ u\left(x,0\right)=u_0& x\in \end{array}& \left(37\right)\end{array}$

where χ denotes the mask function:

$\begin{array}{cc}\chi \left(x\right)=\{\begin{array}{cc}1& x\in \Omega \\ 0& x\in \backslash \Omega =\Gamma \bigcup \Lambda \end{array}& \left(38\right)\end{array}$ $\mathrm{and}:$ $\begin{array}{cc}w\left(x,t\right)=\{\begin{array}{cc}0& x\in \Omega \\ ℊ\left(x,t\right)& x\in \Gamma \\ z\left(x,t\right)& x\in \Lambda \end{array}& \left(39\right)\end{array}$

where g(x,t) is the volume constraint given by the original problem defined in EQ. (36), and z(x,t) is a function describing gap values.

Note that both classical (local) Dirichlet and Neumann boundary conditions can be handled by the volume-constrained problem description provided in EQ. (36). For example, certain profiles of g(x,t) can be chosen such that effectively a classical Dirichlet or a classical Neumann boundary condition is enforced at the domain boundary ∂Ω. The particular type of classical boundary conditions to be enforced can be specified by a proper choice of the g function in EQ. (36). Further examples of techniques that can be used to specify g for Dirichlet or Neumann boundary conditions are provided in Appendix A as EQS. (A.2) and (A.4).

Note that unlike the constraint g, gap values (z) do not directly interact with u on Ω, since Λ and Ω are more than a horizon size apart (e.g., as shown in FIG. 5), and therefore, in principle, they can be chosen arbitrarily (e.g., zero can be chosen, for example). This may lead to non-smoothness at the Γ/Λ interface. While such a lack of smoothness in the extensions may be expected to lead to associated errors (e.g., Gibbs-type phenomenon), numerical experiments, so far, indicate no oscillations are produced. The absence of the Gibbs oscillations were also investigated when non-smooth extensions are used for elastic problems. Additionally, smooth extensions can be generated using approaches available in the literature for the local model (e.g., as described in Canuto et al., “Spectral methods,” Springer, 2006, and Bueno-Orovio et al., “Spectral methods for partial differential equations in irregular domains: the spectral smoothed boundary method,” SIAM Journal on Scientific Computing, 28 (2006) 886-900. Also, note that smoothing effects are expected to be introduced in a convolution setting.

Having replaced the original volume-constrained problem with a periodic problem, the Fourier collocation techniques described above in connection with FIG. 4 can be used to spatially discretize the periodic problem. In some embodiments, the domain can be discretized according to the EQ. (23), then the problem of EQ. (37) can be defined as becomes the following initial value problem for a system of ordinary differential equations (ODEs):

$\begin{array}{cc}\{\begin{array}{cc}{\frac{{\mathrm{du}}_{\mathrm{ijl},t}^N}{\mathrm{dt}}={\chi }_{\mathrm{ijl}}({ℒ}_{\gamma }u❘}_{\mathrm{ijl},t}+r_{\mathrm{ijl},t})+\left(1-{\chi }_{\mathrm{ijl}}\right)\frac{{\mathrm{dw}}_{\mathrm{ijl},t}}{\mathrm{dt}}& x\in ,\\ u_{\mathrm{ijl},t=0}=u_{0,\mathrm{ijl}}& x\in \end{array}& \left(40\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{\chi }_{\mathrm{ijl}}=\{\begin{array}{cc}1& x_{\mathrm{ijl}}\in \Omega \\ 0& x_{\mathrm{ijl}}\in \Gamma \bigcup \Lambda \end{array}& \left(41\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{w}_{\mathrm{ijl},t}=\{\begin{array}{cc}0& x_{\mathrm{ijl}}\in \Omega \\ {ℊ}_{\mathrm{ijl},t}& x_{\mathrm{ijl}}\in \Gamma \\ z_{\mathrm{ijl},t}& x_{\mathrm{ijl}}\in \Lambda \end{array}& \left(42\right)\end{array}$

_γu|_ijl,t is the discrete peridynamic Laplacian operator calculated via Fourier transform as described above in connection with FIG. 4. Note that the particular formula for _γu|_ijl,t can depend on the kernel inside the operator (see EQS. (30) and (34) as two examples).

Using EQ. (41), the system of ODEs in EQ. (40) can be rewritten as:

$\begin{array}{cc}\frac{{\mathrm{du}}_{\mathrm{ijl},t}^N}{\mathrm{dt}}=\{\begin{array}{cc}{{ℒ}_{\gamma }u❘}_{\mathrm{ijl},t}+r_{\mathrm{ijl},t}& x_{\mathrm{ijl}}\in \Omega \\ \frac{{\mathrm{dw}}_{\mathrm{ijl},t}}{\mathrm{dt}}& x_{\mathrm{ijl}}\in \Gamma \bigcup \Lambda \end{array}& \left(43\right)\end{array}$

In some embodiments, the system of ODEs in EQ. (43) can be solved using an ODE solver (e.g., the forward Euler scheme was used). Knowing u=w on \Ω, EQ. (43) can be (approximately) solved using the following forward Euler scheme:

$\begin{array}{cc}{u}_{\mathrm{ijl},t+\Delta t}^N=\{\begin{array}{cc}{u}_{\mathrm{ijl},t}^N+\Delta t\left({{ℒ}_{\gamma }u❘}_{\mathrm{ijl},t}+r_{\mathrm{ijl},t}\right)& x_{\mathrm{ijl}}\in \Omega \\ w_{\mathrm{ijl},t+\Delta t}& x_{\mathrm{ijl}}\in \Gamma \bigcup \Lambda \end{array}& \left(44\right)\end{array}$

Additionally, EQ. (41) can be used to re-write EQ. (44) as:

$\begin{array}{cc}\begin{array}{c}{u}_{\mathrm{ijl},t+\Delta t}^N={\chi }_{\mathrm{ijl}}\left[u_{\mathrm{ijl},t}^N+\Delta t\left({{ℒ}_{\gamma }u❘}_{\mathrm{ijl},t}+r_{\mathrm{ijl},t}\right)\right]\\ +\left(1-{\chi }_{\mathrm{ijl}}\right)w_{\mathrm{ijl},t+\Delta t}\end{array}& \left(45\right)\end{array}$

In some embodiments, techniques for imposing volume constraints described above in connection with EQS. (40) to (45) can be referred to as “embedded constraint” techniques, as the constraints can be embedded into the governing integro-differential equation.

In some embodiments, when using the forward Euler scheme to solve the system of (time dependent) ODEs in EQ. (43), a time step restriction has to be imposed. For example, using a stability analysis described in Jafarzadeh, et al., “Efficient solutions for nonlocal diffusion problems via boundary-adapted spectral methods,” Journal of Peridynamics and Nonlocal Modeling, 2 (2020) 85-110, the stability condition in the absence of a penalization term for linear diffusion results in:

$\begin{array}{cc}\Delta t\le \frac{1}{v\beta }& \left(46\right)\end{array}$

This condition implies that stability does not depend on the spatial discretization resolution, but depends on the horizon size (δ) and the kernel function. This condition is somewhat similar to the Courant-Friedrichs-Lew (CFL) condition in the meshfree method. In the volume penalization technique described in Jafarzadeh et al., (which is summarized in Appendix B, which is hereby incorporated herein by reference), higher accuracy may require a potentially larger penalization parameter which may shrink the maximum allowed time step, and consequently increase the computational cost. An advantage of embedded constraint techniques described herein compared with volume penalization is that the maximum time step size does not depend on a penalization parameter. Note that this is an example specific to the forward Euler scheme and explicit time integration for a time-dependent problem. Other approaches can be used in connection with different techniques for solving the system of ODEs and/or other types of problems, which may utilize a different time step restriction or no restriction. For example, a-stable time integration can be used without imposing a restriction on time step size. As another example, a time step restriction is not needed for problems that are not time dependent (e.g., steady state diffusion, static or quasi-static deformation, etc.). In such an example, a load step (e.g., of any suitable size) may be used. In a particular example, a peridynamic solution to a quasistatic deformation using FCBM techniques is described in Appendix D.

Another advantage of embedded constraint techniques described herein over volume penalization is that embedded constraint techniques can be applied directly to other PD equations (e.g., PD equations of motion), while the volume penalization technique cannot without modification. For example, the specifics of the volume penalization techniques depend on whether PDEs are solved with first-order or second-order partial derivatives in time.

FIG. 6 shows an example 600 of a process for efficient peridynamic modeling with bounded domains in accordance with some embodiments of the disclosed subject matter. At 602, process 600 can receive parameters associated with a peridynamic model of a system to be analyzed. In some embodiments, process 600 can receive the parameters using any suitable technique or combination of techniques. For example, process 600 can receive parameters at 602 via a user interface, such as a graphical user interface (e.g., a graphical user interface presented by a display of a computing device executing at least a portion of process 600, such as display 204 of computing device 16), or other user interface. As another example, process 600 can receive parameters at 602 from another device. In a more particular example, process 600 can receive parameters from another computing device (e.g., computing device 16 can receive parameters from server 120, or server 120 can receive parameters from computing device 16). As another more particular example, process 600 can receive the parameters from a portable storage device (e.g., a portable memory device, such as a USB flash drive, an optical disk, etc.). In such examples, the parameters can be formatted using any suitable format, such as a file (e.g., a .txt. file, a .m file, an XML, file, etc.)

In some embodiments, the parameters received at 602 can include any suitable parameters that can be used to specify characteristics of a simulation (e.g., characteristics of a problem to be solved), such as physical parameters, initial conditions, boundary conditions, number of nodes (e.g., in each dimension), and/or time step size

$\left(e.g.,\Delta t\le \frac{1}{v\beta }\right).$

For example, the parameters can include a type of simulation to be performed (e.g., a simulation of diffusion, a simulation of brittle fracture, etc.). As another example, the parameters can include information identifying one or more materials to be modeled. As yet another example, the parameters can include a shape of a domain(s) to be modeled (e.g., a shape of the domain Ω). As still another example, the parameters can include initial conditions associated with a problem to be solved (e.g., u₀=u(x, t=0)).

As a further example, the parameters can include boundary conditions to be applied (e.g., u_b=u(Ω,t)|_∂ΩD, or q_b=∇u(⋅,t)·n|_∂ΩN). As another further example, the parameters can include information indicative of an external force(s) to be applied to the model (e.g., the parameters can include a bias term and/or information that can be used to define a bias term, such as a source term r(x,t), or an external body force density term b(x,t)).

As yet another further example, the parameters can include information indicative of a function that can be used to simulate interactions between different portions of a model. In a more particular example, the parameters can include a kernel(s) associated with the material(s) to be modeled, such as γ(x) or μ(x), as described above in connection with EQS. (4), (6), (11), and (12)). As another more particular example, the parameters can include a kernel(s) associated with models of deformation and/or fracture, such as C(x) or a(x), as described in Appendix D.

As still another further example, the parameters can include any other suitable physical parameters, such as a scalar diffusivity constant (e.g., ν), a horizon size (e.g., δ), the maximum time of the simulation (e.g., t_max), and/or any other suitable physical parameters. In some embodiments, process 600 can derive additional parameters (e.g., β, an integral of a kernel function) from the received parameters.

At 604, process 600 can retrieve a convolutional function(s) and/or a bias term(s) that can be used to simulate interactions between different portions of a peridynamic model. For example, based on the type of simulation to be performed, one or more material(s) included in a body to be modeled, etc., process 600 can retrieve one or more appropriate kernel functions that can be used to simulate interactions between different portions of the model. In some embodiments, kernel functions that can be used to model various different problems (e.g., involving different materials, involving different processes, etc.) can be developed and made accessible to a device executing process 600 (e.g., as part of a software package that executes process 600, from a server via an API, etc.).

As another example, based on the type of simulation to be performed, one or more material(s) included in a body to be modeled, etc., process 600 can retrieve one or more bias terms (e.g., a source term r(x,t), an external body force density term b(x,t), etc.) that can be used to simulate interactions between an environment and the model. In some embodiments, bias terms that can be used to model various different problems (e.g., involving different materials, involving different processes, etc.) can be developed and made accessible to a device executing process 600 (e.g., as part of a software package that executes process 600, from a server via an API, etc.).

At 606, process 600 can draw a periodic box around the peridynamic model and a boundary region (e.g., Γ) using any suitable technique or combination of techniques. For example, as described above in connection with FIG. 5, process 600 can draw a box around the peridynamic model.

Additionally, in some embodiments, process 600 can initialize the peridynamic model using any suitable technique or combination of techniques. For example, process 600 can calculate grid spacing (e.g., for a 3D model, process 600 can calculate grid spacing long three axes as:

$\Delta x_1=\frac{L_1}{N_1},\Delta x_2=\frac{L_2}{N_2},\Delta x_3=\frac{L_3}{N_3},$

where N_n is the number of nodes along an axis, and L_n is a length of the periodic box along the axis). Process 600 can discretize the box : x_ijl={x_i, x_j, x_i} where x_i=x₁⁰+(i−1)Δx₁ and i=1, 2, . . . , N₁, x_j=x₂⁰+(j−1)Δx₂ and j=1, 2, . . . , N₂, and x_l=x₃⁰+(l−1)Δx₃ and l=1, 2, . . . , N₃. Process 600 can Discretize the kernel function μ_ijl=μ(x_ijl). Process 600 can set μ_ijl|_xijl=0 (e.g., for singular kernels only). Process 600 can shift the discretized periodic kernel function (e.g., if x₁⁰, x₂⁰, or x₃⁰≠0): μ_ijl^s=μ(x_i−x₁⁰, x_j−x₂⁰, x_l−x₃⁰) (note that this can be performed to align the lower bounds of the periodic box with zero, for example, to align the domain to zero for DFT formulations that are configured for such domains). Process 600 can discretize the initial condition and the source term or external body force density term (e.g., μ_ijl⁰=u₀(x_ijl, 0); r_ijl⁰=r(x_ijl, 0); b_ijl⁰=b(x_ijl, t)). Process 600 can perform a Fast Fourier transform of μ_ijl^s and u_ijl⁰:_k1k2k3=FFT(μ_ijl^s) and _k1k2k3=FFT(u_ijl⁰). Process 600 can construct the discretized mask function:

${\chi }_{\mathrm{ijl}}=\{\begin{array}{cc}1& x_{\mathrm{ijl}}\in \Omega \\ 0& x_{\mathrm{ijl}}\in \backslash \Omega =\Gamma \bigcup \Lambda \end{array}.$

Process 600 can initialize a step counter to n=0. Process 600 can initialize a current time to tⁿ=0.

At 608, process 600 can generate periodic boundary conditions based on the parameters received at 602 using any suitable technique or combination of techniques. For example, process 600 can generate volume constraints based on received local boundary conditions using techniques described above in connection with EQS. (36) to (39), and in Appendix A. In a more particular example, process 600 can construct the variable containing volume constraints and gap values (w), for example, by allocating w_ijl=0, calculating and assigning volume constraints on F based on given boundary conditions and techniques described in Appendix A: w_ijl|_Γ,t=0, and calculating and assign values on the gap (if Λ≠Ø): w_ijl|_Λ,t=0.

At 610, process 600 can perform a peridynamic simulation of the model based on parameters received at 602, the kernel(s) received at 602 and/or retrieved at 604, and the periodic boundary conditions generated at 608. For example, as described above in connection with EQS. (14) to (34), process 600 can use Fourier-based techniques to perform a peridynamic simulation from an initial time to a final time (e.g., advancing a time step Δt).

In a more particular example, process 600 can update a current time based on the time step, and update a solution for the new time until the time has reached the maximum time (e.g., tⁿ⁺¹=tⁿ+Δt) Process 600 can update a solution to the peridynamic equation at the new time (e.g., in the 3D example described below in connection with FIGS. 13-16 the solution at time step n+1 can be determined using the relationship u_ijlⁿ⁺¹=χ_ijl(u_ijlⁿ+Δt{νFFT⁻¹[Δx₃Δx₂Δx₁{tilde over (μ)}_k1k2k3FFT(u_ijlⁿ)]−νβⁿu_ijlⁿ+r_ijl})+(1−χ_ijl)w_ijlⁿ). Process 600 can update the source term or force density term (e.g., in the 3D example described below in connection with FIGS. 13-16 the source term at time step n+1 can be determined using the relationship r_ijlⁿ⁺¹=r_ijl(x_ijl, tⁿ⁺¹). Process 600 can re-apply volume constraints and gap values by updating the term w_ijlⁿ⁺¹ (e.g., based on EQ. 39). Process 600 can update the step counter (e.g., n=n+1). Process 600 can repeatedly update the solution (e.g., using the preceding techniques) until a maximum time is reached).

At 612, process 600 can present results of the peridynamic analysis using any suitable technique, and/or cause results of the peridynamic analysis to be presented. For example, process 600 can present a static state of the simulation at a particular time step or time steps. As another example, process 600 can present a series of states of the simulation at a series of time steps (e.g., as an animation). In some embodiments, process 600 can cause the results to be presented by a computing device executing process 600 (e.g., if computing device 16 is executing process 600, process 600 can cause the results to be presented by a display device associated with computing device 16). Additionally or alternatively, in some embodiments, process 600 can cause the results to be presented by another computing device (e.g., that is not executing process 600, or is executing only a portion of process 600). For example, if server 120 is executing process 600, process 600 can cause the results to be presented by computing device 16.

FIG. 7 shows examples comparing numerical solutions to a 1D nonlocal diffusion problem with two Dirichlet (non-periodic) boundary conditions determined using mechanisms described herein and the exact solution for the problem at various times.

Described below in connection with FIGS. 7-16, the accuracy and convergence of mechanisms described herein for efficient peridynamic modeling with bounded domains is numerically evaluated for various examples, including a 1D peridynamic linear diffusion example with non-periodic boundary conditions, a 2D volume-constrained peridynamic nonlinear diffusion problem, and a 3D diffusion problem with non-periodic local boundary conditions using mechanisms described herein.

To evaluate the performance of mechanisms described herein in terms of accuracy and convergence, mechanisms described herein were used to solve a PD problem for which the analytical solution is known (based on the method of manufactured solutions). Consider the following 1D PD linear diffusion problem on the bounded domain

$\Omega =\left[-\frac{L}{2},\frac{L}{2}\right]:$ $\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}=v{\int }_{{\mathscr{H}}_x}\mu \left(x^{\prime }-x\right)\left[u\left(x^{\prime }\right)-u\left(x\right)\right]{\mathrm{dx}}^{\prime }+r\left(x,t\right)& \left(47\right)\end{array}$ $\mathrm{where}:$ $\begin{array}{cc}\mu \left(x\right)=\{\begin{array}{cc}\frac{12}{{\delta }^3}\left(1-\frac{❘x❘}{\delta }\right);& ❘x❘\le \delta \\ 0;& ❘x❘>\delta \end{array}& \left(48\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\begin{array}{c}r\left(x,t\right)=v\{\frac{6L^2}{{\delta }^4{\pi }^2}\left[\cos \left(\frac{2\mathrm{\pi \delta }}{L}\right)-1\right]+\frac{12}{{\delta }^2}\\ -1\}{e}^{-\mathrm{vt}}\sin \left(\frac{2\pi x}{L}\right),\end{array}& \left(49\right)\end{array}$

with the initial condition:

$\begin{array}{cc}u\left(x,0\right)=\frac{2x}{L}+\sin \left(\frac{2\pi x}{L}\right),& \left(50\right)\end{array}$

and subjected to Dirichlet boundary conditions:

$\begin{array}{cc}u\left(-\frac{L}{2},t\right)=-1& \left(51\right)\end{array}$ $\begin{array}{cc}u\left(\frac{L}{2},t\right)=1& \left(52\right)\end{array}$

The exact analytical solution to the above problem is:

$\begin{array}{cc}u\left(x,t\right)=\frac{2x}{L}+e^{-\mathrm{vt}}\sin \left(\frac{2\pi x}{L}\right)& \left(53\right)\end{array}$

In order to use mechanisms described herein the bounded domain

$\Omega =\left[-\frac{L}{2},\frac{L}{2}\right]$

is extended by a length δ on both ends to create the boundary volumes (Γ). In the extended domain

$=[-\frac{L}{2}-\delta ,\frac{L}{2}+\delta ),$

volume constraints should be specified on

${\Gamma }_1=[-\frac{L}{2}-\delta ,-\frac{L}{2})\mathrm{and}{\Gamma }_2=(\frac{L}{2},\frac{L}{2}+\delta ]$

as the left and the right boundary layers. In this example, =Ω∪Γ is considered to be periodic; therefore no gap region is included (i.e., Λ=Ø). There are multiple alternatives for applying the local boundary conditions given in EQS. (51) and (52). For example, the exact volume constraints can be used (because the exact, manufactured solution for the nonlocal problem is known). As another example, a technique for constructing volume constraints based on the given local boundary conditions (e.g., as described in Appendix A). In order to compare with the results generated using the volume penalization technique described in Jafarzadeh et al., the second approach was used. This example also demonstrates how local boundary conditions can be applied to a general peridynamic problem, even when the nonlocal analytical solution is not known.

In this example, L=2, ν=0.2, δ=0.2, the total diffusion time t_max=15, and the expression =[−1.2, 1.2) can describe the periodic extended domain with −1.2 identified with 1.2. Additionally N=2⁹ total nodes were used with Δt=1.67×10⁻² as the maximum time step allowed by EQ. (46). FIG. 7 shows the numerical results generated using techniques described herein (labeled as FCBM-EC) versus the exact solution. As shown in FIG. 7, boundary conditions were correctly enforced using mechanisms described herein.

FIG. 8 shows example convergence studies in terms of spatial resolution for three time step sizes based on results determined using mechanisms described herein and results determined using another technique.

A convergence test was performed on this example to study the influence of spatial and temporal resolutions on the maximum relative error:

$\underset{0<t<t_{\max }}{\max }\frac{{u_{i,t}-u_{i,t}^N}_{\infty }}{{u_{i,t=0}}_{\infty }}.$

In peridynamics the ratio

$m=\frac{\delta }{\Delta x}$

is usually referred to as “m factor” which is a measure for family node density. For convergence tests in terms of spatial discretization, the spatial discretization is refined while keeping the horizon size fixed. This is sometimes referred to as “m-convergence”. The number of nodes was varied from N=2⁶ to 2¹³ (m varied from 5.3 to 682.7 accordingly). Results are presented in FIG. 8, panel (a). A similar convergence study using the volume penalization technique described in Jafarzadeh et al., and the results are shown in FIG. 8, panel (b).

Comparing results in FIG. 8, panels (a) and (b), the simulation performed using mechanisms described herein are more accurate than the volume penalization approach for the same time step size. Note that the temporal resolution in the volume penalization technique is tied to penalization parameter (E) magnitude. For each time step, the smallest possible E was used to test the optimum performance of the two approaches. For both approaches, the rate of convergence is O(Δx²) for sufficiently small Δt.

Note that in most cases, Fourier-based techniques described herein show spectral accuracy for smooth functions. From the convergence study shown in FIG. 8, mechanisms described herein produced quadratic accuracy, not “spectral” accuracy. However, the mechanisms described herein can still be referred to as “spectral” given the use of global and smooth functions in the Fourier series representations of the kernel and the unknown functions. Techniques using such basis functions are generally termed as “spectral”.

FIG. 9 shows an example convergence study in terms of spatial resolution based on results determined using mechanisms described herein for a 1D diffusion example with non-periodic boundary conditions. Although FIG. 8, panel (a) shows the error for three different time steps, a separate temporal convergence test was conducted with spatial resolution fixed at N=2¹³ (m=682.7). A linear rate of convergence, O(Δt), can be observed, which is expected from the Forward Euler time integrator. Methods with higher order accuracy (e.g., explicit 4^th order Runge-Kutta) can be used for temporal integration. In that case, the stability restriction on time step size for explicit methods would be different from that shown in EQ. (46).

FIG. 10 shows examples of computational complexity of a meshfree technique and techniques described herein for solving a 1D diffusion example with non-periodic boundary conditions. To check the computational complexity for the meshfree method using direct quadrature, and mechanisms described herein using the FFT, the 1D problem described above in connection with FIG. 7 was solved, considering various spatial discretization sizes. The meshfree method was performed using N=2⁹ to 2²⁰ and a simulation utilizing mechanisms described herein was performed with N=2¹⁶ to 2²⁸. The time step in all these tests was Δt=1.67×10⁻². Simulations were conducted using a Dell-Precision T7810 workstation PC, Intel® Xeon® CPU E5-2687 W v3@3.10 GHz logical processors, and 64 GB of installed memory. We utilized only one CPU for these computations. FIG. 7 shows the computational run times for each method, relative to the problem size.

As expected, the complexity of meshfree method with direct quadrature was O(N²) (e.g., as shown in FIG. 10, panel (a)), while for the simulation performed using mechanisms described herein FCBM-EC we find it to be O(N log₂ N) (e.g., as shown in FIG. 10, panel (b)).

FIG. 11 shows examples of numerical solutions of a 2D nonlinear volume-constrained peridynamic transient diffusion problem determined using mechanisms described herein for efficient peridynamic modeling with bounded domains. The images in the left column of FIG. 11 show the solution to the 2D problem described below generated using mechanisms described herein at three times (t=1, 5, and 15), and the right column shows the corresponding relative error distribution compared with the exact nonlocal solution.

The method of manufactured solutions was used to construct a 2D nonlinear volume-constrained peridynamic transient diffusion problem. Consider the following function:

u(x,t)=e^−0.2t(1−x₁²)(1−x₂²),x={x₁,x₂}∈² and t≥0  (54)

EQ. (54) is the analytical solution to the following nonlinear problem:

$\begin{array}{cc}\{\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}={ℒ}_{\gamma }u\left(x,t\right)+r\left(x,t\right)& x\in \Omega ={\left[-1,1\right]}^2,t>0\\ u\left(x,0\right)=\left(1-x_1^2\right)\left(1-x_2^2\right)& x\in \Omega ={\left[-1,1\right]}^2\\ u\left(x,t\right)=e^{-0.2t}\left(1-x_1^2\right)\left(1-x_2^2\right)& x\in \Gamma ={\left[-1.2,1.2\right]}^2{\text{[}-1,1]}^2,t\ge 0\end{array}& \left(55\right)\end{array}$ $\mathrm{Where}$ $\begin{array}{cc}{ℒ}_{\gamma }u\left(x,t\right)=\gamma \left(x^{\prime },x\right)\left[u\left(x^{\prime },t\right)-u\left(x,t\right)\right]{\mathrm{dx}}^{\prime }& \left(56\right)\end{array}$ $\mathrm{with}$ $\begin{array}{cc}\gamma \left(x^{\prime },x\right)=u\left(x,t\right)u\left(x^{\prime },t\right)\omega \left(x^{\prime }-x\right)& \left(57\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\omega \left(x\right)=0.2\left(1-\frac{x}{\delta }\right).& \left(58\right)\end{array}$

Note that this is a special case of EQ. (13) where f=h=u. If U is temperature for example, with this kernel, EQ. (55) is a nonlinear PD problem with temperature dependent diffusivity. In EQ. (55):

$\begin{array}{cc}\begin{array}{c}r\left(x,t\right)=\frac{-0.2\pi e^{-0.6t}}{443520}\left(1-x_1^2\right)\left(1-x_2^2\right)[189{\delta }^{10}\\ +{\delta }^8\left(4620x_1^2+4620x_2^2-3080\right)\\ +{\delta }^6(7920x_1^4+92400x_1^2{x}_2^2-44880x_1^2\\ +7920x_2^4-44880x_2^2+23760)\\ +{\delta }^4(110880x_1^4{x}_2^2-22176x_1^4+110880x_1^2{x}_2^4\\ -443520x_1^2{x}_2^2+155232x_1^2-22176x_2^4\\ +155232x_2^2-44352)+44352e^{0.4t}]\end{array}& \left(59\right)\end{array}$

Assuming δ=0.4, consider Ω∪Γ=[−1.4,1.4]² can be the 2D periodic box to solve the problem using mechanisms described herein. Note that for this example there is also no gap region (i.e., Λ=Ø). The domain is discretized with uniform grid spacing in both directions: N₁=N₂=2⁷ which results in an m factor of 36.6. The numerical solution using techniques described herein is calculated from:

u _ij,t+Δt ^N=χ_ij[u_ij,t^N+Δt(_γu|_ij,t+r_ijl,t)]+(1−χ_ij)w_ij,t+Δt  (60)

where

_γ u| _ij,t=0.2Δx₂Δx₁{u_ij,t^NFFT⁻¹[FFT(ω_ij,t^N)FFT((u_ij,t^N)²)]−(u_ij,t^N)²FFT⁻¹[FFT(ω_ij,t^N)FFT(u_ij,t^N)]}  (61)

EQ. (61) was obtained when substituting f and h with U in EQ. (34). EQ. (61) was computed to solve this nonlinear transient problem up to t=15 with Δt=0.01. Simulation results for three snapshots at t=1, 5 and 15 are shown in the left column of FIG. 11. The relative error (shown in the right column of FIG. 11) was computed from:

$\begin{array}{cc}\mathrm{error}\left(x_{\mathrm{ij}},t\right)=\frac{❘u_{\mathrm{ij},t}^N-u_{\mathrm{ij},t}^{\mathrm{exact}}❘}{{u_{\mathrm{ij},0}}_{\infty }}& \left(62\right)\end{array}$

This example shows that mechanisms described herein are capable of solving nonlinear volume-constrained PD diffusion problems in higher dimensions.

FIG. 12 shows examples of convergence studies in terms of spatial and temporal resolutions for the nonlinear 2D volume-constrained peridynamic transient diffusion problem determined using mechanisms described herein for efficient peridynamic modeling with bounded domains.

A spatial and temporal convergence study for the nonlinear 2D example described above in connection with FIG. 11 was also performed. For spatial convergence, a fixed Δt=10⁻⁵ (for 0<t≤5) was used, and m-factor was varied from 4.6 to 73.1 (N₁=N₂ from 2⁴ to 2⁸). For temporal convergence, a fixed spatial resolution of N₁=N₂=2⁸ (m=73.1) was used, and Δt was varied from 0.1 to 10⁻⁵ (for 0<t≤5). The results for these tests are shown in FIG. 12, with the results for spatial convergence shown in FIG. 12, panel (a), and the results for temporal convergence shown in FIG. 12, panel (b). Here, the relative error was computed by:

$\underset{0<t<5}{\max }\frac{{u_{\mathrm{ij},t}-u_{\mathrm{ij},t}^N}_{\infty }}{{u_{\mathrm{ij},t=0}}_{\infty }}.$

Similar to the 1D example, the convergence rate is quadratic in terms of Δx and linear in terms of Δt (due to the Forward Euler time integrator) for this 2D nonlinear example.

FIG. 13 shows example schematics of a cube with an insulated cut-out, and an extended domain with nonlocal boundary layers for imposing boundary conditions using mechanisms described herein for efficient peridynamic modeling with bounded domains.

In this example, mechanisms described herein and the meshfree (one-point Gaussian integration) spatial discretization were used to solve a 3D PD transient diffusion example with non-periodic local boundary conditions in a domain with a non-trivial geometry. The performance of the two approaches was compared for discretization sizes ranging from 2 million (m-factor of 5.8) to more than 1 billion nodes (m-factor of 46.5), in 3D. Note that in practical problems, there is rarely a reason for using m-factor values larger than 5-10. However, when high accuracy is required, or for problems where nonlocality is relatively large compared to the sample's size, larger m-values (such as the values used for this example) may be needed. The quadratic convergence rate for the spatial discretization discussed in connection with FIGS. 8 and 12 motivates using such large m-values.

The example in this section is a 3D PD diffusion in a 2×2×2 cube with a cutout, having fixed (but different) U values imposed on the top and bottom surfaces (Dirichlet boundary conditions), and zero flux conditions imposed on all four sides of the cube (Neumann boundary conditions). As shown in FIG. 13, the cube has a horizontally aligned disk-shaped insulated cutout/obstacle at the center. The cutout acts as a diffusion barrier with zero flux boundary conditions on all of its surfaces. The mathematical layout for this problem is:

$\begin{array}{cc}\frac{\partial u\left(x,t\right)}{\partial t}=v\mu \left(x^{\prime }-x\right)\left[u\left(x^{\prime }\right)-u\left(x\right)\right]{\mathrm{dx}}^{\prime }& \left(63\right)\end{array}$

on Ω=[−1,1]³ with the kernel function:

$\begin{array}{cc}\mu \left(x\right)=\{\begin{array}{cc}\frac{9}{2\pi {\delta }^3{x}^2};& x\le \delta \\ 0;& x>\delta \end{array},& \left(64\right)\end{array}$

subjected to the initial condition:

u(x,0)=0,  (65)

and the local boundary conditions:

$\begin{array}{cc}\begin{array}{cc}u\left(x,t\right)=1,x=\left\{x_1,x_2,1\right\}& \left(a\right)\\ u\left(x,t\right)=0,x=\left\{x_1,x_2,-1\right\}& \left(b\right)\\ \nabla u\left(x,t\right)·n_1=0,x=\left\{1,x_2,x_3\right\}\mathrm{and}{n}_1=\left\{1,0,0\right\}& \left(c\right)\\ \nabla u\left(x,t\right)·n_2=0,x=\left\{-1,x_2,x_3\right\}\mathrm{and}{n}_2=\left\{-1,0,0\right\}& \left(d\right)\\ \nabla u\left(x,t\right)·n_3=0,x=\left\{x_1,1,x_3\right\}\mathrm{and}{n}_3=\left\{0,1,0\right\}& \left(e\right)\\ \nabla u\left(x,t\right)·n_4=0,x=\left\{x_1,-1,x_3\right\}\mathrm{and}{n}_4=\left\{0,-1,0\right\}& \left(f\right)\\ \nabla u\left(x,t\right)·n_5=0,x=\left\{x_1,x_2,0.1\right\}\mathrm{and}\sqrt{x_1^2+x_2^2}\le 0.7,\mathrm{and}{n}_5=\left\{0,0,1\right\}& \left(g\right)\\ \nabla u\left(x,t\right)·n_6=0,x=\left\{x_1,x_2,-0.1\right\}\mathrm{and}\sqrt{x_1^2+x_2^2}\le 0.7,\mathrm{and}{n}_6=\left\{0,0,-1\right\}& \left(h\right)\\ \nabla u\left(x,t\right)·n_7=0,\sqrt{x_1^2+x_2^2}=0.7\mathrm{and}-0.1\le x_3\le 0.1,\mathrm{and}{n}_7=\left\{\frac{x_1}{0.7},\frac{x_2}{0.7},0\right\}& \left(i\right)\end{array}& \left(66\right)\end{array}$

EQS. (66)(a) and (b) are the top and bottom Dirichlet boundary conditions, (c) to (f) are zero flux conditions on the cube sides, and (g), (h), and (i) are the zero-flux boundary conditions on the cutout surfaces. For this example, the cutout radius is 0.7, and δ=0.1. Similar to the 1D example in FIG. 7, in order to impose local boundary conditions, the fictitious domain and the formulas described in Appendix A were used. For the external surfaces, the cube was extended by δ in all six directions, and volume constraints were defined on these extensions. For the top and bottom (with conditions of EQS. (66) (g) and (h)) surfaces of the cutout, the extension ends up being the actual cutout volume, of thickness 2δ. FIG. 13 shows the original and the domain extended by the constrained volume: Ω∪Γ.

Since the extended domain (Ω∪Γ) is already a cube, similar to the 1D example in FIG. 7, the periodic box can be defined to be Ω∪Γ without any gap (i.e., Λ=Ø).

Note that here, the boundary condition in EQ. (66)(i) in the PD model was not enforced, because, given the non-convex characteristics of the geometry, it would overlap regions of fictitious nodes associated with the boundary conditions of EQS. (66) (g), (h), and (i). That would create difficulties in the Fourier space because different points in the domain would be linked to values at the same fictitious node. The cutout thickness of 2δ was chosen to avoid this complication for the domain extension corresponding to the top and bottom surfaces of the cutout. Note that this issue may not appear when using a different method to enforce the local boundary conditions. Here, for simplicity the fictitious nodes method was used. However, as described below in connection with FIG. 15, the error introduced by not enforcing the condition in EQ. (66) (i) is relatively small. Note that this issue is only present when enforcing a local boundary condition in a PD formulation; it does not exist for volume constraints (nonlocal boundary conditions).

FIG. 14 shows examples of 3D peridynamic diffusion simulation results for a cube with an insulated cutout determined using mechanisms described herein for efficient peridynamic modeling with bounded domains.

For discretization, N₁=N₂=N₃=2⁷ was used first, resulting in m-factor of 5.8 and N=2²¹ as the total number of nodes. The time step used was Δt=5.5×10⁻⁵ to satisfy the stability condition in EQ. (46). A general algorithm for 3D diffusion using mechanisms described herein is described in Appendix C, which is hereby incorporated by reference herein. FIG. 14 shows simulation results in three snapshots at t=0.1, 0.3 and 1. Two perpendicular cross-sectional views of the 3D domain are plotted to make the insulated discontinuity (cutout) visible in the figure.

FIG. 15 shows examples comparing the 3D solution in a cross-section determined using mechanisms described herein for efficient peridynamic modeling with bounded domains, and results generated using conventional techniques.

The temperature/concentration distributions in FIG. 14 shows that the boundary conditions were enforced and the insulated cutout acted as an obstacle to the diffusive transport process. The solution generating using mechanisms described herein for the peridynamic problem was compared with a finite element solution for the corresponding classical boundary value diffusion problem. For the finite element analysis, Abaqus/Standard 6.14-2 (note that this is an implicit solver, compared with the mechanisms described herein which can provide an explicit solver) was used with over 10⁷ DC3D4 type elements and more than 1.8×10⁶ nodes, which is about the same number of degrees of freedom used in the simulation performed in accordance with mechanisms described herein.

As shown in FIG. 15, the nonlocal solution generated using mechanisms described herein is very close to the finite element local solution generated with Abaqus. Although contour levels are set to be the same for two top legends, the contour colors are slightly different since the corresponding colors from Abaqus graphics are slightly different from those in Tecplot used in plotting the results obtained using mechanisms described herein. The bottom row in FIG. 15 shows the absolute relative error distribution, along with its mean value. Relative error here was computed as

$\frac{❘u_{\mathrm{ij},t}^{\mathrm{FCBM}}-u_{\mathrm{ij},t}^{\mathrm{Abaqus}}❘}{{u_{\mathrm{ij},0}^{\mathrm{Abaqus}}}_{\infty }}$

on the displayed cross-sectional plane.

Note that there is a small discrepancy between the two solutions near the sides of the rectangular obstacle. This is likely due to not enforcing the exact insulated condition on the short lateral sides of the cutout. For the Abaqus solution, all boundary conditions in EQ. (66) were enforced.

In order to demonstrate the superior efficiency that can be achieved using mechanisms described herein, the same example was solved via the conventional meshfree method with one-point Gaussian quadrature, and the computational times for these approaches was compared for various spatial resolutions. Both approaches were implemented in MATLAB R2014b. Simulations were performed on a supercomputer machine in Holland Computing Center of the University of Nebraska-Lincoln, with Intel Xeon E5-2670 2.60 GHz CPUs, up to 512 GB RAM per CPU and a Tesla V100 GPU with 32 GB memory. MATLAB's built-in parallel FFT solvers use multithreading and GPU computations, and this was used in the computational tests. TABLE 1 shows the computational time for the 3D example above using the meshfree with Gaussian quadrature on a single CPU versus using mechanisms described herein on a single CPU, on 8 or 16 CPUs, and on a GPU, for four spatial discretizations. The time steps are identical for all simulations since the stability condition is independent of the discretization size (it only depends on ν, μ and δ, which are the same for all tests). TABLE 1 shows computational time for solving the 3D example up to t=1, using mechanisms described herein (labeled FCBM-EC) and the meshfree method with Gaussian quadrature (GQ).

	TABLE 1
	m-factor
	5.8	11.6	23.3	46.5
	Number of nodes
	221 =	224 =	227 =	230 =
	2,097,152	16,777,216	134,217,728	1,073,741,824
Meshfree with GQ (1 CPU)	4.6	hrs	12	days	2.1	yrs (!)	1.4	centuries (!)
FCBM-EC (1 CPU)	10.7	min	1.9	hrs	13.7	hrs	5.1	days
FCBM-EC (8 CPUs)	6	min	36.6	min	4.5	hrs	1.4	days
FCBM-EC (16 CPUs)	2.7	min	34.7	min	3.9	hrs	1.3	days
FCBM-EC (GPU)	29	min	38.3	min	2.1	hrs	out of memory

As observed the diffusion example with over 1 billion nodes is solved in a few days using mechanisms described herein, while the same computation is intractable with the meshfree method using one-point Gaussian quadrature.

Note that the computational time for the meshfree method in the cases with m=11.6, 23.3, and 46.5 (N=2²⁴, 2²⁷, and 2³⁰) are estimated based on the assumption that it scales as O(N²). Although the computational time for the meshfree method in the case with m=11.6 (N=2²⁴) would be reasonable (several days), the available memory was not sufficient to perform this test. The simulations performed using mechanisms described herein, in contrast, did not run out of memory even for 1 billion nodes. A reason is likely due to peridynamic meshfree solvers storing family nodes for each node at the beginning to increase efficiency when doing the Gaussian quadrature for the convolution integral. However, there is no need to store family nodes when using mechanisms described herein, as the convolution integrals are replaced by multiplication of the Fourier coefficients, which do not require the family nodes for each node. Note that while the meshfree method can also be performed without storing nodal families, the computational time increases substantially when one needs to re-compute them at every time step.

The results in TABLE 1 show efficiency gains that can be realized using mechanisms described herein versus the meshfree discretization. Notice also the correlation of the gains with the m-factor value. The reason for the larger gains at higher m-values is related to the same point mentioned above: the convolution integrals are replaced by multiplication of the Fourier coefficients, which do not require the family nodes for each node. In other words, a nonlocal computation in the physical space can be transformed into a local computation in Fourier space.

FIG. 16 shows examples of efficiency gains relative to m-factor values for the 3D solution determined using mechanisms described herein for efficient peridynamic modeling with bounded domains over the meshfree technique, and the speedup versus problem size (N) using mechanisms described with parallel fast Fourier transform solvers.

Efficiency gains (computed using EQ. (67)) and speedup (computed using EQ. (68)) are plotted versus the problem size (degrees of freedom) and m-factor values in FIG. 13.

$\begin{array}{cc}\mathrm{Efficiency}\mathrm{gain}=\frac{\mathrm{computation}\mathrm{time}\mathrm{for}\mathrm{FCBM}}{\mathrm{computation}\mathrm{time}\mathrm{for}\mathrm{meshfree}\mathrm{method}\mathrm{with}\mathrm{one}\mathrm{CPU}}& \left(67\right)\end{array}$ $\begin{array}{cc}\mathrm{Speedup}=\frac{\mathrm{computation}\mathrm{time}\mathrm{for}\mathrm{FCBM}\mathrm{with}\mathrm{multiple}\mathrm{CPUs}\mathrm{or}\mathrm{GPU}}{\mathrm{computation}\mathrm{time}\mathrm{for}\mathrm{FCBM}\mathrm{with}\mathrm{one}\mathrm{CPU}}& \left(68\right)\end{array}$

These computational experiments (via Matlab's parallel FFT functions) show that little was gained from using more than 8 CPUs with multi-threading. They also show that GPU computations using Matlab's intrinsic functions was superior to Matlab's multi-threading for problems with higher spatial resolution (here when the N>2²⁴≅16.8×10⁶), as long as GPU computations are feasible from the memory point of view.

The FE and FCBM-EC simulations shown in FIG. 15 were performed on the same machine using 16 CPUs, and both have a similar number of degrees of freedom (about 2 million). The Abaqus simulation took about 2.5 hours using an implicit solver (note that there is no explicit solver in Abaqus for diffusion problems) to complete while the simulation using mechanisms described herein took only 2.7 minutes. Note that general-purpose software, such as the diffusion solver in Abaqus, is not necessarily optimized for the particular problem setup here. Additionally, an explicit-type solver was used for the PD results, while Abaqus used an implicit solver. This is not necessarily a disadvantage for Abaqus, as an implicit solver implemented for mechanisms described herein can be expected to lead to even more efficient performance. Additionally, the implementation of the mechanisms described herein used to perform the simulation is not necessarily an optimal one, as Matlab's libraries of functions were used.

Further Examples Having a Variety of Features

Implementation examples are described in the following numbered clauses:

1. A method comprising: receiving a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, wherein the plurality of parameters includes a shape of a domain Ω, a horizon size δ, and one or more boundary constraints; drawing a periodic box T that includes the domain Ω and a boundary layer Γ extending from the domain Ω; imposing the one or more boundary constraints on the boundary layer Γ; performing a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and causing a result of the simulation to be presented.

2. The method of clause 1, wherein the non-local model comprises a peridynamic model.

3. The method of any one of clauses 1 or 2, wherein the plurality of parameters includes one or more initial conditions.

4. The method of any one of clauses 1 to 3, wherein performing the simulation comprises: discretizing the non-local model using one or more Fourier transforms, wherein integral operators associated with the model are expressed in terms of convolutions; and solving the discretized non-local model on T.

5. The method of any one of clauses 1 to 4, wherein the non-local model comprises a model of diffusion.

6. The method of clause 5, wherein the model of diffusion is non-linear.

7. The method of any one of clauses 1 to 4, wherein the non-local model comprises a model of deformation and fracture.

8. The method of clause 7, wherein the model of deformation and fracture is non-linear.

9. The method of any one of clauses 1 to 4, wherein the non-local model comprises a model of dissolution damage.

10. The method of any one of clauses 1 to 4, wherein the non-local model comprises a model of corrosion damage.

11. The method of any one of clauses 1 to 10, wherein the domain Ω comprises a three-dimensional (3D) shape.

12. The method of any one of clauses 1 to 11, further comprising, receiving a time step size Δt, wherein the time step size satisfies a constraint.

13. The method of clause 12, wherein the time step size Δt satisfies Δt≤1/vβ, where v is a scalar diffusivity constant associated with the non-local model, and β corresponds to an integral of a kernel function associated with the non-local model over a region H_0 of the domain Ω, H_0 having a size based on the horizon size δ and the dimensionality of the domain Ω.

14. The method of any one of clauses 1 to 13, further comprising, receiving a number of nodes N, in each dimension.

15. The method of any one of clauses 1 to 14, wherein causing the result of the simulation to be presented comprises: causing a visual representation of at least a portion of the non-local model at a particular time step to be presented.

16. The method of any one of clauses 1 to 12 or 14, wherein causing the result of the simulation to be presented comprises: causing a visual state of at least a portion of the non-local model at a particular load step to be presented.

17. The method of any one of clauses 1 to 16, wherein receiving the plurality of parameters associated with the peridynamic model to be simulated, comprises: receiving the plurality of parameters from a remote computing device over a wide area network (WAN).

18. A system for peridynamic modeling, the system comprising: at least one processor that is configured to: receive a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, wherein the plurality of parameters includes a shape of a domain Ω, a horizon size δ, and one or more boundary constraints; draw a periodic box T that includes the domain Ω and a boundary layer Γ extending from the domain Ω; impose the one or more boundary constraints on the boundary layer Γ; perform a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and cause a result of the simulation to be presented.

19. A non-transitory computer readable medium containing computer executable instructions that, when executed by a processor, cause the processor to perform a method for peridynamic modeling, the method comprising: receiving a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, including a shape of a domain Ω, a horizon size δ, and one or more boundary constraints; drawing a periodic box T that includes the domain Ω and a boundary layer Γ extending from the domain Ω; imposing the one or more boundary constraints on the boundary layer Γ; performing a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and causing a result of the simulation to be presented.

20. A system comprising: at least one hardware processor that is configured to: perform a method of any one of clauses 1 to 17.

21. A non-transitory computer readable medium containing computer executable instructions that, when executed by a processor, cause the processor to perform a method of any one of clauses 1 to 17.

In some embodiments, any suitable computer readable media can be used for storing instructions for performing the functions and/or processes described herein. For example, in some embodiments, computer readable media can be transitory or non-transitory. For example, non-transitory computer readable media can include media such as magnetic media (such as hard disks, floppy disks, etc.), optical media (such as compact discs, digital video discs, Blu-ray discs, etc.), semiconductor media (such as RAM, Flash memory, electrically programmable read only memory (EPROM), electrically erasable programmable read only memory (EEPROM), etc.), any suitable media that is not fleeting or devoid of any semblance of permanence during transmission, and/or any suitable tangible media. As another example, transitory computer readable media can include signals on networks, in wires, conductors, optical fibers, circuits, or any suitable media that is fleeting and devoid of any semblance of permanence during transmission, and/or any suitable intangible media.

It should be noted that, as used herein, the term mechanism can encompass hardware, software, firmware, or any suitable combination thereof.

It should be understood that the above-described steps of the processes of FIG. 6 can be executed or performed in any order or sequence not limited to the order and sequence shown and described in the figures. Also, some of the above steps of the processes of FIG. 6 can be executed or performed substantially simultaneously where appropriate or in parallel to reduce latency and processing times.

Although the invention has been described and illustrated in the foregoing illustrative embodiments, it is understood that the present disclosure has been made only by way of example, and that numerous changes in the details of implementation of the invention can be made without departing from the spirit and scope of the invention, which is limited only by the claims that follow. Features of the disclosed embodiments can be combined and rearranged in various ways.

Claims

1. A method comprising: receiving a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, wherein the plurality of parameters includes a shape of a domain Ω, a horizon size δ, and one or more boundary constraints;drawing a periodic box that includes the domain Ω and a boundary layer Γ extending from the domain Ω;imposing the one or more boundary constraints on the boundary layer Γ;performing a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; andcausing a result of the simulation to be presented.

2. The method of claim 1, wherein the non-local model comprises a peridynamic model.

3. The method of claim 1, wherein the plurality of parameters includes one or more initial conditions.

4. The method of claim 1, wherein performing the simulation comprises: discretizing the non-local model using one or more Fourier transforms, wherein integral operators associated with the model are expressed in terms of convolutions; andsolving the discretized non-local model on .

5. The method of claim 1, wherein the non-local model comprises a model of diffusion.

6. The method of claim 5, wherein the model of diffusion is non-linear.

7. The method of claim 1, wherein the non-local model comprises a model of deformation and fracture.

8. The method of claim 7, wherein the model of deformation and fracture is non-linear.

9. The method of claim 1, wherein the non-local model comprises a model of dissolution damage.

10. The method of claim 1, wherein the non-local model comprises a model of corrosion damage.

11. The method of claim 1, wherein the domain Ω comprises a three-dimensional (3D) shape.

12. The method of claim 1, further comprising, receiving a time step size Δt, wherein the time step size satisfies a constraint.

13. The method of claim 12, wherein the time step size Δt satisfies

14. The method of claim 1, further comprising, receiving a number of nodes N, in each dimension.

15. The method of claim 1, wherein causing the result of the simulation to be presented comprises: causing a visual representation of at least a portion of the non-local model at a particular time step to be presented.

16. The method of claim 1, wherein causing the result of the simulation to be presented comprises: causing a visual state of at least a portion of the non-local model at a particular load step to be presented.

17. The method of claim 1, wherein receiving the plurality of parameters associated with the peridynamic model to be simulated, comprises: receiving the plurality of parameters from a remote computing device over a wide area network (WAN).

18. A system for peridynamic modeling, the system comprising: at least one processor that is configured to: receive a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, wherein the plurality of parameters includes a shape of a domain Ω, a horizon size δ, and one or more boundary constraints;draw a periodic box that includes the domain Ω and a boundary layer Γ extending from the domain Ω;impose the one or more boundary constraints on the boundary layer Γ;perform a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; andcause a result of the simulation to be presented.

19. A non-transitory computer readable medium containing computer executable instructions that, when executed by a processor, cause the processor to perform a method for peridynamic modeling, the method comprising: receiving a plurality of parameters associated with a non-local model that has a convolutional structure to be used for one or more simulations, including a shape of a domain Ω, a horizon size δ, and one or more boundary constraints;drawing a periodic box that includes the domain Ω and a boundary layer Γ extending from the domain Ω;imposing the one or more boundary constraints on the boundary layer Γ;performing a simulation of a physical process using the non-local model and convolutional properties of Fourier transforms; and causing a result of the simulation to be presented.