A fast collocation method for a static bond-based linear peridynamic model

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303 www.elsevier.com/locate/cma

A fast collocation method for a static bond-based linear peridynamic model Xuhao Zhang a , Hong Wang b,∗ a School of Mathematics, Shandong University, Jinan, Shandong 250100, China b Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA

Received 9 October 2015; received in revised form 28 June 2016; accepted 17 August 2016 Available online 25 August 2016

Abstract We develop a fast collocation method for a steady-state bond-based linear peridynamic model in two space dimensions. The method reduces the computational work from O(N 2 ) per Krylov subspace iteration in a traditional collocation method to O(N log N ) and the memory requirement from O(N 2 ) in a collocation method to O(N ), where N is the number of unknowns in the discrete system. Furthermore, the method reduces the computational work of evaluating and assembling the stiffness matrix, which often constitutes a major portion of the overall computational work, from O(N 2 ) in a collocation method to O(N ). The significant reduction of the CPU times and storage in the fast collocation method is achieved by carefully exploring the structure of the stiffness matrix of the collocation method without any lossy compression involved. In other words, the fast method is evaluated in an equivalent but more efficient manner. Therefore, the fast method generates identical numerical solutions as the traditional collocation method does and naturally inherits the stability and convergence properties that were already proved for a traditional collocation method. Numerical results are presented to show the utility of the fast method. c 2016 Elsevier B.V. All rights reserved. ⃝

Keywords: Fast collocation method; Nonlocal diffusion model; Peridynamic model

1. Introduction The classical theory of continuum solid mechanics assumes that all internal forces act effectively through a zero distance and yields mathematical models that are expressed in terms of partial differential equations. These models have difficulties to describe problems with evolving discontinuities, due to its differentiability assumption on displacement fields. The peridynamic models [1,2], as a reformulation of classical continuum solid mechanics, yield nonlocal mathematical formulations that are based on long-range interactions. Constitutive models in the peridynamics

∗ Corresponding author. Fax: +1 803 777 6527.

E-mail address: [email protected] (H. Wang). http://dx.doi.org/10.1016/j.cma.2016.08.020 c 2016 Elsevier B.V. All rights reserved. 0045-7825/⃝

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

281

depend on finite deformation vectors, instead of deformation gradients in classical constitutive models. Consequently, the peridynamic models are natural for the representation of discontinuities in displacement fields and the description of cracks and their evolution in materials. To date, the peridynamic theory and related nonlocal models have been successfully applied in many important applications, including failure and damage in composite laminates [3–6], crack propagation and branching [7–9], crack nucleation [10], phase transformations in solids [11], impact damage [12–16], polycrystal fracture [17], crystal plasticity [18], damage in concrete [19], and geomaterial fragmentation [20]. Mathematically, the peridynamic models give rise to integro-differential equations with singular kernels defined on circular material horizons, where spatial integration is employed to compute the contribution of internal forces in a body to the material response. As the governing equations in peridynamics are continuum models, different numerical approximation techniques can be used to discretize these equations which differ in computational expense, memory requirements, and implementational effort as well as accuracy, convergence, and stability of the numerical solutions [21]. For example, collocation methods and meshfree methods apply directly to the strong form of the governing equations in the peridynamic models, and are relatively simple to implement [22–24,15,25,26]. On the other hand, finite element discretizations of the governing equations in the peridynamic models are based on Galerkin weak forms, and enjoy high convergence rates and are suitable for a wider range of singular kernels [27,28,5,25,29,30]. However, finite element methods double the number of spatial dimensions that need to be discretized, which means a significantly increased computational work and memory requirement. Numerical methods for the peridynamic models present new numerical difficulties that were not encountered in those for the classical continuum solid mechanical models: (i) The accuracy and convergence of the numerical methods for the peridynamic models depends heavily on the accurate evaluation of the singular integrals, which are usually defined on a disk, in the entries of the stiffness matrices. Consequently, different numerical quadratures have been developed to effectively evaluate these singular integrals, which often constitute a major portion of the overall computational work [12,13,31–33,22,15,34]. (ii) Because of their nonlocality, numerical methods for the peridynamic models usually generate dense stiffness matrices in which the bandwidths increase to infinity as the mesh size decreases to zero. Direct solvers are widely used in the peridynamic modeling, which have O(N 2 ) memory requirement to store the stiffness matrix and O(N 3 ) computations to find the numerical solutions. In this paper we develop fast numerical methods for the peridynamic models via an alternative approach. Namely, we do not intend to develop a more accurate numerical discretization or numerical quadrature for these models. Rather, we explore the structure of the stiffness matrices of existing numerical methods in order to significantly reduce the computational work to evaluate and assemble the stiffness matrix and to compute the numerical solutions, as well as the memory requirement to store the stiffness matrix. The idea applies to different numerical methods and numerical quadratures in the literature. For definiteness, in this paper we develop a fast collocation method for a static bond-based linear peridynamic model in two space dimensions. The rest of the paper is organized as follows. In Section 2 we present a piecewise bilinear collocation method for model problem (1). In Section 3 we analyze the structure of the stiffness matrix and based on the properties of the displacement vector and the stiffness matrix we derive an appropriate decomposition. In Section 4 we prove the appropriately decomposed submatrices of the stiffness matrix are block-Toeplitz–Toeplitz-block (BTTB). In Section 5 we develop an efficient storage mechanism for the stiffness matrix that reduces the storage from O(N 2 ) to O(N ) and a fast Krylov subspace iterative method that reduces each iteration from O(N 2 ) computations to O(N log N ) computations. In Section 6 we carry out numerical experiments to show the utility of the fast method. In Section 7 we draw concluding remarks. 2. A static bond-based linear peridynamic model and its collocation discretization We consider the following static bond-based linear peridynamic model in two space dimensions [1,35,2]  C(x′ , x)(u(x′ ) − u(x))dx′ = b(x), x ∈ Ω , Bδ (x)

u(x) = g(x), x ∈ Ωc .

(1)

282

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

Here C is the micromodulus tensor of the form  (x′ − x) (x′ − x)  C(x′ , x) = σ (∥x′ − x∥) ⊗ . ∥x′ − x∥ ∥x′ − x∥

(2)

In (1) and (2), x := (x, y)T , Ω = (0, x R ) × (0, y R ) is a rectangular domain, the material horizon Bδ (x) is assumed to be a closed disk of radius δ > 0 centered at x, where δ is referred to as the peridynamic horizon. Ωc denotes a boundary zone surrounding Ω of width δ. σ (∥x∥) is the kernel function of the integral operator. u(x) = [v(x), w(x)]T  T represents the displacement of the material, b(x) = bv (x), bw (x) represents the external force density, and  T g(x) = g v (x), g w (x) is the prescribed nonlocal boundary data imposed on the domain Ωc . Let the positive integers N x and N y denote the numbers of intervals in the discretized domain in the x and y directions, respectively. We define a spatial partition in Ω = [0, x R ] × [0, y R ] by xi := i h x for i = 0, 1, . . . , N x and y j := j h y for j = 0, 1, . . . , N y , where h x := x R /N x and h y := y R /N y are the mesh sizes in the x and y directions. To handle the discretization of the boundary zone Ωc , we extend the partition to the nodes (xi , y j ) for i = −K + 1, . . . , −1, 0, 1, . . . , N x − 1, N x , N x + 1, . . . , N x + K − 1 and j = −L + 1, . . . , −1, 0, 1, . . . , N y − 1, N y , N y + 1, . . . , N y + L − 1. Here K := ⌈δ/ h x ⌉,

(3)

L := ⌈δ/ h y ⌉

are the ceilings of δ/ h x and δ/ h y , respectively. Let ψ(ξ ) = 1 − |ξ | for ξ ∈ [−1, 1] or zero otherwise. The two dimensional pyramid functions φi j (x, y) centered at (xi , y j ) can be expressed as x − x  y − y  j i ψ , (x, y) ∈ Ω ∪ Ωc , hx hy −K + 1 6 i 6 N x + K − 1, −L + 1 6 j 6 N y + L − 1. φi j (x, y) = ψ

(4)

Then the x-component v(x, y) and the y-component w(x, y) of the displacement can be expressed as v(x, y) = w(x, y) =

+L−1 Nx +K −1 N y i ′ =−K +1 j ′ =−L+1 +L−1 Nx +K −1 N y

vi ′ , j ′ φi ′ , j ′ (x, y), (5) wi ′ , j ′ φi ′ , j ′ (x, y).

i ′ =−K +1 j ′ =−L+1

For −K + 1 ≤ i ′ ≤ 0, N x ≤ i ′ ≤ N x + K − 1 and −L + 1 ≤ j ′ ≤ 0, N y ≤ j ′ ≤ N y + L − 1, vi ′ , j ′ = g v (xi ′ , y j ′ ) and wi ′ , j ′ = g w (xi ′ , y j ′ ) are specified in (5). Hence, we insert (5) into (1) and enforce the resulting equations at the collocation points (xi , y j ) for i = 1, 2, . . . , N x − 1 and j = 1, 2, . . . , N y − 1 to obtain the following scheme  Bδ (xi ,y j )

+L−1 Nx +K −1 N y   σ (∥(x ′ − xi , y ′ − y j )∥)  ′ 2 ′ ′ ′ , j ′ φi ′ , j ′ (x , y ) (x − x ) v − v i i, j i ((x ′ − xi )2 + (y ′ − y j )2 ) i ′ =−K +1 j ′ =−L+1

 + (x − xi )(y − y j ) wi, j − ′

′

+L−1 Nx +K −1 N y

 wi ′ , j ′ φi ′ , j ′ (x ′ , y ′ ) d x ′ dy ′ = bv (xi , y j ),

i ′ =−K +1 j ′ =−L+1

 Bδ (xi ,y j )

+L−1 Nx +K −1 N y   σ (∥(x ′ − xi , y ′ − y j )∥)  ′ 2 (y − y ) w − wi ′ , j ′ φi ′ , j ′ (x ′ , y ′ ) i, j j ′ 2 ′ 2 ((x − xi ) + (y − y j ) ) i ′ =−K +1 j ′ =−L+1

 + (x ′ − xi )(y ′ − y j ) vi, j −

+L−1 Nx +K −1 N y i ′ =−K +1 j ′ =−L+1

 vi ′ , j ′ φi ′ , j ′ (x ′ , y ′ ) d x ′ dy ′ = bw (xi , y j ).

(6)

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

283

The collocation scheme (6) can be written in a matrix form A2N u2N = b2N .

(7)

Here N := (N x − 1)(N y − 1) refers to the number of grid points at which the values of v and w need to be computed. The 2N -dimensional unknown displacement vector u2N and the right-hand side loading vector b2N of (6) are labeled in their natural order. At each node (xi , y j ), the x component vi, j or bi,v j is put first and the y component wi, j or bi,w j second:  u2N := v1,1 , w1,1 , . . . , v Nx −1,1 , w Nx −1,1 , v1,2 , w1,2 , . . . , v Nx −1,2 , w Nx −1,2 , . . . , v1,N y −1 , w1,N y −1 , . . . , v Nx −1,N y −1 , w Nx −1,N y −1

T

,

(8)

v w v b2N := [b1,1 , b1,1 , . . . , bvNx −1,1 , bw N x −1,1 , b1,2 , w b1,2 , . . . , bvNx −1,2 , v bw N x −1,2 , . . . , b1,N y −1 , T w b1,N , . . . , bvNx −1,N y −1 , bw N x −1,N y −1 ] . y −1

The 2N -by-2N stiffness matrix A2N can be written as the following N -by-N block matrix

A2N

 A 1,1  A2,1  . =  ..  A N −1,1 A N ,1

A1,2 A2,2 .. . A N −1,2 A N ,2

... ... .. . ... ...

A1,N −1 A2,N −1 .. .

A1,N A2,N .. .

A N −1,N −1 A N ,N −1

A N −1,N A N ,N

   ,  

(9)

where each entry of A2N is a 2-by-2 matrix  Am,n =

Av,v m,n

Av,w m,n

Aw,v m,n

Aw,w m,n

 ,

(10)

with the entries in Am,n being defined by

Bδ (xi ,y j )

σ (∥(x ′ − xi , y ′ − y j )∥) ′ (x − xi )2 (δm,n − φi ′ , j ′ (x ′ , y ′ ))d x ′ dy ′ , ((x ′ − xi )2 + (y ′ − y j )2 )

Bδ (xi ,y j )

σ (∥(x ′ − xi , y ′ − y j )∥) ′ (x − xi )(y ′ − y j )(δm,n − φi ′ , j ′ (x ′ , y ′ ))d x ′ dy ′ , ((x ′ − xi )2 + (y ′ − y j )2 )

Av,v m,n :=



Av,w m,n :=



Aw,w m,n

 := Bδ (xi ,y j )

(11)

σ (∥(x ′ − xi , y ′ − y j )∥) ′ (y − y j )2 (δm,n − φi ′ , j ′ (x ′ , y ′ ))d x ′ dy ′ . ((x ′ − xi )2 + (y ′ − y j )2 )

v,w Here δm,n = 1 for m = n or 0 otherwise, and Aw,v m,n = Am,n . The global indices m and n are related to the indices ′ ′ (i, j) and (i , j ) by

m = ( j − 1) × (N x − 1) + i, 1 6 i 6 N x − 1, 1 6 j 6 N y − 1, n = ( j ′ − 1) × (N x − 1) + i ′ , 1 6 i ′ 6 N x − 1, 1 6 j ′ 6 N y − 1.

(12)

284

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

The entries of the right-hand side vector b2N are given by    v v bi, j = b (xi , y j ) + −K +1≤i ′′ ≤0 −L+1≤ j ′′ ≤0 N x ≤i ′′ ≤N x +K −1 N y ≤ j ′′ ≤N y +L−1

Bδ (xi ,y j )

σ (∥(x ′ − xi , y ′ − y j )∥) ((x ′ − xi )2 + (y ′ − y j )2 )

× (x ′ − xi )2 g v (xi ′′ , y j ′′ )φi ′′ , j ′′ (x ′ , y ′ )d x ′ dy ′    σ (∥(x ′ − xi , y ′ − y j )∥) + ′ 2 ′ 2 Bδ (xi ,y j ) ((x − x i ) + (y − y j ) ) −K +1≤i ′′ ≤0 −L+1≤ j ′′ ≤0 N x ≤i ′′ ≤N x +K −1 N y ≤ j ′′ ≤N y +L−1

bi,w j

× (x ′ − xi )(y ′ − y j )g w (xi ′′ , y j ′′ )φi ′′ , j ′′ (x ′ , y ′ )d x ′ dy ′ ,    σ (∥(x ′ − xi , y ′ − y j )∥) = bw (xi , y j ) + ′ 2 ′ 2 Bδ (xi ,y j ) ((x − x i ) + (y − y j ) ) −K +1≤i ′′ ≤0 −L+1≤ j ′′ ≤0

(13)

N x ≤i ′′ ≤N x +K −1 N y ≤ j ′′ ≤N y +L−1

× (x ′ − xi )(y ′ − y j )g v (xi ′′ , y j ′′ )φi ′′ , j ′′ (x ′ , y ′ )d x ′ dy ′    σ (∥(x ′ − xi , y ′ − y j )∥) + ′ 2 ′ 2 Bδ (xi ,y j ) ((x − x i ) + (y − y j ) ) −K +1≤i ′′ ≤0 −L+1≤ j ′′ ≤0 N x ≤i ′′ ≤N x +K −1 N y ≤ j ′′ ≤N y +L−1

× (y ′ − y j )2 g w (xi ′′ , y j ′′ )φi ′′ , j ′′ (x ′ , y ′ )d x ′ dy ′ ,

1 6 i 6 N x − 1, 1 6 j 6 N y − 1.

The indices (i ′′ , j ′′ ) refer to the nodes in Ωc . Because of the nonlocal property of the peridynamic model (1), the stiffness matrix A2N is dense since K = O(δ/ h x ) and L = O(δ/ h y ) tend to infinity as h x and h y tend to zero. Hence, the widely used direct solvers for the linear system (7) require O(N 2 ) memory and O(N 3 ) computational work. Moreover, because of the singularity of the kernel function in the peridynamic model (1), the numerical integration of the nonzero entries of the stiffness matrix A2N is involved and expensive. The evaluation and assembly of the O(N 2 ) nonzero entries of the stiffness matrix A2N often constitutes a major portion of the overall computations. 3. Analysis of A2N u2N for any vector u2N ∈ R2N A conventional evaluation of the matrix–vector multiplication A2N u2N in any Krylov subspace iterative method has a computational complexity of O(N 2 ) and a memory requirement of O(N 2 ), while all other operations have a linear computational complexity. Therefore, to develop a fast Krylov subspace iterative method we need to design an efficient storage mechanism for the stiffness matrix A2N and a fast matrix–vector multiplication algorithm to evaluate A2N u2N for any vector u2N ∈ R2N . Due to the coupling between the x component v and the y component w of the displacement u2N introduced by the tensor operator C(x′ , x) in (2), the stiffness matrix A2N in (9) is not BTTB. To develop a fast method, we analyze the displacement vector u2N and the stiffness matrix A2N . We first decompose the displacement vector u2N as the following two auxiliary vectors, which replace the even or odd entries of u2N by 0, respectively  T uo2N := v1,1 , 0, . . . , v Nx −1,1 , 0, v1,2 , 0, . . . , v Nx −1,2 , 0, . . . , v1,N y −1 , 0, . . . , v Nx −1,N y −1 , 0 ,  T ue2N := 0, w1,1 , . . . , 0, w Nx −1,1 , 0, w1,2 , . . . , 0, w Nx −1,2 , . . . , 0, w1,N y −1 , . . . , 0, w Nx −1,N y −1 .

(14)

We also introduce two N -dimensional vectors v N and w N as follows  T v N := v1,1 , . . . , v Nx −1,1 , v1,2 , . . . , v Nx −1,2 , . . . , v1,N y −1 , . . . , v Nx −1,N y −1 ,  T w N := w1,1 , . . . , w Nx −1,1 , w1,2 , . . . , w Nx −1,2 , . . . , w1,N y −1 , . . . , w Nx −1,N y −1 .

(15)

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

285

It is clear that u2N = uo2N + ue2N .

(16)

Moreover, eliminating all the zero entries in the 2N -dimensional vectors uo2N and ue2N naturally yields the N dimensional vectors v N and w N , respectively. We accordingly introduce the following two 2N -by-2N matrices Ao2N and Ae2N , which are obtained by replacing the even column entries or odd column entries of A2N by zeros, respectively. Namely, Ao1,1  Ao2,1   :=  ...  o A N −1,1 AoN ,1

Ao1,2 Ao2,2 .. . AoN −1,2 AoN ,2



Ao2N

... ... .. . ... ...

Ao1,N −1 Ao2,N −1 .. . AoN −1,N −1 AoN ,N −1

 Ao1,N Ao2,N    .. , .  o A N −1,N  AoN ,N

(17)

where each entry of Ao2N is a 2-by-2 matrix of the form  Aom,n :=

Av,v m,n

 0

Aw,v m,n

0

,

(18)

and Ae1,1  Ae2,1   :=  ...  e A N −1,1 AeN ,1 

Ae2N

Ae1,2 Ae2,2 .. . AeN −1,2 AeN ,2

... ... .. . ... ...

Ae1,N −1 Ae2,N −1 .. . AeN −1,N −1 AeN ,N −1

 Ae1,N Ae2,N    .. , .  e A N −1,N  AeN ,N

(19)

where each entry of Ae2N is a 2-by-2 matrix of the form Aem,n

:=

 0

Av,w m,n

0

Aw,w m,n

 .

(20)

v,w w,v w,w Finally, we introduce the N -by-N matrices Av,v by replacing each 2-by-2 block Am,n in N , A N , A N , and A N v,v v,w w,v the matrix A2N in (9) by the corresponding entry Am,n , Am,n , Am,n , and Aw,w m,n , respectively. That is,

Av,v N

 Av,v 1,1  v,v  A2,1   . :=   ..   Av,v  N −1,1 Av,v N ,1

Av,w N

 Av,w 1,1  v,w A  2,1   . :=   ..   Av,w  N −1,1 Av,w N ,1

Av,v 1,2

...

Av,v 1,N −1

Av,v 2,2

...

Av,v 2,N −1

.. .

..

.. .

Av,v N −1,2

...

Av,v N −1,N −1

Av,v N ,2

...

Av,v N ,N −1

Av,w 1,2

...

Av,w 1,N −1

Av,w 2,2

...

Av,w 2,N −1

.. .

..

.. .

Av,w N −1,2

...

Av,w N −1,N −1

Av,w N ,2

...

Av,w N ,N −1

.

.

 Av,v 1,N  Av,v 2,N    .. , .   v,v A N −1,N  

(21)

Av,v N ,N

 Av,w 1,N  Av,w 2,N    .. , .   v,w A N −1,N   Av,w N ,N

(22)

286

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

Aw,w N

 Aw,w 1,1  w,w  A2,1   . :=   ..   Aw,w  N −1,1

Aw,w 1,2

...

Aw,w 1,N −1

Aw,w 2,2

...

Aw,w 2,N −1

.. .

..

.. .

Aw,w N −1,2

...

Aw,w N −1,N −1

Aw,w N ,2

...

Aw,w N ,N −1

Aw,w N ,1

.

 Aw,w 1,N  Aw,w 2,N    .. , .   w,w A N −1,N  

(23)

Aw,w N ,N

v,w and Aw,v N = AN . With all these preparations, we can decompose the matrix–vector multiplication A2N u2N as follows

A2N u2N = A2N uo2N + A2N ue2N = Ao2N uo2N + Ae2N ue2N .

(24)

In addition, if we collect all the odd numbered rows and all the even numbered rows of Ao2N uo2N , we obtain the following matrix–vector multiplications Av,v N vN ,

Aw,v N vN ,

(25)

respectively. Similarly, if we collect all the odd numbered rows and all the even numbered rows of Ae2N ue2N , we obtain the following matrix–vector multiplications Av,w N wN ,

Aw,w N wN ,

(26)

respectively. In summary, to evaluate the matrix–vector multiplication A2N u2N efficiently, we need only to efficiently evaluate w,v v,w w,w the matrix–vector multiplications Av,v N v N , A N v N , A N w N , and A N w N . Similarly, to efficiently store and assemble w,v v,w w,w the stiffness matrix A2N , we need only to efficiently store and assemble the submatrices Av,v N , A N , A N , and A N . We will address their efficient evaluation and storage in Section 5. v,v

v,w

w,w

4. Structure of matrices A2N , A N , A N , and A N

v,w w,w In this section, we analyze the structure of the matrices A2N , Av,v N , A N and A N . We begin with the analysis of the structure of the stiffness matrix A2N .

Theorem 1. Let K and L be defined by (3). The stiffness matrix A2N can be expressed as the following block-banded matrix of order (N y − 1) with a block bandwidth 2L + 1 B1,1 · · · B1,L+1 0  . . . . .. .. ..  ..   ..  L+1,1 . . . B L+1,L+1 . B   .. ..  0 . . B L+2,L+2   . .. .. .. = . . .  ..   .. ..  0 . . 0   .. ..  0 . . 0    .. .. .. .. . . .  . 

A2N

0

···

0

0 ′

··· .. . .. . .. . .. . .. . .. . .. .

0 .. .

0 .. . .. .

0 .. . 0 .. .. . . .. .. . . .. . B N y −L−1,N y −L−1 .. .. . .

··· 0

 ··· 0  .. ..  . .   ..  .  0   ..  . 0   .. .. . . .    ..  . 0   ..  N −L−1,N −1 y y . B   ..  .. .  .

B N y −1,N y −L−1 · · ·

(27)

B N y −1,N y −1 ′

Furthermore, each matrix block B j, j is of order 2(N x − 1). For | j − j ′ | 6 L, all the matrix blocks B j, j can be expressed as

287

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303



B j, j

′

j, j ′

B1,1  ..  .   j, j ′ B  K +1,1    0  .. =  .    0    0   ..  . 0

j, j ′

0 · · · B1,K +1 .. .. .. . . . .. .. j, j ′ . B K +1,K +1 . ′ .. .. j, j . . B K +2,K +2 .. .. .. . . . .. .. . . 0 .. .. . . 0 .. .. .. . . . ··· 0 0

··· 0 0 .. .. .. . . . .. .. . 0 . .. .. . . 0 .. .. .. . . . .. .. .. . . . .. .. j, j ′ . . B N −K −1,N −K −1 x x .. .. .. . . . j, j ′ · · · 0 B Nx −1,Nx −K −1

··· .. . .. . ..

.

..

.

.. .. ..

. .

. ···

0 .. .



     0     0  . ..  .    0   j, j ′ B Nx −K −1,Nx −1    ..  .

(28)

j, j ′

B Nx −1,Nx −1

Finally, for | j ′ − j| 6 L and |i ′ − i| 6 K , each 2-by-2 matrix block satisfies j, j ′

Bi,i ′ = Am,n , where the indices i,

(29) i ′,

j, and

j′

and the indices m and n are related by (12).

Proof. We first rewrite the stiffness matrix as the following form   1,1 B B1,2 ··· B1,N y −1   2,1 ..  B . B2,N y −1  B2,2 .  A2N =   .. .. .. ..   . . . . N y −1,1 N y −1,2 N y −1,N y −1 B B ··· B

(30)

′

Here each matrix block B j, j of order 2(N x − 1) represents the interaction of row j and row j ′ in the discrete system ′ for 1 6 j 6 N y − 1 and 1 6 j ′ 6 N y − 1. Furthermore, each matrix block B j, j can be expressed as  j, j ′  j, j ′ j, j ′ B1,1 B1,2 ··· B1,Nx −1   .. j, j ′ j, j ′  B j, j ′ . B2,2 B2,Nx −1    2,1 j, j ′ (31) B = . .. .. .. ..   . .   . . j, j ′ j, j ′ j, j ′ B Nx −1,1 B Nx −1,2 · · · B Nx −1,Nx −1 j, j ′

Here each matrix block Bi,i ′ of order 2 represents the interaction between the ith node in the j row of the grid and the i ′ th node in the j ′ row of the grid for 1 6 i 6 N x − 1 and 1 6 i ′ 6 N x − 1. In addition, we have j, j ′

Bi,i ′ = Am,n ,

(32)

where the indices i, i ′ , j, and j ′ and the indices m and n are related by (12). Thus to prove Theorem 1, it remains ′ j, j ′ to prove the matrix blocks B j, j with | j − j ′ | > L + 1 are zeros and the matrix blocks Bi,i ′ with | j − j ′ | 6 L and |i − i ′ | > K + 1 are zeros. We conclude from the expressions of the entries of the matrix A2N in (11) and the relations between indices in (12) that the matrix blocks Am,n = 0 if supp(φi ′ , j ′ ) ∩ Bδ (xi , y j ) = ∅.

(33) | j′

We deduce from the definition L in (3) that (33) holds for − j| > L + 1. Thus, the matrix (30) can be expressed in the form of (27). We next investigate the structure of each block matrix Bi, j with | j ′ − j| 6 L. We again use the definition K in (3) to conclude that (33) holds for |i ′ − i| > K + 1 and thus the matrix (31) can be expressed in the form of (28). Eq. (29) can be directly obtained from (32) for |i − i ′ | 6 K and | j − j ′ | 6 L. 

288

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

v,w w,w Theorem 2. Each of the matrices Av,v can be expressed as an (N y − 1)-by-(N y − 1) block-banded N , A N and A N Toeplitz matrix with a block bandwidth 2L + 1

T0I  ..  .   I T−L    0   I =  ... AN    0    0   .  .. 0 

··· .. . .. . .. . .. . .. . .. . .. . ···

T LI .. .

0 .. .

T0I .. .

.

..

.

0 .. . ..

. 0

..

T0I .. . ..

.

0 .. . 0

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. .

0 .. .

0 .. .

.

..

.

T0I .. . ..

. 0

..

0 .. . .. . T0I .. . I T−L

··· .. . .. . .. . .. . .. . .. . .. . ···

 0 ..  .    0   0  ..  . .    0   I TL  ..  .  T0I

(34)

where the superscript I = (v, v), (v, w) or (w, w). Furthermore, each matrix block T Ij , with −L 6 j 6 L, is an (N x − 1)-by-(N x − 1) banded Toeplitz matrix with a bandwidth 2K + 1 t0,I j  ..  .   I t  −K , j   0   . T Ij =   ..    0    0   .  .. 0 

··· .. . .. . .. . .. . .. . .. . .. . ···

t KI , j .. . t0,I j .. . .. . 0 .. . .. . 0

0 .. . .. . t0,I j .. . .. . 0 .. . 0

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. . 0 .. . .. . t0,I j .. . .. . 0

0 .. . .. . 0 .. . .. . t0,I j .. . I t−K ,j

··· .. . .. . .. . .. . .. . .. . .. . ···

 0 ..  .    0    0   ..   . .   0    t KI , j   ..  . 

(35)

t0,I j

v,w w,w Proof. We need only to prove (34)–(35) for Av,v are similar. By Theorem 1, A2N N . The proofs for A N and A N is a (2L + 1) block-banded, (2K + 1) block-banded block matrix in which each matrix block is a matrix of order 2. Following the proof of Theorem 1 and the definition (21), the matrix Av,v N is a (2L + 1) block-banded, (2K + 1) banded matrix. More precisely, Av,v is of the form N

B1,1 · · · B1,L+1 0 v v  .. .. .. ..  . . . .  ..  L+1,1 . . L+1,L+1 . Bv .  Bv   .. ..  0 . . BvL+2,L+2  . .. .. .. =  .. . . .   .. ..  0 . . 0  .. ..  . .  0 0  . .. .. ..  . . . . . 0 ··· 0 0 

Av,v N

··· 0 0 .. .. .. . . . .. .. . 0 . .. .. . . 0 .. .. .. . . . .. .. .. . . . . . . . N y −L−1,N y −L−1 . . Bv .. .. .. . . . N y −1,N y −L−1 · · · 0 Bv

··· .. . .. . ..

.

..

.

.. .. ..

. .

. ···

0 .. .



     0    0  , ..  .    0  N y −L−1,N y −1   Bv  ..  . N y −1,N y −1 Bv

(36)

289

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303 j, j ′

where each matrix block Bv 

′ Bvj, j

j, j ′

b1,1

  ..  .    b j, j ′  K +1,1     0   =  ..  .    0     0   .  .  . 0

is a banded matrix with a bandwidth 2K + 1 j, j ′

···

b1,K +1

0

··· 0

0

···

0

..

..

..

.

..

.

..

.

.. .

..

.

.. .. . . .. . 0

..

.

..

.

. . j, j ′ b K +2,K +2 . . . .

0

..

.

..

.

..

.

0

.. .. .. . . . .. .. .. . . . . . . . j, j ′ . . b N −K −1,N −K −1 x x

..

.

..

.. .. . .

..

..

.

··· 0

j, j ′ b Nx −1,Nx −K −1

..

.

.

′

j . b Kj,+1,K +1

..

.

..

.

..

.

..

.

..

.

..

.

0 .. .

..

.

..

···

.

0

..

.

..

.

.

0

.

···



      0     0    ..   .    0    j, j ′ b Nx −K −1,Nx −1    ..   .

(37)

j, j ′

b Nx −1,Nx −1

and j, j ′

bi,i ′ = Av,v m,n ,

|i − i ′ | 6 K , | j − j ′ | 6 L ,

(38)

where the indices i, i ′ , j, and j ′ , and the indices m and n are related by (12). By the following translations ξ1 = x ′ − xi ,

ξ2 = y ′ − y j ,

the pyramid function φi ′ , j ′ (x ′ , y ′ ) defined in (4) can be expressed as  x ′ − x ′   y′ − y ′  j i ψ hx hy ξ − x ′  ξ − y ′  2 j −j 1 i −i =ψ ψ = φi ′ −i, j ′ − j (ξ1 , ξ2 ). hx hy

φi ′ , j ′ (x ′ , y ′ ) = ψ

(39)

Therefore, the first equation in (11) can be deduced to Av,v m,n =

 Bδ (0,0)

  ξ1 2 σ (∥(ξ1 , ξ2 )∥) δm,n − φi ′ −i, j ′ − j (ξ1 , ξ2 ) dξ1 dξ2 . 2 2 ξ1 + ξ2

(40) j1 , j1′

To prove the block-Toeplitz property of the matrix Av,v N , we consider the matrix blocks Bv are on the same diagonal. Namely, their superscripts j1 and j1′ as well as j2 and j2′ satisfy j1′ − j1 = j2′ − j2 = l

j2 , j2′

and Bv

in (36) that

(41)

for some −L 6 l 6 L. Let m 1 := ( j1 − 1)(N x − 1) + i, n 1 :=

( j1′

− 1)(N x − 1) + i , ′

1 6 i 6 N x − 1, 1 6 j1 6 N y − 1, 1 6 i ′ 6 N x − 1, 1 6 j1′ 6 N y − 1,

m 2 := ( j2 − 1)(N x − 1) + i,

1 6 i 6 N x − 1, 1 6 j2 6 N y − 1,

n 2 := ( j2′ − 1)(N x − 1) + i ′ ,

1 6 i ′ 6 N x − 1, 1 6 j2′ 6 N y − 1.

(42)

290

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

Then, we use (38) and (41) to prove that for 1 6 i, i ′ 6 N x − 1 j1 , j1′

bi,i ′

= Av,v m 1 ,n 1  = Bδ (0,0)

  ξ1 2 σ (∥(ξ1 , ξ2 )∥) δm 1 ,n 1 − φi ′ −i, j ′ − j1 (ξ1 , ξ2 ) dξ1 dξ2 2 1 ξ1 + ξ2

Bδ (0,0)

  ξ1 2 σ (∥(ξ1 , ξ2 )∥) δm 2 ,n 2 − φi ′ −i, j ′ − j2 (ξ1 , ξ2 ) dξ1 dξ2 2 2 2 ξ1 + ξ2

 =

2

j2 , j ′

2 = Av,v m 2 ,n 2 = bi,i ′ .

Therefore, we have proved that j1 , j1′

Bv

j2 , j2′

= Bv

, j1 , j ′

j2 , j ′

if the matrix blocks Bv 1 and Bv 2 are on the same diagonal in (36). Accordingly, we have proved that Av,v N is a block-banded Toeplitz matrix with a block bandwidth 2L + 1 and the matrix (36) can be expressed in the form of (34). j, j ′ j, j ′ j, j ′ Next, we fix each block matrix Bv with | j ′ − j| 6 L and look at their entries bi ,i ′ and bi ,i ′ that are on the same diagonal. Namely, their subscripts i 3 and i 3′ as well as i 4 and i 4′ satisfy

3 3

4 4

i 3′ − i 3 = i 4′ − i 4 = k for some −K 6 k 6 K . Let m 3 := ( j − 1)(N x − 1) + i 3 ,

1 6 i 3 6 N x − 1, 1 6 j 6 N y − 1,

− 1)(N x − 1) + i 3′ ,

1 6 i 3′ 6 N x − 1, 1 6 j ′ 6 N y − 1,

m 4 := ( j − 1)(N x − 1) + i 4 ,

1 6 i 4 6 N x − 1, 1 6 j 6 N y − 1,

− 1)(N x − 1) + i 4′ ,

1 6 i 4′ 6 N x − 1, 1 6 j ′ 6 N y − 1.

n 3 := ( j n 4 := ( j

′

′

We again apply (38) to obtain j, j ′ ′ 3 ,i 3

bi

= Av,v m 3 ,n 3

Bδ (0,0)

  ξ1 2 σ (∥(ξ1 , ξ2 )∥) δm 3 ,n 3 − φi ′ −i3 , j ′ − j (ξ1 , ξ2 ) dξ1 dξ2 2 2 3 ξ1 + ξ2

Bδ (0,0)

  ξ1 2 σ (∥(ξ1 , ξ2 )∥) δm 4 ,n 4 − φi ′ −i4 , j ′ − j (ξ1 , ξ2 ) dξ1 dξ2 2 2 4 ξ1 + ξ2

 =  =

j, j ′ ′. 4 ,i 4

= Av,v m 4 ,n 4 = bi

j, j ′

Hence we conclude that each matrix block Bv in (36) with | j − j ′ | 6 L is a (2K + 1) banded Toeplitz matrix and that the matrix (37) can be expressed in the form of (35).  5. Efficient storage and assembly of the stiffness matrix and a fast Krylov subspace iterative method v,w w,w We utilize the structure of the matrices Av,v proved in Section 4 to develop an efficient storage N , A N , and A N and assembly mechanism for the stiffness matrix A2N . v,w w,w can be accomplished by Theorem 3. The evaluation and storage of each of the matrices Av,v N , A N , and A N evaluating and storing only (2K +1)(2L +1) entries. More precisely, we need only to evaluate and store the following

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

291

matrices of order (2K + 1)(2L + 1)

G v,v

 v,v t−K ,−L  ..  .  v,v :=   t0,−L  .  .. t Kv,v,−L

 v,w t−K ,−L  ..  .  v,w :=   t0,−L  .  ..

G v,w

t Kv,w ,−L

··· .. . ··· .. . ···

v,v t−K ,0 .. . v,v t0,0 .. . t Kv,v,0

··· .. . ··· .. . ···

··· .. . ··· .. . ···

v,w t−K ,0 .. . v,w t0,0 .. . v,w t K ,0

··· .. . ··· .. . ···

··· .. . ··· .. . ···

w,w t−K ,0 .. . w,w t0,0 .. . t Kw,w ,0

··· .. . ··· .. . ···

v,v  t−K ,L ..  .  v,v  , t0,L  ..  .  v,v t K ,L v,w  t−K ,L ..  .  v,w  t0,L  , ..  .  t Kv,w ,L

and

G w,w

 w,w t−K ,−L  ..  .  w,w :=   t0,−L  .  .. t Kw,w ,−L

w,w  t−K ,L ..  .  w,w  . t0,L  ..  .  w,w t K ,L

v,w w,w Proof. We need only to prove the theorem for the matrix Av,v are similar. By the N . The proofs for A N and A N v,v expression (34) in Theorem 2, to evaluate and store the matrix A N we need only to evaluate and store the 2L + 1 matrices Tlv,v for −L 6 l 6 L. By the expression (35) in Theorem 2, to evaluate and store each matrix Tlv,v , we need v,v only to evaluate and store the entries tk,l for −K 6 k 6 K and −L 6 l 6 L. 

Remark 1. The indices (k, l) represent the relative positions of points (xi , y j ) and (xi ′ , y j ′ ) for |i − i ′ | 6 K and j, j ′

I = tI I | j − j ′ | 6 L, i.e., k = i ′ − i and l = j ′ − j. Therefore, we have tk,l i ′ −i, j ′ − j = bi,i ′ = Am,n with I = (v, v), (v, w), and (w, w), where the second equality is obtained from the fact that (35) is equal to (37), the third equality is obtained by (38) or (29) and the indices i, i ′ , j, and j ′ as well as the indices m and n are related by (12). v,w Remark 2. The mechanism in Theorem 3 reduces the evaluation and storage of the stiffness matrices Av,v N , AN w,w 2 2 and A N from O(N ) to O((2K + 1)(2L + 1)) = O(δ /(h x h y )) = O(N x N y ) = O(N ) for a fixed δ. Hence, the stiffness matrix A2N can be stored in O(N ) memory instead of O(N 2 ) memory and can be assembled by evaluating only O(N ) entries instead of O(N 2 ) entries. v,w w,w Theorem 4. For any vector d N ∈ R N , the matrix–vector multiplications Av,v N d N , A N d N , and A N d N can be evaluated in O(N log N ) operations. v,w v,v Proof. We need only to prove that Av,v N d N can be evaluated in O(N log N ) operations. A N d N and A N d N can be v,v v,v computed by the similar way. By Theorem 2, A N is BTTB. Then each Toeplitz matrix T j of order N x − 1 in (34) can be embedded into a 2(N x − 1)-by-2(N x − 1) circulant matrix [36,37]

Cv,v j

 :=

Tv,v j

T˜ v,v j

˜ v,v T j

Tv,v j

 ,

(43)

292

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

where the matrix T˜ v,v j is defined by

T˜ v,v j

 0  ..  .    0    0   :=  ...    0   t v,v  K, j  .  . . v,v t1, j

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. .

0 .. .

0

0

0 .. .

0 .. . .. .

0 .. . .. .

t Kv,v, j

0 .. . 0

v,v t−K ,j .. . .. .

0 .. .

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. . .. . .. . .. . .. . 0

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. . .. . 0 .. . 0

v,v  t−1, j ..  .   v,v  t−K , j    0   ..  . .    0    0   ..   . 0

(44)

v,v Replacing each Tv,v j in (34) by the corresponding circulant matrix C j , we obtain a block-Toeplitz–circulant-block v,v (BTCB) matrix D

Cv,v 0  ..  .   v,v C−L    0   :=  ...    0    0   .  .. 0 

D

v,v

··· .. . .. . .. . .. . .. . .. . .. . ···

Cv,v L .. . Cv,v 0 .. . .. . 0 .. .

..

. 0

0 .. .

..

.

Cv,v 0 .. . .. . 0 .. . 0

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. .

0 .. .

..

0 .. . ..

0 .. . .. .

.

Cv,v 0 .. . .. . 0

.

Cv,v 0 .. . v,v C−L

··· .. . .. . .. . .. . .. . .. . .. . ···

 0 ..  .    0    0   ..  . .    0    v,v  CL  ..  .  v,v C0

(45)

Then Dv,v can be embedded into a 2(N y − 1)-by-2(N y − 1) circulant block matrix Cv,v as follows [36,37] Cv,v :=



Dv,v ˜ v,v D

 ˜ v,v D , Dv,v

(46)

where 0  ..  .    0    0   :=  ...    0   Cv,v  L  .  .. Cv,v 1 

˜ v,v D

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. .

0 .. .

0

0

0 .. .

0 .. . .. .

0 .. . .. . Cv,v L

0 .. . 0

··· .. . .. . .. . .. . .. . .. . .. . ···

0 .. . 0 .. . .. . .. . .. . .. . 0

Cv,v −L .. . .. . 0 .. . .. . 0 .. . 0

··· .. . .. . .. . .. . .. . .. . .. . ···

 Cv,v −1 ..  .     Cv,v −L   0   ..  . .    0    0   ..  .  0

(47)

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

293

Note that Cv,v is a 2(N y − 1)-by-2(N y − 1) block-circulant–circulant-block (BCCB) matrix with 2(N x − 1)-by2(N x − 1) circulant blocks. Let cv,v be its first column vector and cˆ v,v be the Fourier transform of cv,v   cˆ v,v := F2(N y −1) ⊗ F2(Nx −1) cv,v . (48) Here F2(N y −1) ⊗ F2(Nx −1) represents the two-dimensional Fourier matrix of order 4N . It is known that Cv,v has the following diagonalization [36,37],  −1    Cv,v := F2(N y −1) ⊗ F2(Nx −1) diag cˆ v,v F2(N y −1) ⊗ F2(Nx −1) . (49) Then we take the following procedures to carry out the fast matrix–vector multiplication Av,v N d N for any vector dN ∈ RN : 1. Let the vector d N be denoted by d N = [d(1) , d(2) , . . . , d(N y −1) ]T = [d1,1 , . . . , d Nx −1,1 , d1,2 , . . . , d Nx −1,2 , . . . , d1,N y −1 , . . . , d Nx −1,N y −1 ]T . Then expand this vector into a 2N -dimensional vector d2N by d2N = [d(1) , 0, d(2) , 0, . . . , d(N y −1) , 0]T , where 0 represents a N x − 1-dimensional zero vector. 2. Define d4N by d4N = [d2N , 0]T , where 0 represents a 2N -dimensional zero vector.  3. Compute the two-dimensional fast Fourier transform cˆ v,v by (48) and w1 = F2(N y −1) ⊗ F2(Nx −1) d4N . Evaluate the Hadamard product w2 = cˆ v,v · w1 and then carry out the two-dimensional inverse fast Fourier transform  −1 w3 = F2(N y −1) ⊗ F2(Nx −1) w2 . 4. Keep only the first 2N entries of the vector w3 to obtain a vector w4 and remove every other I − 1 entries of the vector w4 to obtain the vector Av,v N dN . We can conclude from the above procedures that the main computational work required for Av,v N d N is the computations related to the fast Fourier transform in step 3. It needs O(4N log(4N )) = O(N log N ) operations via the fast Fourier transform.  The efficient evaluation and storage of the stiffness matrix and the fast matrix–vector multiplications proved in this section can be applied in the development of any fast Krylov subspace iterative methods. Here, we take the standard conjugate gradient squared (CGS) method for the matrix equation (7) as an example to demonstrate the idea [38,39]. Note that in the conventional conjugate gradient squared (CGS) method the matrix–vector multiplication requires O(N 2 ) operations, while other operations can be performed in O(N ) operations. Furthermore, the conventional evaluation and storage of the stiffness matrix A2N require O(N 2 ) memory and computations. The efficient mechanism developed in this section reduces the storage and assembly of the stiffness matrix A2N to O(N ). w,w v,w We also proved that for any d N ∈ R N the matrix–vector multiplications Av,v N d N , A N d N and A N d N can be evaluated in O(N log N ) operations. Finally, we address the fast matrix–vector multiplication A2N u2N for any u2N ∈ R2N : (i) We first express u2N in the form of (8). (ii) We collect odd and even entries of u2N , respectively, to obtain two N -dimensional vectors v N and w N as defined in (15). (iii) We apply Theorem 4 to evaluate Av,v N vN , v,w w,v w,w w,v w,w v,v w and A v +A v , and A w in O(N log N ) operations. (iv) We embed A v +A Av,w w , A N N N N N N N N N wN N N N N into the odd and even entries of a 2N -dimensional vector, respectively, to obtain A2N u2N . 6. Numerical experiments Note that the fast collocation method is merely a traditional collocation method, which significantly reduces the number of entries of the stiffness matrix to be evaluated and assembled as well as the computational complexity of the matrix–vector multiplication in any Krylov subspace iterative methods, and so naturally inherits the stability, convergence, and other numerical advantages of traditional collocation methods. In this section we carry out numerical experiments for relatively simple model problems to investigate the improvements of the fast collocation method which is solved by the fast CGS (FCGS) method over the traditional collocation method which is solved with widely used Gaussian elimination (Gauss) and the conventional CGS method. All the solvers are implemented using Matlab and all the numerical experiments are performed on a computer with 8-GB memory.

294

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

In the numerical experiments, we consider the peridynamic model (1) with the kernel function [1,35] σ (|(x, y)|) = 

1 x 2 + y2

1+s ,

s < 1/2.

(50)

The domain Ω = (0, 1)×(0, 1) and the neighborhood Bδ (x) is assumed to be a closed disk of radius δ = 1/8 centered at x, where δ is referred to as the peridynamic horizon. For simplicity, we choose h x = h y = h. In the numerical experiments we use Gauss numerical quadrature to compute the matrix entries in (11). Furthermore, we use a linear regression to fit the convergence rate α and the associated constant Cα in the error estimate ∥u − uh ∥ L 2 (Ω ) 6 Cα h α .

(51)

6.1. A peridynamic model with a continuous displacement The exact solution u(x, y) = (v(x, y), w(x, y)) is chosen to be v(x, y) = w(x, y) = x(1 − x)y(1 − y),

(x, y) ∈ Ω

which are also used to define g v and g w on the boundary zone Ωc . The right-hand side external forcing term b(x, y) can be computed analytically by polar coordinates transformation as follows: 3π δ 2−2s π δ 2−2s π δ 4−2s (y − y 2 ) + (3x + 2y − 4x y − x 2 − 1) − 8 − 8s 8 − 8s 32 − 16s 2−2s 2−2s 3π δ πδ π δ 4−2s bw (x, y) = (x − x 2 ) + (2x + 3y − 4x y − y 2 − 1) − . 8 − 8s 8 − 8s 32 − 16s bv (x, y) =

(52)

Example 1. We consider the kernel function (50) with s = 3/8. We present the L 2 errors of the Gauss, CGS, and FCGS numerical solutions, the CPU time consumed by these methods, and the number of iterations by CGS and FCGS for the mesh sizes from 1/24 to 1/29 in Table 1. We observe that all the three solvers generate identical numerical solutions as we have anticipated, since the fast collocation method is merely the traditional one with a much efficient evaluation and assembly of the stiffness matrix and a fast evaluation of the matrix–vector multiplication. The numerical results in the table justify that the fast method has a significant improvement of computational complexity over the traditional one. The direct solver uses the most CPU time and runs out of memory for h = 1/27 , while both CGS and FCGS still work for the mesh size h = 1/27 , 1/28 , and 1/29 . In fact, h = 1/27 is the finest mesh size for the CGS to work on this particular computer. However, for h = 1/28 and h = 1/29 , the FCGS still works because of less memory requirement. In principle, the FCGS still applies to the problem with much finer mesh sizes. The numerical results also show that the FCGS has a significantly reduced CPU time as anticipated. We can also conclude from Table 1 that the convergence rate α of these methods is sublinear. For clarity, we present a representative figure of the collocation approximations for the displacement field components v(x, y) and w(x, y) on a 64 × 64 mesh in Fig. 1. Example 2. We consider the case s = 0 in (50), which corresponds to a borderline case between a nonintegrable kernel and an integrable kernel. We present the numerical results in Table 2. We virtually have the same observations as in Example 1, except that we observe better accuracy and convergence rates of the numerical solutions and less number of iterations as a reflection of improved condition number of the stiffness matrix. The corresponding numerical solutions for the displacement field components v(x, y) and w(x, y) are shown in Fig. 2. Example 3. We choose s = −1/2 in (50). In this case the kernel σ ∈ L 1 (Ω ) but σ ̸∈ L 2 (Ω ). We present numerical results in Table 3. We again have the same observations as in Examples 1 and 2, except that a further improved accuracy and a convergence rate of almost 2 of the numerical solutions and further reduced number of iterations. Furthermore, the corresponding numerical solutions for the displacement field components v(x, y) and w(x, y) on a 64 × 64 mesh are shown in Fig. 3.

295

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303 Table 1 Performance of Gauss, CGS, and the FCGS in Example 1.

Gauss

Nx = N y

∥uh − u∥ L 2

24

3.2905 × 10−2

25 26 27 28 29

2.3835 × 10−2 1.4337 × 10−2

# of iter.

CPUs 1.23 s 1 min 30 s 16 h 25 min

Out of memory Out of memory Out of memory Cα = 0.18, α = 0.60

CGS

24 25 26 27 28 29

3.2905 × 10−2 2.3835 × 10−2 1.4337 × 10−2 7.4180 × 10−3

69 138 216 338 Out of memory Out of memory

0.12 s 2.07 s 50 s 1 d 12 h

73 138 225 340 478 685

0.44 s 2.3 s 31.8 s 3 min 46 s 21 min 18 s 3 h 42 min

Cα = 0.26, α = 0.72

FCGS

24 25 26 27 28 29

3.2905 × 10−2 2.3835 × 10−2 1.4337 × 10−2 7.4180 × 10−3 5.7085 × 10−3 2.5619 × 10−3 Cα = 0.28, α = 0.73

Fig. 1. The collocation approximations for v(x, y)(left) and w(x, y)(right) on a 64 × 64 mesh in Example 1.

6.2. A peridynamic model with a discontinuous displacement We simulate the peridynamic model (1) with a discontinuous displacement  x + y, x < 0.5 v(x, y) = x 2 + y 2 , x > 0.5, w(x, y) = x(1 − x)y(1 − y). The computational nodes are placed along the discontinuity 0.5. We consider the kernel function (50) with s = −1/2 as in Example 3. The right-hand side term b is computed by the numerical integration. We present the numerical results in Table 4. We have similar observations as in Example 3, except a greatly reduced convergence rate O(h 1/2 ).

296

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303 Table 2 Performance of Gauss, CGS, and FCGS in Example 2.

Gauss

Nx = N y

∥uh − u∥ L 2

24

1.6387 × 10−2

25 26 27 28 29

6.8624 × 10−3 2.3722 × 10−3

# of iter.

CPUs 1.22 s 1 min 30 s 16 h 45 min

Out of memory Out of memory Out of memory Cα = 0.81, α = 1.39

CGS

24 25 26 27 28 29

1.6387 × 10−2 6.8624 × 10−3 2.3722 × 10−3 9.4000 × 10−4

48 67 82 121 Out of memory Out of memory

0.07 s 1.26 s 22.3 s 14 h 26 m

48 71 87 121 133 151

0.42 s 1.31 s 12 s 1 min 55 s 12 min 3 s 1 h 15 min

Cα = 0.80, α = 1.39

FCGS

24 25 26 27 28 29

1.6387 × 10−2 6.8624 × 10−3 2.3722 × 10−3 9.4000 × 10−4 4.3047 × 10−4 2.2423 × 10−4 Cα = 0.48, α = 1.25

Fig. 2. The collocation approximations for v(x, y)(left) and w(x, y)(right) on a 64 × 64 mesh in Example 2.

Because the true solutions both lack enough regularity, this convergence rate is the same as the Galerkin finite element method had in [27]. For clarity, we show the numerical solutions for the displacement field components v(x, y) and w(x, y) on a 64 × 64 mesh in Fig. 4. From this figure, we can see the numerical solutions at the points where there is a discontinuity have a poor performance. 7. Conclusions In this paper we develop a fast collocation method for a static bond-based linear peridynamic model (1) in two space dimensions. The method reduces the computational work from O(N 2 ) per Krylov subspace iteration to O(N log N ) and the memory requirement from O(N 2 ) to O(N ). In addition, the method reduces the computational work of

297

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303 Table 3 Performance of Gauss, CGS, and FCGS in Example 3.

Gauss

Nx = N y

∥uh − u∥ L 2

24

9.6360 × 10−3

25 26 27 28 29

3.0229 × 10−3 8.2596 × 10−4

# of iter.

CPUs 1.23 s 1 min 31 s 16 h 30 min

Out of memory Out of memory Out of memory Cα = 1.34, α = 1.77

CGS

24 25 26 27 28 29

9.6360 × 10−3 3.0229 × 10−3 8.2596 × 10−4 1.9064 × 10−4

41 48 51 63 Out of memory Out of memory

0.12 s 0.94 s 15.7 s 6 h 18 min

41 48 51 63 71 80

0.33 s 0.95 s 7.6 s 43.6 s 5 min 1 s 35 min 13 s

Cα = 1.93, α = 1.89

FCGS

24 25 26 27 28 29

9.6360 × 10−3 3.0229 × 10−3 8.2596 × 10−4 1.9064 × 10−4 7.2990 × 10−5 2.2827 × 10−5 Cα = 1.27, α = 1.77

Fig. 3. The collocation approximations for v(x, y)(left) and w(x, y)(right) on a 64 × 64 mesh in Example 3.

evaluating and assembling the stiffness matrix, which often constitutes a major portion of the overall computational work, from O(N 2 ) to O(N ). In fact, the fast collocation method is nothing but the conventional one, and so naturally inherits the stability, convergence, and other numerical advantages of conventional collocation methods [40,41,23, 42,15,25,43]. Furthermore, any numerical techniques developed for traditional collocation methods, such as the numerical quadratures for the accurate evaluation of singular integrals [12,13,31–33,22,15,34], can be applied to the fast collocation method. In the fast collocation method we take the full advantage of the structure of its stiffness matrix in its storage, evaluation, and assembly and in the efficient solution of the corresponding matrix equation. This would significantly reduce the computational complexity and storage of the fast collocation method over conventional ones, without

298

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303 Table 4 Performance of Gauss, CGS, and FCGS in the example in Section 6.2.

Gauss

Nx = N y

∥uh − u∥ L 2

24

7.2376 × 10−2

25 26 27 28 29

4.1886 × 10−2 2.7860 × 10−2

# of iter.

CPUs 1.25 s 1 min 27 s 17 h 5 min

Out of memory Out of memory Out of memory Cα = 0.48, α = 0.69

CGS

24 25 26 27 28 29

7.2376 × 10−2 4.1886 × 10−2 2.7860 × 10−2 1.9248 × 10−2

61 71 76 83 Out of memory Out of memory

0.14 s 1.25 s 20.2 s 8 h 51 min

61 71 76 83 85 88

0.53 s 1.22 s 10 s 56.4 s 5 min 15 s 36 min 4 s

Cα = 0.40, α = 0.63

FCGS

24 25 26 27 28 29

7.2376 × 10−2 4.1886 × 10−2 2.7860 × 10−2 1.9248 × 10−2 1.3491 × 10−2 9.5472 × 10−3 Cα = 0.32, α = 0.57

Fig. 4. The collocation approximations for v(x, y)(left) and w(x, y)(right) on a 64 × 64 mesh in the example in Section 6.2.

using any lossy compression. For instance, we would use the same numerical quadratures for conventional methods to evaluate the singular integrals in the stiffness matrix, except that we only need to evaluate O(N ) of them instead of O(N 2 ) of them. Numerical results are presented to show the utility of the fast method. Technically, the contributions of this paper are summarized as follows: (i) The development of the fast collocation method for the peridynamic model (1) is motivated by that of the fast finite difference methods for FPDEs [44,45, 39,46–48]. However, the peridynamic model (1) gives rise to the coupling in all the directions, in contrast to FPDEs in which the fractional derivatives are expressed only in the coordinate directions, and so represents a much stronger coupling than FPDEs do. In this paper we analyze and utilize the structure of the stiffness matrix to develop a fast collocation method with efficient matrix assembly and storage. (ii) The fast numerical methods developed for FPDEs

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

299

[39,46] and nonlocal diffusion models [29,26] are primarily for scalar problems. However, the peridynamic model (1) is a vector equation. Due to the coupling between the x component v and the y component w of the displacement vector u introduced by the tensor product kernel in (2), the stiffness matrix in the collocation method is not BTTB as that for nonlocal diffusion models [26]. Thus, the fast numerical methods developed for FPDEs and nonlocal diffusion models do not apply to the peridynamic model (1). In this paper we analyze the coupling effect of the peridynamic model (1) on the stiffness matrix of the collocation method and utilize the structure in the development of a fast collocation method. We end this section by discussing possible extensions and limitations of the current method: (i) The fast collocation method developed in this paper can be extended to time-dependent problems in a straightforward manner. In fact, the condition number of the stiffness matrix is improved greatly due to the impact of the time step size. (ii) The extension to three-dimensional problem is conceptually straightforward as the stiffness matrix should retain a similar structure, but with one more layer of matrix blocks. However, the practical implementation requires tremendous effort in terms of evaluating the three-fold integrals with singular kernels. This will be investigated in the near future. (iii) The extension to state-based peridynamic models [2] is much more involved, as the corresponding peridynamic model doubles the folds of integrals. For example, one needs to calculate a four-fold integral for each entry of the stiffness matrix for a two-dimensional problem. It is not clear whether the stiffness matrix still possesses a somewhat similar structure in this context. This issue will be investigated in the future. (iv) The peridynamic models for problems with material heterogeneities and anisotropy contain a variable-coefficient inside the integral operator [49]. It is not clear whether the idea in the current paper carries over. This issue will be investigated in the future. Acknowledgments The authors would like to express their sincere thanks to Dr. Pablo Seleson for the very helpful discussions about numerical approximations of peridynamic models. The authors would also like to express their sincere thanks to the referees for their very helpful comments and suggestions, which greatly improved the quality and readability of this paper. This work was partially supported by the OSD/ARO MURI Grant W911NF-15-1-0562, by the National Science Foundation under Grants DMS-1216923 and DMS-1620194, and the National Natural Science Foundation of China under Grants 91130010, 11471194, and 11571115, and the Taishan research project of Shandong Province (10000082961001). Appendix In this section, we give a detailed description of the fast method by taking an example. We consider a small system by choosing N x = N y = 5, L = K = 1, and h x = h y = 1 in the discrete system discussed in Section 2. The kernel function σ (x, y) is also chosen as (50). In this small system, we have 16 collocation points at which the values of v(x, y) and w(x, y) need to be computed, which are illustrated in Fig. 5. Then we only consider the evaluation and construction of the matrix Tv,v ∈ R16,16 for j = 0, −1, and 1. The j v,w w,w evaluations and constructions of the matrices T j and T j are carried out similarly. By the proof of Theorem 2, the value of each entry in Tv,v j is only determined by the relative position of (x i , y j ) and (x i ′ , y j ′ ), where 1 ≤ i, j ≤ 4 ′ ′ and 1 ≤ i , j ≤ 4. Thus we can choose arbitrary collocation point to compute all the entries in the matrix Tv,v j for j = 0, −1, and 1. Without loss of generality, we choose the point (x2 , y3 ) as the reference point as depicted in Fig. 5. To construct v,v the matrix Tv,v j , we need to compute the matrix entries Am,n defined in (11), where m and n are related to the indices ′ ′ (i, j) and (i , j ) by m = 4( j − 1) + i, i = 2, j = 3, n = 4( j ′ − 1) + i ′ , 1 6 i ′ 6 3, 2 6 j ′ 6 4.

(53)

In fact, these entries can be computed analytically. For example, by the transformations ξ1 = x ′ − x2 and ξ2 = y ′ − y3 , we have  (x ′ − x2 )2 ′ ′ ′ ′ Av,v = −  2+s (1 + x − x3 )(1 + y − y4 )d x dy 10,15 Ω1 (x ′ − x 2 )2 + (y ′ − y3 )2

300

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

Fig. 5. A reference grid point.



1



1−ξ22

=− 0



(ξ12 + ξ22 )2+s

0 1  π/2

=− 0

ξ13 ξ2

dξ1 dξ2

r 1−2s cos3 θ sin θdθ dr

0

1 =− , 8 − 8s

(54)

  where Ω1 = B1 (x2 , y3 ) ∩ (x2 , x3 ) × (y3 , y4 ) and we have used polar coordinates transformation. Other entries can be computed similarly. v,v Let α := i ′ − 2 and β := j ′ − 3 represents the relative position of the points (xi , y j ) and (xi ′ , y j ′ ). Let tα,β = Av,v m,n , ′ ′ where m and n are related to the indices (i, j) and (i , j ) by (53). Then we can store the 9 nonzero entries by the following matrix of order 3

Gv,v

 v,v t−1,−1  t v,v =  0,−1 v,v t1,−1

v,v t−1,0 v,v t0,0

v,v t1,0

v,v t−1,1



v,v  t0,1 . v,v t1,1

v,v as follows: Then we generate the matrix Tv,v j for j = −1, 0, and 1 using the matrix G



v,v t0,−1

v,v t1,−1

0

0



 v,v  v,v v,v t1,−1 0  t−1,−1 t0,−1   Tv,v = v,v v,v v,v  , −1  0 t−1,−1 t0,−1 t1,−1   v,v v,v 0 0 t−1,−1 t0,−1  v,v  v,v t0,0 t1,0 0 0  v,v  v,v v,v t1,0 0  t−1,0 t0,0 v,v   T0 =  v,v v,v v,v  , 0 t t t  −1,0 0,0 1,0  v,v v,v 0 0 t−1,0 t0,0

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

 Tv,v 1

v,v t0,1

 v,v t−1,1 =  0  0

v,v t1,1

0

v,v t0,1

v,v t1,1

0

v,v t−1,1

v,v t−1,1

v,v t0,1

0

301



 0   v,v  . t1,1  v,v t0,1

Then we construct the BTTB matrix Av,v 16 in the following way  v,v  T0 Tv,v 0 0 1  v,v  v,v v,v T1 0  T−1 T0 v,v   A16 =  v,v v,v v,v  . 0 T T T  −1 0 1  v,v 0 0 Tv,v T −1 0 16 Next we give a detailed description of the matrix–vector multiplication Av,v N d for any d ∈ R . It is known that the v,v v,v matrix T j can be embedded into a 8 × 8 circulant matrix C j for j = −1, 0, or 1 as follows:  v,v  Tj T˜ v,v j v,v Cj = , ˜ v,v Tv,v T j j

where 

T˜ v,v j

0  0 =  0 v,v t1, j

0 0 0 0 0 0 0 0

v,v  t−1, j 0  . 0  0

v,v v,v v,v Then we obtain the BTCB matrix Dv,v by replacing the matrix Tv,v −1 , T0 , and T1 in the matrix A16 with the v,v v,v v,v circulant matrix C−1 , C0 , and C1 , respectively, i.e.,  v,v  C0 Cv,v 0 0 1  v,v  Cv,v 0  C−1 Cv,v 0 1 v,v  . D = v,v  Cv,v Cv,v  0 −1 C0 1  v,v 0 0 Cv,v −1 C0

The matrix Dv,v can also be embedded into a BCCB matrix Cv,v of order 32 as follows:   v,v ˜ v,v D D Cv,v = ˜ v,v , D Dv,v where 

˜ v,v D

0  0 =  0 Cv,v 1

 0 0 Cv,v −1 0 0 0  . 0 0 0  0 0 0

Let cv,v be the first column vector of the matrix Cv,v . Let F8 ⊗ F8 be the two-dimensional discrete Fourier transform matrix. The Fourier transform of cv,v is given by c˜ v,v = (F8 ⊗ F8 )cv,v . Then the matrix Cv,v has the following diagonalization Cv,v = (F8 ⊗ F8 )−1 diag(˜cv,v )(F8 ⊗ F8 ).

302

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

Then we can carry out the fast matrix–vector Av,v N d by the procedures in the proof of Theorem 3. Not that we use the Matlab functions fft2 and ifft2 to compute the two-dimensional fast Fourier transform and the two-dimensional inverse fast Fourier transform, respectively. References [1] S. Silling, Reformulation of elasticity theory for discontinuous and long-range forces, J. Mech. Phys. Solids 48 (2000) 175–209. [2] S. Silling, M. Epton, O. Wecker, J. Xu, E. Askari, Peridynamic states and constitutive modeling, J. Elasticity 88 (2007) 151–184. [3] W. Hu, Y.D. Ha, F. Bobaru, Peridynamic model for dynamic fracture in unidirectional fiberreinforced composites, Comput. Methods Appl. Mech. Engrg. 217–220 (2012) 247–261. [4] B. Kilic, A. Agwai, E. Madenci, Peridynamic theory for progressive damage prediction in centercracked composite laminates, Compos. Struct. 90 (2009) 141–151. [5] E. Oterkus, E. Madenci, O. Weckner, S. Silling, P. Bogert, Combined finite element and peridynamic analyses for predicting failure in a stiffened composite curved panel with a central slot, Compos. Struct. 94 (2012) 839–850. [6] J. Xu, A. Askari, O. Weckner, S. Silling, Peridynamic analysis of impact damage in composite laminates, J. Aerosp. Eng., SPECIAL ISSUE: Impact Mechanics of Composite Materials for Aerospace Application 21 (2008) 187–194. [7] Y.D. Ha, F. Bobaru, Studies of dynamic crack propagation and crack branching with peridynamics, Int J. Fract. 162 (2010) 229–244. [8] Y.D. Ha, F. Bobaru, Characteristics of dynamic brittle fracture captured with peridynamics, Eng. Fract. Mech. 78 (2011) 1156–1168. [9] B. Kilic, E. Madenci, Prediction of crack paths in a quenched glass plate by using peridynamic theory, Int. J. Fract. 156 (2009) 165–177. [10] S. Silling, O. Weckner, E. Askari, F. Bobaru, Crack nucleation in a peridynamic solid, Int. J. Fract. 162 (2010) 219–227. [11] K. Dayal, K. Bhattacharya, Kinetics of phase transformations in the peridynamic formulation of continuum mechanics, J. Mech. Phys. Solids 54 (2006) 1811–1842. [12] F. Bobaru, Y.D. Ha, Adaptive refinement and multiscale modeling in 2D peridynamics, Int’l J. Multiscale Comput. Engrg. 9 (2011) 635–659. [13] F. Bobaru, Y.D. Ha, W. Hu, Damage progression from impact in layered glass modeled with peridynamics, Central European J. Engrg. 2 (2012) 551–561. [14] P. Seleson, M.L. Parks, On the role of the infuence function in the peridynamic theory, Int’l J. Multiscale Comput. Engrg. 9 (2011) 689–706. [15] S. Silling, E. Askari, A meshfree method based on the peridynamic model of solid mechanics, Comput. Struct. 83 (2005) 1526–1535. [16] M.R. Tupek, J.J. Rimoli, R. Radovitzky, An approach for incorporating classical continuum damage models in state-based peridynamics, Comput. Methods Appl. Mech. Engrg. 263 (2013) 20–26. [17] M. Ghajari, L. Iannucci, P. Curtis, A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media, Comput. Methods Appl. Mech. Engrg. 276 (2014) 431–452. [18] S. Sun, V. Sundararaghavan, A peridynamic implementation of crystal plasticity, Int’l J. Solids and Struct. 51 (2014) 3350–3360. [19] W. Gerstle, N. Sau, S. Silling, Peridynamic modeling of concrete structures, Nucl. Eng. Des. 237 (2007) 1250–1258. [20] X. Lai, B. Ren, H. Fan, S. Li, C.T. Wu, R.A. Regueiro, L. Liu, Peridynamics simulations of geomaterial fragmentation by impulse loads, Int. J. Numer. Anal. Methods Geomech. 39 (2015) 1304–1330. [21] E. Emmrich, O. Weckner, The peridynamic equation and its spatial discretisation, Math. Model. Anal. 12 (2007) 17–27. [22] P. Seleson, Improved one-point quadrature algorithms for two-dimensional peridynamic models based on analytical calculations, Comput. Methods Appl. Mech. Engrg. 282 (2014) 184–217. [23] P. Seleson, Q. Du, M.L. Parks, On the consistency between nearest-neighbor peridynamic discretizations and discretized classical elasticity models, Comput. Methods Appl. Mech. Engrg. (2016). http://dx.doi.org/10.1016/j.cma.2016.07.039. [24] P. Seleson, D. Littlewood, Convergence studies in meshfree peridynamic simulations, Comput. Math. Appl. 71 (2016) 2432–2448. [25] X. Tian, Q. Du, Analysis and comparison of different approximations to nonlocal diffusion and linear peridynamic equations, SIAM J. Numer. Anal. 51 (2013) 3458–3482. [26] H. Wang, H. Tian, A fast and faithful collocation method with efficient matrix assembly for a two-dimensional nonlocal diffusion model, Comput. Methods Appl. Mech. Engrg. 273 (2014) 19–36. [27] X. Chen, M. Gunzburger, Continuous and discontinuous finite element methods for a peridynamics model of mechanics, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1237–1250. [28] Q. Du, L. Ju, L. Tian, K. Zhou, A posteriori error analysis of finite element method for linear nonlocal diffusion and peridynamic models, Math. Comp. 82 (2013) 1889–1922. [29] H. Wang, H. Tian, A fast Galerkin method with efficient matrix assembly and storage for a peridynamic model, J. Comput. Phys. 231 (2012) 7730–7738. [30] K. Zhou, Q. Du, Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary condition, SIAM J. Numer. Anal. 48 (2010) 1759–1780. [31] S.F. Henke, S. Shanbhag, Mesh sensitivity in peridynamic simulations, Comput. Phys. Comm. 185 (2014) 181–193. [32] M.L. Parks, R.B. Lehoucq, S.J. Plimpton, S.A. Silling, Implementing peridynamics within a molecular dynamics code, Comput. Phys. Comm. 179 (2008) 777–783. [33] M.L. Parks, P. Seleson, S.J. Plimpton, S.A. Silling, R.B. Lehoucq, Peridynamics with LAMMPS: A User Guide v0.3 Beta, SAND Report 2011-8523, Sandia National Laboratories, Albuquerque, NM, Livermore, CA, 2011. [34] K. Yu, X.J. Xin, K.B. Lease, A new adaptive integration method for the peridynamic theory, Modelling Simul. Mater. Sci. Eng. 19 (2011) 045003. [35] S. Silling, Linearized theory of peridynamic states, J. Elasticity 99 (2010) 85–111. [36] P.J. Davis, Circulant Matrices, Wiley-Intersciences, New York, 1979. [37] R.M. Gray, Toeplitz and circulant matrices: a review, Found. Trends Commun. 2 (2006) 155–239.

X. Zhang, H. Wang / Comput. Methods Appl. Mech. Engrg. 311 (2016) 280–303

303

[38] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, Rhode Island, 2003. [39] H. Wang, T.S. Basu, A fast finite difference method for two-dimensional space-fractional diffusion equations, SIAM J. Sci. Comput. 34 (2012) A2444–A2458. [40] G. Evangelatos, P. Spanos, A collocation approach for spartial discretization of stochastic peridynamic modeling of fracture, J. Mech. Mater. Struct. 6 (2011) 1171–1195. [41] B. Kilic, E. Madenci, Structural stability and failure analysis using peridynamic theory, Int. J. Non-Linear Mech. 44 (2009) 845–854. [42] P. Seleson, M.L. Parks, M. Gunzburger, R.B. Lehoucq, Peridynamics as an upscaling of molecular dynamics, Multiscale Model. Simul. 8 (2009) 204–227. [43] X. Zhang, M. Gunzburger, L. Ju, Nodal-type collocation methods for hypersingular integral equations and nonlocal diffusion problems, Comput. Methods Appl. Mech. Engrg. 299 (2016) 401–420. [44] S. Chen, F. Liu, X. Jiang, I. Turner, V. Anh, A fast semi-implicit difference method for a nonlinear two-sided space-fractional diffusion equation with variable diffusivity coefficients, Appl. Math. Comput. 257 (2015) 591–601. [45] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, J. Comput. Appl. Math. 166 (2004) 209–219. [46] H. Wang, K. Wang, T. Sircar, A direct O(N log2 N ) finite difference method for fractional diffusion equations, J. Comput. Phys. 229 (2010) 8095–8104. [47] Q. Yang, I. Turner, F. Liu, M. Ilis, Novel numerical methods for solving the time-space fractional diffusion equation in 2D, SIAM Sci. Comput. 33 (2011) 1159–1180. [48] Q. Yang, I. Turner, T. Moroney, F. Liu, A finite volume scheme with preconditioned lanczos method for two-dimensional space-fractional reaction–diffusion equations, Appl. Math. Model. 38 (2014) 3755–3762. [49] T. Mengasha, Q. Du, Nonlocal constrained value problems for a linear peridynamic Navier equation, J. Elasticity 116 (2014) 27–51.