Particle nite element method (PFEM) is an e ective numerical tool for solv- ing large-deformation problems in geomechanics. By incorporating the node integration technique with strain smoothing into the PFEM, we proposed the smoothed particle nite element method (SPFEM). This paper extends the SPFEM to three-dimensional cases and presents a SPFEM executed on graphics processing units (GPUs) to boost the computational eciency. The detailed parallel computing strategy on GPU is introduced. New compu- tation formulations related to the strain smoothing technique are proposed to save memory space in the GPU parallel computing. Several benchmark problems are solved to validate the proposed approach and to evaluate the GPU acceleration performance. Numerical examples show that with the new formulations not only the memory space can be saved but also the computa- tional eciency is improved. The computational cost is reduced by 70% for the double-precision GPU parallel computing with the new formulations. Corresponding author. E-mail: wei.wu@boku.ac.at Preprint submitted to Computers and Geotechnics October 14, 2020This is a preprint of the paper published in Computers and Geotechnics Accepted citation: Zhang, W., Zhong, Z.H., Peng, C., Yuan, W.H. and Wu, W., 2021. GPU-accelerated smoothed particle finite element method for large deformation analysis in geomechanics. Computers and Geotechnics, 129, p.103856. DOI: https://doi.org/10.1016/j.compgeo.2020.103856 Compared with the sequential CPU simulation, the GPU-accelerated sim- ulation results in a signi cant speedup. The overall speedup ranges from 8.21 to 11.17 for double-precision simulations. Furthermore, the capability of the GPU-accelerated SPFEM in solving large-scale complicated problems is demonstrated by modelling the progressive failure of a long slope with strain-softening soil. Keywords: Particle nite element method; Parallel computing; GPU; Speedup; Geomechanics; Large deformation 1. Introduction Many problems in geotechnical engineering involve large deformation, e.g., sampling, in-situ penetration testing, pile installation, landslides and debris ow, etc. Large-deformation numerical simulations can provide a cru- cial perspective for these problems. However, they often pose some diculties for conventional numerical method, e.g. the nite element method (FEM). The mesh distortion associated with large deformation may lead to reduction of accuracy or even termination of the numerical calculations. Many numerical frameworks have been presented to solve large-deformation problems in geomechanics so far, such as Arbitrary Lagrangian-Eulerian (ALE) method [1], remeshing and interpolation technique with small strain (RITSS) method [2], Smoothed particle hydrodynamics (SPH) [3], mate- rial point method (MPM) [4], particle nite element method (PFEM) [5], etc. Among them, the PFEM attracts more and more attentions in re- cent years. It inherits the solid theoretical foundation of FEM which can guarantee the accuracy and convergence, whilst it possesses the exibility of 2 mesh-free methods for large deformation analysis. In the PFEM, the nodes in traditional FEM are treated as Lagrangian particles. The Delaunay tri- angulation technique is used to re-generate the mesh when it is too distorted due to large deformation. Therefore, the PFEM is actually an updated La- grangian FEM approach with remeshing technique to overcome the mesh distortion problem. In the eld of geomechanics, the PFEM has been devel- oped and applied to solve many di erent types of problems, including the ground excavation process [6, 7], granular ow [8{10], landslide and landslide- generated wave [11{13], soil penetration problem [14, 15], hydro-mechanical coupled problem [15{17], and progressive slope failure [18{20]. In the original PFEM, the numerical integration is carried out at Gasuss points while the eld information is stored at nodes/particles, As a result, during the numerical computation, it is required to transfer information be- tween Gasuss points and nodes/particles frequently, which inevitably intro- duces errors and increases complexity. In view of this, the authors proposed the smoothed particle nite element method (SPFEM) recently [21, 22]. In the SPFEM, a strain smoothing technique for node integration is incorpo- rated. Thanks to this node integration technique, all the eld variables are computed and stored at nodes/particles, making the method more like a La- grangian particle-based method. There are also some extra advantages, such as utilization of low order elements without volumetric locking and insensi- tivity to mesh distortion. The former can assure the computational eciency while the latter is especially bene cial for large deformation analysis. Very recently, the SPFEM has been extended to solve the uid dynamics prob- lem [23] and uid{structure interaction problem [24], which has been termed 3 'PFEM with Nodal Integration' in their studies [25]. In recent years, using graphic processing units (GPUs) to achieve high- performance computing has gained popularity rapidly. GPUs were initially designed for three-dimensional image rendering that is not complex but highly computationally intensive. In 2007, NVIDIA released the compute uni ed device architecture (CUDA) programming platform which is designed for general-purpose computing. Since then, GPUs have been successfully used in many scienti c elds for high-performance computing (e.g. [26{29]). GPUs, which feature much more cores, lower thread-scheduling cost and higher memory bandwidth than ordinary CPUs, are capable of executing computations with thousands of independent threads, operating in Single- Instruction-Multiple-Data (SIMD) mode. Here, SIMD is a terminology used to describe a class of parallel computers with multiple processing elements that perform the same operation on multiple data points simultaneously. In the computational process of PFEM/SPFEM, there are many procedures in which the same calculations need to be performed for all elements or nodes/particles. Meanwhile, data dependence is local and usually only ex- tends to a few surrounding elements or nodes. Therefore, most calculations can be distributed into independent GPU threads and executed in paral- lel. Compared with the PFEM, the computation process of SPFEM is much simpli ed due to the adoption of node integration scheme. The computation process of the explicit SPFEM [22] is especially simple because the accumu- lation of global sti ness matrix in the implicit SPFEM is avoided. These simpli cations can undoubtedly greatly facilitate GPU parallel computing. However, to the best of the authors' knowledge, GPU parallel computing for 4 PFEM has not been reported yet. This study aims to develop a GPU-accelerated SPFEM for large defor- mation analysis in geomechanics based on CUDA, with focus on the explicit version of SPFEM [22]. The explicit SPFEM is rst extended to three- dimensional cases. Then the GPU parallel computing strategy is presented after a brief view on the basic theory of SPFEM. To reduce the memory con- sumption during the GPU parallel computing, which is crucial for increasing the scale of problems analysed, new computation formulations related to the strain smoothing technique are used. Several benchmark examples are solved with sequential CPU simulations and GPU-accelerated simulations respec- tively. The speedup of the presented GPU-accelerated SPFEM is carefully evaluated. 2. Smoothed particle nite element method 2.1. Strain smoothing technique for node integration In PFEM, the problem domain is discretized into a set of Lagrangian particles, which carry mass, velocity and other state variables (e.g. stress, hardening parameters related to constitutive models). The FEM computa- tion mesh is generated based on this set of Lagrangian particles using the Delaunay triangulation technique combined with the alpha shape algorithm. Thus, nodes in traditional FEM correspond to particles in PFEM. The terms 'node' and 'particle' can be used interchangeably in PFEM. In SPFEM, the strain smoothing technique for node integration is in- corporated. The rst step of this technique is to construct smoothing cells associated with particles. In two-dimensional cases, the smoothing cell as- 5 sociated with particle kis created by connecting sequentially the mid-edge points to the centroid points of the surrounding triangular elements of the particle. This approach can be easily extended to three-dimensional cases, as shown in Fig. 1. Note that it is not required to explicitly calculate the speci c geometries during the computation process, which will be shown in the following. With this approach, the whole problem domain is discretized into nite strain smoothing cells associated with particles. Thus, all the integration over the whole problem domain can be performed by accumulating over the smoothing cells associated with particles. Furthermore, all the eld variables can be calculated in these smoothing cells. In SPFEM, each particle repre- sents a smoothing cell associated with it, which makes the SPFEM more like a particle-based method. In the smoothed strain technique, the volume of the smoothing cell asso- ciated with node kis calculated by Vk=EkX i=11 4Vi(1) whereEkis the number of tetrahedral elements related to node kandViis the volume of the ith element. The smoothed strain matrix is calculated by ~Bk=1 VkEkX i=11 4ViBi(2) whereBiis the strain gradient matrix used in the standard FEM for the ith element. It can be seen that the smoothed strain is actually the average of 6 the strains of the elements related to node kweighted by volume. For more details about the strain smoothing technique, readers are suggested to refer to [31, 32]. 2.2. Computational formulations of SPFEM Similar to the FEM, there are two numerical solution strategies of SPFEM with di erent time integration schemes: the implicit SPFEM [21] and the ex- plicit SPFEM [22]. In this study, the explicit SPFEM is considered, which has simpler formulations and thus is more suitable for GPU parallel comput- ing. In the computation domain, the motion of a continuum can be described as a=r+b (3) where,is the material density, ais the acceleration, is the Cauchy stress tensor andbis the speci c body force density. Considering the principle of virtual displacement and the divergence theorem, the weak form is expressed asZ uad =Z Su SdS+Z ubd Z u:d (4) whereuis the displacement vector, represents the con guration domain, Srepresents the boundary, and Sis the prescribed traction. After the node-based discretization, the above equation reads TX k=1akVk=TX k=1Z SNk SdS+TX k=1bkVkTX k=1~BT kkVk (5) whereTis the total number of nodes in the computation domain. The above 7 Algorithm 1 Computational cycle of SPFEM 1. Generate mesh using Delaunay triangulation and alpha shape technique 2. Get basic data of elements and nodes 3. Calculate smoothed strains of nodes 4. Update stresses of nodes through constitutive integration 5. Calculate internal forces of nodes 6. Update velocities and positions of nodes discrete form can be then written in the vector or matrix form as Ma =FextFint(6) in which Fext=TX k=1Z SNk SdS+TX k=1bkVk (7) Fint=TX k=1~BT kkVk (8) M=TX k=1Vk (9) whereFextandFintare termed as the external and internal forces, respec- tively, andMis the diagonal mass matrix. A typical computational cycle of SPFEM is shown in Algorithm 1. More details are available in Ref. [21, 22]. 2.3. Adaptive critical time step size In utilizing the explicit time integration scheme, the incremental time step size should satisfy the Courant stability restriction to ensure numerical 8 stability. The incremental time step size is generally calculated by tcr= lmin c(10) where tcris the critical time step size, lminis the characteristic length of the mesh,cis the sound speed in the material and is a scale factor. As the characteristic length of the mesh varies during the computation, an adaptive time step size is realized. The determination of time step sizes greatly in uences the computation cost of a numerical method using explicit time integration. In the original PFEM, the determination of time step sizes is the same as the standard FEM. For example, the characteristic length of the mesh in Eq. (10) can be taken to be the minimum radius of spheres inscribed in the tetrahedral elements [33]. As stated in [34], the Delaunay triangulation and the triangulation satisfying the max-min angle criterion are identical in two-dimensional cases, while in three-dimensional cases they are no longer identical in general and the Delaunay triangulation in three-dimensional cases does not seem to satisfy any optimal angle condition. Actually, many so-called silver elements may occur after Delaunay triangulation [35]. As a result, with the criterion of the minimum radius of spheres inscribed in the tetrahedrons, the time step size may reduce by several orders of magnitude during the computation [33]. In this study, as the node integration technique of strain smoothing is adopted, the integration domain for each node is the smoothing cell asso- ciated with this node. The smoothing cells generally have a more regular shape (See Fig. 1) and their volumes are generally much larger than those of the silver elements, because their volumes are the summation of the volumes 9 of elements related to the nodes divided by four (See Eq. (1)). The charac- teristic length of a smoothing cell can be obtained in analogy to an element in the standard FEM. In this study, one half of the minimum distance be- tween the node and its related nodes is used. This value is generally much larger than the value determined by the criterion of the minimum radius of spheres inscribed in the tetrahedrons, and thus the computational eciency can be greatly increased as a result. This can be considered as a newly dis- covered advantage of incorporating the strain smoothing technique for node integration into an explicit PFEM. 2.4. Contact strategies in the explicit SPFEM In the explicit SPFEM, as the solution strategy is similar to that of the explicit FEM, the contact algorithm in the explicit FEM can be easily in- corporated. Speci cally, the slave node to master segment contact scheme is utilized. A simple illustration is shown in Fig. 2. The frictional contact condition can be written as: gn0;n0;ngn= 0 jtjfn0(11) wheregnis the gap between node and segment, nis the normal contact force which is positive corresponding to compression, tis the tangential contact force, andfis the friction coecient. The penalty method is utilized to ful l the frictional contact condition. For the explicit time integration scheme, when a contact pair is detected, i.e. 10 gn<0, the contact force can be calculated as Rcs=; for slave node Rcm;i=Ni; for master nodes(12) whereis the vector of contact force, iis the node number of contact segment andNiis the corresponding shape function of contact segment. For the vector of contact force , the normal part nis obtained as n=gnwhile the tangential part tis obtained as t=min(gt;n) when the Coulomb friction model is used. The above contact forces are added into Fextin Equation (6) and then the accelerations, velocities and positions of nodes can be updated correspondingly. In case that the node to rigid wall contact is considered, the numerical implementation can be further simpli ed. When the penetration between nodekand a rigid wall occurs, i.e. gn<0, the acceleration and velocity of nodekcan be modi ed as follows to make the normal penetration gnto be zero, anew=aold(aoldn)n vnew=vold(voldn)n(13) wherenis the unit normal vector of the rigid wall. When the Coulomb friction model is used, the tangential contact force can be obtained as t=min(mkpvnewvnew=t; n) (14) in whichmkis the mass of node k, tis the time increment, and the normal 11 contact force ncan be calculated by n=mkp vn oldvn old=t (15) where vn oldis the normal velocity of node kbefore modi cation. 2.5. Mesh rebuilding technique in the explicit SPFEM An important part of SPFEM is the mesh rebuilding technique to avoid mesh distortion when large deformations occur in the continuum domain. The mesh rebuilding technique in the explicit SPFEM is basically the same with that in the original PFEM [5], which is generally based on the Delaunay triangulation technique combined with alpha shape algorithm. When the mesh rebuilding is required, the Delaunay triangulation tech- nique, which is widely used in FEM mesh generation [36], is rst performed to rebuild the element connectivity on the basis of a cloud of nodes. Another task of mesh rebuilding is the identi cation of computational domain. There is generally no unique solution to this problem. The solution originally pro- posed by Idelsohn [37] and subsequently adopted as a standard feature of the PFEM is to use the alpha-shape method previously developed by Edels- brunner [38] for computer graphics applications. The basic principle of the alpha-shape method is that for a cloud of points with a characteristic spacing hand a prede ned value of parameter , all nodes on an empty sphere with a radius greater than hare considered as boundary nodes. When the alpha-shape method is combined with the Delaunay triangu- lation technique, a simple two-step algorithm for mesh rebuilding is avail- able [39]: First, Delaunay triangulation is performed to generate the convex 12 domain on the basis of a cloud of nodes, and thus a mesh consisting of tetrahedrons is obtained in three-dimensional case. Then, the radius of the circumsphere of each tetrahedron is examined and the tetrahedrons with a radius greater than hare deleted, and hence the remain tetrahedrons cor- respond to the computational domain. Note that in the SPFEM, thanks to the node integration technique of strain smoothing, more distorted mesh can be used without serious loss of accuracy, i.e. the requirement of mesh quality is not as strict as the original PFEM. Nevertheless, mesh smoothing technique is utilized here to improve the mesh quality. To be speci c, the Laplacian smoothing technique [40, 41] is applied to a few nodes when they are too close to their neighbours. Some other mesh smoothing technique can also be used (e.g. [33, 42]), however, we choose the Laplacian smoothing technique due to its simplicity and high eciency. Moreover, with the Laplacian smoothing technique, the critical time step size can be increased remarkably, which is useful for improving the computational eciency. It should be pointed out that in some PFEM simulations, new particles were introduced during the simulation to further improve the mesh qual- ity (e.g. [14, 43, 44]). However, for the present approach, it seems that introducing new particles is not necessary and thus no new particle is intro- duced during the simulation for simplicity. 13 3. GPU parallel computing strategy 3.1. GPU acceleration scheme The GPU parallel computing is implemented based on the Compute Uni- ed Device Architecture (CUDA), a toolkit developed by NVIDIA to perform parallel computing using their graphics cards. The parallelisation scheme is shown in Fig. 3. The whole computation consists of a CPU part and a GPU part. The CPU part works as a process controller, which is used for loading and saving data, initializing GPU data, and controlling the GPU computing kernels to nish designated computation tasks, while the GPU part executes the actual computing. Ideally, all the steps are expected to be parallelised on GPU to maximise the speedup of GPU parallel computing. However, in the SPFEM, the De- launay triangulation technique is required to re-generate the computation mesh. In this study, a state-of-the-art library for Delaunay triangulation, CGAL [45], is adopted to obtain a robust Delaunay triangulation. So the Delaunay triangulation is executed on CPU. However, the workload of De- launay triangulation is minimal, because it is performed intermittently only when the mesh is too distorted. Hence the CPU execution of the Delaunay triangulation has only a small in uence on the overall speedup. Except for the Delaunay triangulation, all other calculations are executed on GPU in a parallelised manner. To reduce memory transfer as much as possible, es- pecially the transfer between CPU and GPU, all the data are stored in the GPU global memory during the entire simulation, and they are transferred to the CPU memory only when the Delaunay triangulation is required and when the occasional data output is required. 14 It can be seen from Algorithm 1 that the computation cycle of SPFEM is generally quite simple, which greatly ensures the applicability of GPU parallel computing. Especially for Steps 2, Step 4 and Step 6 in Algorithm 1, the GPU parallel computing is straightforward. Because in these steps, the computations at nodes are independent to each other, they can be executed on GPU by parallelising over nodes with no race condition. The so-called race condition is de ned as the situation that 2 or more threads try to change the same shared or global memory location at the same time. However, for Step 3 and Step 5, the smoothed strain and the internal force contribution of a node are related to its adjacent nodes. As a result, the race condition occurs, which should be paid more attention to. Note that Step 3 and Step 5 are both related to the strain smoothing technique. 3.2. Parallelisation of the strain smoothing technique A crucial feature of the SPFEM is the adoption of strain smoothing tech- nique. From a computation perspective, the core of this technique lies in the calculation of the smoothed strain matrix ~Bfor all nodes. It is the premise of the calculation of smoothed strains and internal forces of nodes. In three- dimensional cases, the number of related nodes for a node is much higher than that in two-dimensional cases, which leads to a big smoothed strain matrix for each node and greatly increases the memory requirement during the com- putation of SPFEM. For example, in the kernel function for calculating the internal forces of nodes, we need either to calculate the smoothed strain ma- trix~Bor to read it calculated before. In both approaches, the smoothed strain matrices can be quite memory-consuming, which is a bottleneck for the application to large-scale problems. 15 To reduce the GPU memory consumption related to the strain smoothing technique, new formulations are used in this study. The meaning of reducing memory consumption is twofold. First, the main function of GPU parallel computing is to solve large-scale problems. Reducing memory consumption can increase the scale (i.e. number of nodes) of the problem to be simulated with the same hardware. Please note that the memory of a GPU is constant and not easy to be extended as a CPU. Meanwhile, reducing memory con- sumption is also meaningful for improving eciency. As the global memory used to store the smoothed strain matrices is a high-latency memory, accesses to the global memory are time-consuming and should be reduced as much as possible to improve the eciency. Let's start from the original formulation. The smoothed strain of a node kis calculated by ~"k=~Buk (16) whereukis the vector assembled by the displacements of all the nodes in the smoothing cell associated with node k. By substituting Eq. (2) into Eq. (16), the smoothed strain of the node kis calculated as ~"k= (1 VkEkX i=11 4ViBi)uk (17) When calculating the smoothed strains of nodes with the original formu- lation (Eq. (17)), parallelization over nodes is performed. In the thread of a speci c node k, looping over the elements related to this node is performed. Noting the equivalence between the vector assembled by the displace- ments of the nodes related to node kand that assembled by the displace- 16 ments of the elements related to node k,the above equation can be re-written as ~"k=1 VkEkX j=11 4VjBjuj(18) whereEkis the number of elements related to node k,Vjis the volume of elementj,Bjis the strain matrix of element jin standard FEM, and uj is the vector assembled by the displacements of element j.Please note the di erence between the original formulation (Eq. (17)) and the new formula- tion (Eq. (18)). With the new formulation (Eq. (18)), parallelization over elements can be easily performed to obtain the smoothed strains of all the nodes. The new formulation for the calculation of internal forces related to the strain smoothing technique is a bit more complicated. In the original for- mulation of SPFEM, the internal forces applied on its neighbour nodes for a nodekcan be calculated as, Fint k=NkX i=1~BT kkVk (19) whereNkis the number of the neighbour nodes of node k. By substituting Eq. (2) into Eq. (19), we have Fint k=NkX i=1EkX j=11 4VjBjk (20) whereEkis the number of elements related to node k. Since the two sum- mations in Eq. (20) are both carried out on the elements related to node k, 17 Eq. (20) can be re-written as Fint k=EkX j=14X m=11 4VjBj mk (21) in whichBj mis the strain matrix related to node mof elementj. Let us make a comparison between the original formulations (Eqs. (17) and (19)) and the new ones (Eqs. (18) and (21)) from the perspective of GPU memory consumption. In three-dimensional cases, the strain matrix of nodekrelated to a speci c node consists of 18 oating-point numbers. Hence, using the original formulation, a total of ( Nk+ 1)18 oating-point numbers should be stored to obtain the smoothed strain matrix of a partic- ular nodek. As the number of the neighbour nodes Nkof a speci c node is undetermined, it is sensible to allocate a memory of equal size which is slightly larger than needed. As stated above, the possible maximum value of Nkin three-dimensional cases is much larger than that in two-dimensional cases. In this study, Nkis taken to be 64 to guarantee robustness. Bearing in mind that the smoothed strains and internal forces of all the nodes are calculated in a parallelised way. As a result, M(Nk+1)18 oating-point numbers should be saved during the computation, where Mis the concur- rent thread number in GPU, which generally ranges from several thousands to tens of thousands for a common GPU. One may save memory by only saving the derivatives of shape functions instead of the strain matrix. This approach will surely reduce the memory requirement, but it is not an es- sential improvement. However, with the new formulations, if we calculate and store the strain matrices of all the elements before, which is sensible to 18 avoid unnecessarily repeated computation, only the strain matrix variable Bj mis required to be loaded. Thus, in the kernel function for calculating the smoothed strains of nodes and the kernel function for calculating the inter- nal forces of nodes, the memory occupation related to the strain matrices is greatly reduced from M(Nk+ 1)18 oating-point numbers to M18 oating-point numbers. The memory occupation can be further reduced by only saving the derivatives of shape functions. Note that as parallel computing is performed with GPU threads, the so- called race condition may occur. With the new formulations, the smoothed strain of a particular node may be changed by all the threads of the elements related to it at the same time, and the internal force of a particular node may also be changed by all the threads corresponding to its neighbour nodes at the same time. So race condition occurs for both the calculations of the smoothed strains and the internal forces. With the original formulations, the former race condition does not occur while the latter race condition exists. To deal with the race condition problem, the so-called atomic operations, which are provided by CUDA, are used to calculate the smoothed strains and the internal forces of nodes. With the atomic operations, we are capable of reading, modifying, and writing a value back to the same memory address while avoid interferences between di erent threads. Please note that the time cost of atomic operations is generally notable. However, in the present approach, because atomic operations are only per- formed to accumulate the smoothed strains and the internal forces of nodes whereas the intensive computations of the smoothed strains and the inter- nal forces are still parallelized, the time cost of waiting in atomic operations 19 is not so notable. Taking the computation of the smoothed strains as an example, although atomic operations are performed with the new formula- tion, the time cost is still slightly lower than that with the original formula- tion. The main reason is that for each thread the computational workload is much heavier with the original form although atomic operations are avoided. With the original formulation, parallelization over nodes is performed. In the thread of a speci c node k, looping over the elements related to node k is required. Moreover, the identi cation of the sequential number of node k in its surrounding elements is also required. On the contrary, with the new formulation, parallelization over elements is performed with a rather simple procedure for each thread. 3.3. Parallelisation for getting the indices of elements and nodes related to nodes Another crucial step in the SPFEM is to get the indices of elements and nodes related to all the nodes, although this step is not so time-consuming in the entire simulation. This step should be performed in Step 2 in Algorithm 1, before the calculation of smoothed strains. A naive method is to loop over all the elements for each node to determine if the element is related to this node. However, this method is both time-consuming and memory- consuming. A more ecient algorithm is used instead in this study. As shown in Fig. 4, a kernel function in GPU is executed to perform the same calculation for all the elements with parallelised GPU threads. For the elementi, its indexiis recorded respectively to the four nodes of the element. A key step is to accumulate the serial number of current element related to each node. The so-called atomic operations are used again to do this. As 20 the sequence of the related elements for each node does not matter in the SPFEM, desired results can be obtained by the parallelisation over elements with GPU threads. It is more straightforward to get the indices of nodes related to nodes once the indices of elements related to nodes are obtained. Parallelisation over nodes with GPU threads can be performed without the so-called race condition. For each node, all the nodes of the elements related to this node are picked out and the unique indices of nodes are records as the indices of nodes related to nodes. 4. Numerical examples To validate the proposed approach and to evaluate the speedup of the GPU parallel computing, three benchmark examples are numerically solved. The rst example is the vibration of an elastic continuum bar. This prob- lem has an analytical solution and is suitable for validating the proposed approach, especially for investigating the new formulations related to the strain smoothing technique in terms of correctness and eciency. In the sec- ond example, the laboratory experiment of soil collapse is simulated. The experiment was carried out by Bui [3] to validate their SPH approach for geo-materials. The third example is the progressive failure of a long strain- softening soil slope. Since the geometry of the problem is complicated during the whole simulation, and a strain softening constitutive model is used, it is suitable to validate the application capability of the proposed approach on an actual engineering problem. The GPU parallel computing was conducted on a personal computer with 21 an NVIDIA GeForce GTX 1080Ti GPU. On this card, 3584 cores are avail- able and the global memory is 11 GB. For comparison purpose, the rst and the second example were also solve by sequential CPU simulation. The CPU is Intel i7-8650U with a frequency of 1.90 GHz. Di erent meshes with various number of nodes are used for the rst and the second examples to evaluate the speedup thoroughly. In the rst example, both single-precision and double- precision simulations are performed for comparison. For other examples, only double-precision simulations are performed, because single-precision simula- tion is not suitable for complicated problems, such as the strain-softening [46] and other speci c behaviour of soil (e.g. [47{53]). Furthermore, the results of the rst example show that the eciency of single-precision simulation is only slightly higher than that of double-precision simulation. 4.1. Axial vibration of an elastic continuum bar The rst example considers the axial vibration of an elastic continuum bar. As shown in Fig. 5,the length, width and height of the continuum bar are 10 m, 1 m and 1 m respectively. The bar is linearly elastic and the material parameters are: Young's modulus E= 100 Pa, Poisson's ratio = 0 and the density = 1kg=m3. The left end (x = 0 m) of the bar is xed while all other surfaces of the bar are free. Free vibration of the bar along the axial direction (i.e. x-direction) is considered. This problem may be solved using separation of variables [54]. The bar is found to oscillate in modes that are dependent on the initial conditions. A speci c mode can be excited by using initial conditions which are multiples of this natural mode of the system. Here, the rst mode is considered, and 22 thus the analytical solution of this problem is given by [54, 55]: v(x;t) =v0cos(!1t)sin( 1x) u(x;t) =v0sin(!1t)sin( 1x)=!1(22) wherev(x;t) is the velocity at time tandu(x;t) is the displacement at time t. The eigenvalue is 1==2L, and the frequency of oscillation related to the eigenvalue is !1= 1cwherec=p E= is the speed of elastic wave. The corresponding initial velocity and displacements are given by v(x;0) =v0sin( 1x) u(x;0) = 0(23) A total of 5 meshes with di erent numbers of nodes (12221, 88641, 289261, 674081, 1303101) are utilized in the simulations. For each mesh, both se- quential CPU simulation and GPU accelerated simulation are conducted, and both single-precision and double-precision simulations are performed. Moreover, the original formulations and the new formulations related to the strain smoothing technique are used respectively. A total time of 10.0 s is simulated. The incremental time step size remains constant during each sim- ulation. However, for di erent mesh density, the incremental time step size is di erent to maintain numerical stability. Correct results are obtained with all the simulations. The numerical results obtained by the single-precision GPU simulation with the coarsest mesh and the new formulations are compared with the analytical solution, as shown in Fig. 6. They are almost identical to each other. 23 Figure 7shows the comparison of the computational time cost spent with the original formulations and with the new ones. The time cost ratio, which is de ned as the ratio between the time cost with the original formulations and that with the new formulations, is de ned here to make a clear compar- ison. It can be seen that less time cost is spent with the new formulation during all the simulations. The ratio between the time cost with the original formulation and that with the new formulation ranges from 1.71 to 2.54 for single-precision GPU computation and ranges from 3.69 to 4.20 for double- precision GPU computation. It is clear that adopting the new formulation can not only save memory occupation but also greatly improve computa- tional eciency. Therefore, the new formulations are adopted throughout all the following simulations due to their superiority. The computational costs for di erent meshes with di erent simulation measures are listed in Table 1. The double-precision CPU simulation with the largest number of nodes is taken as a reference to evaluate the distribution of the workload across the steps of the present approach. As shown in Fig. 8, for this example, most of the computational e ort is allocated to Step 3 and Step 5 (please see Algorithm 1), which respectively account for 29.89% and 64.87% of the total workload. The computational e ort of Step 1 and Step 2 is minimal because no re-meshing is performed. The computational e ort of Step 4 is low because a simple elastic constitutive model is used. The speedup ratios between the GPU simulation and CPU simulation are plotted in Fig. 9. The overall speedup ranges from 7.82 to 10.89 for single- precision computation and ranges from 8.40 to 12.30 for double-precision computation. It is clear that with the GPU-accelerated approach, a speedup 24 around 10 can be achieved. Moreover, the speedup with ner mesh is gen- erally slightly higher than that with coarser mesh. As shown in Fig 9, the performance of GPU simulation varies for di erent steps in Algorithm 1. The speedup of Step 4 is the largest because all the threads can be executed with- out mutual e ect. For Step 3 and Step 5, the speedup is smaller. However, the speedup performance is still rather obvious. The speedup for Step 6 is the lowest, which may be due to the algorithm complexity related to bound- ary condition application and the relatively low computation intensity. Note that since Step 1 is executed on the CPU, no speedup is recorded. Moreover, Step 2 only takes a minimal time even for the nest mesh, which increases the randomness of the value of speedup, so the speedup of Step 2 is also not plotted in Fig 9. 4.2. Laboratory experiment of soil collapse Bui et al. [3] conducted a simple laboratory experiment of soil col- lapse to validate their SPH approach for large-deformation analysis of geo- materials. Besides the SPH simulation [3], the experiment was also simulated by MPM [56] and two-dimensional SPFEM [22], so it can serve as a good benchmark. In the laboratory experiment, the soil continuum was modelled by a large number of small aluminium bars, which is a common approach to model cohesiveless soil physically. Two types of aluminium bars were used with a diameter of 1 mm and 1.5 mm respectively and a length of 50 mm. These aluminium bars were piled upon a at surface in parallel, and thus a plane-strain condition is mimicked. The soil continuum was modelled by arranging the aluminium bars into a rectangular area (200 mm x 100 mm) which is generated by standing two at solid walls on the at surface. The 25 soil collapse was triggered by quickly removing the right supporting wall. A reference grid was set behind the soil continuum, out of contact with the aluminium bars, to measure surface con guration of soil continuum after collapse. Moreover, square grids (20 mm x 20 mm) were plotted on the soil continuum to visualize the failure pattern. A critical problem is to determine the mechanical parameters of the soil continuum modelled by aluminium bars that can be used in the following numerical simulations. In the laboratory experiment, a shear box test appa- ratus was designed to perform tests on a group of aluminium bars, and then four experiments under di erent normal loading conditions were conducted to obtain the parameters of the soil continuum. A total of 5 meshes with di erent number of nodes (10999, 25535, 50338, 75720, 100072) are utilized for simulation (Fig. 10). For each mesh, double- precision GPU simulation and sequential CPU simulation are conducted re- spectively. The soil is treated as an elastic perfect plastic material with Mohr-Coulomb yield criterion. A zero dilation angle is assumed, correspond- ing to a non-associated ow rule. The following model parameters are used according to [3]: shear modulus G = 0.7 MPa, Poisson's ratio = 0.3, density = 2650kg=m3, cohesion strength c = 0 kPa and friction angle '= 19.8. Before the simulation of soil collapse, the initial stress due to gravity is gen- erated by a dynamic relaxation simulation, in which the bottom boundary is fully xed and the four side boundaries are constrained at their normal directions. Then the soil collapse is triggered by removing the restraint of the right boundary. A Coulomb contact model between soil and ground is used, in which no tangential relative motion is assumed. A total time of 1.0 26 s is simulated, after which a stable slope is formed. Adaptive time step size is utilized during the simulation, and the scale factor in Eq. (10) is taken to be 0.9. The nal surface con gurations after collapse obtained with di erent mesh densities are compared with the experimental result, as shown in Fig. 11(a). All numerical results show a similar response, and they all match well with the experimental results. With the increase of the number of nodes, the runout distance increases slightly. The failure line, which is de ned as the line between motionless soil and the moved soil, is also compared against the experimental result and good agreement can also be found. As shown in Fig. 11(b), the numerical results are compared with the results obtained by other numerical methods (i.e. SPH [3] and MPM [56]), as well as two- dimensional SPFEM [22]. All the numerical methods can well capture the experimental result. The run-out distance obtained by this study is slightly shorter than the two-dimensional SPFEM [22]. The reason may be that the nodes in the three-dimensional case interact with more nodes in comparison to two-dimensional cases. The nal mesh con guration with the coarsest mesh is chosen to make a visual comparison with the experimental result, as shown in Fig. 12. It is clear that the three-dimensional SPFEM is capable of reproducing the experimental result. The computational costs of double-precision GPU simulation and sequen- tial CPU simulation are listed in Table 2. The latter with the largest number of nodes is taken as a reference to evaluate the distribution of the workload across the steps of the present approach. As shown in Fig. 13, the time- consuming steps are Step 2, Step 3 and Step 5, which respectively account 27 for 19.33%, 20.84% and 36.66% of the total workload, while the computa- tional e ort of Step 6 is minimal. The above workload distribution may be considered as a typical case for geotechnical applications. As shown in Fig. 14, an overall speedup ranging from 8.21 to 11.17 is achieved by GPU parallel computing. A speedup from 35.31 to 66.49 can be achieved for the constitutive integration associated with plastic problems, which is the highest in all the steps, because the parallel operations over nodes are independent of other nodes. The speedup of Step 6 is the lowest among all the steps, which is similar to the rst example, i.e. the axial vibra- tion of a continuum bar. Compared with the rst example, more Delaunay triangulation operations are performed in Step 1, which are not accelerated by GPU parallel computing. Meanwhile, the constitutive integration in Step 4 achieves a higher speedup. These two factors together lead to an overall speedup similar to the rst example. 4.3. progressive failure of a long strain-softening soil slope This example considers the progressive failure of a long strain-softening soil slope. As shown in Fig. 15, the slope is initially 25 m long and 5 m high and has a slope angle of 45. A width of 5 m is chosen to generate a three-dimensional computation domain. The soil is considered to be a strain-softening soil to mimic the progressive failure of a sensitive clay slope, which is commonly found in nature. A main feature of sensitive clay is the decrease of undrained shear strength with the increase of plastic strain. This feature can be captured by a simple strain-softening Tresca model. As shown in Fig. 16, the cohesion equals to its peak value ( cp) at the beginning and decreases proportionally to the equivalent plastic strain invariant de ned by 28 "=p 2=3ep:ep. whereepis the deviatoric part of the plastic strain tensor. When the invariant achieve a certain value ( "pr), the cohesion reduces to its residual value ( cr). This example, feature by the complicated geometry during the whole simulation process and the strain-softening constitutive model, has been well simulated by two-dimensional MPM [57], PFEM [18] and SPFEM [13]. For the long strain-softening soil slope, the failure mode, run-out dis- tance and retrogression distance are greatly in uenced by the degree of the strength reduction (i.e. cr=cp). As stated in Ref. [59], lower residual strength corresponds to more times of retrogression failure, longer run-out distance and longer retrogression distance. In this study, a relatively severe strength reduction (i.e. 0.2) is considered to generate multiple retrogression failure. It should be pointed out that the utilization of the strain-softening model may result in the solution becoming mesh dependent. This stems from the fact that the related boundary value problem is no longer elliptic in statics or hyperbolic in dynamics [58]. Special regularisation techniques, such as nonlocal theory or gradient plasticity, have to be incorporated to overcome this issue. As this study focuses on the GPU acceleration of SPFEM, the regularisation is not involved. Double-precision GPU simulation is conducted to validate the capability of the GPU-accelerated SPFEM on solving large-scale geotechnical problems. A total number of 250584 nodes are used in the simulation. The following model parameters are used: Young's modulus E = 1.0 MPa, Poisson's ratio = 0.33, density = 2000kg=m3, peak cohesion cp= 20 kPa, residual cohesioncr= 4 kPa, and the softening modulus H = -50.0 kPa. A Coulomb 29 contact model between soil and ground is used and the friction coecient is taken to be 0.3. As shown in Fig. 15, the bottom boundary is fully xed and the left side boundary is constrained at x-direction. All the motions along y-direction are constrained to mimic a plane strain response. The soil is rst assumed to be elastic and the initial stress due to gravity is generated by a dynamic relaxation simulation. After the initial stress is generated, the soil plasticity is considered. Since the slope can not maintain stability with these material parameters, the landslide is triggered. A total of 10 seconds are simulated to reproduce the whole progressive failure process. Adaptive time step size is utilized and the scale factor in Eq. (10) is taken to be 0.9. The simulation is completed in around 5.66 hours. The initial vertical stress is shown in Fig. 17(a). The development process of the progressive failure is shown in Fig. 17(b)-(g), with coloured contours representing the equivalent plastic strain invariant  ". Many shear bands occur during the whole simulation process due to the adoption of the strain-softening soil con- stitutive model. As shown in Fig. 17(b), the rst major shear band (MS1) initiate from the bottom and propagate towards the top surface. During the sliding, two crossed shear bands (S1, S2) occur in the middle of the landslide body. One (S1) propagates towards the top surface and the other (S2) prop- agates towards the front inclined surface, resulting in a graben (Fig. 17(c)). With the advancement of the landslide body, a new major shear band (MS2) initiate from the bottom and propagate towards the top surface due to lack of support, leading to the second retrogressive failure (Fig. 17(d)). And then, a new shear band (S3) forms in the middle of the second landslide body, which initiates from the bottom and propagate towards the front inclined 30 surface (Fig. 17(e)). Then a new shear band (S4) initiates from the bottom of the landslide body and propagates towards the top surface (Fig. 17(f)). The retrogressive failure nally stops with the help of the friction between the landslide bodies and the ground. Fig. 17(g) shows the nal con guration of the slope. Eventually, the retrogressive failure results in a deposit with a run-out distance of 12.0 m and a retrogression distance of 13.6 m (Fig. 17(g)). A crucial factor a ecting the run-out distance and the retrogression dis- tance is the friction coecient between soil and ground. To investigate this, a new simulation is conducted with a friction coecient of 0.1. The com- parison of the nal con guration of the slope with di erent friction coef- cients is shown in Fig. 18. A smaller friction coecient corresponds to larger run-out distance and retrogression distance. The nal run-out dis- tance is 21.0 m while the retrogression distance is 18.4 m. Taken together, the GPU-accelerated SPFEM is capable of solving large-scale complicated large-deformation problem in geomechanics. 5. Conclusions PFEM is an attractive numerical framework for large deformation anal- ysis in geomechanics. By incorporating the strain smoothing technique for node integration, the authors proposed an improved PFEM, i.e. SPFEM. The primary advantages of SPFEM is that all the eld viables are computed and stored at nodes/particles. Some additional advantages include no volu- metric locking when low order elements are used, and insenstivity to element distortion. The objectives of the present work are twofold. First, the explicit SPFEM is extended to three dimensional cases. Then, a GPU accelerated 31 SPFEM is developed. A GPU parallelisation strategy for the three-dimensional explicit SPFEM is presented in detail. All the steps are executed on GPU except that the De- launay triangulation is executed on CPU with a state-of-the-art library. Since the numerical implementation of the operations related to strain smoothing technique is crucial to memory consumption, new formulations to accumu- late the smoothed strains and the internal forces of nodes are proposed. The computational time cost ratio between the original formulation and the new formulation ranges from 1.71 to 2.54 for single-precision GPU computation and from 3.69 to 4.20 for double-precision computation. Compared with the sequential CPU simulation, the GPU-accelearted simulation results in a signi cant speedup. The overall speedup ranges from 7.82 to 10.89 for single- precision computation for elastic problems and from 8.40 to 12.30 for double- precision computation. A speedup from 35.31 to 66.49 can be achieved for the constitutive integration associated with plastic problems, which is the highest in all types of calculations in the SPFEM. Finally, the developed GPU-accelerated SPFEM is employed to model the progressive failure of a long strain-softening soil slope. It is found that the presented method can model large-scale complex problems eciently and delivers reliable results. Acknowledgements The research is supported by the Natural Science Foundation of Guang- dong Province (Grant No. 2018A030310346), the H2020 Marie Sk lodowska- Curie Actions RISE 2017 HERCULES (778360) and FRAMED (734485), the Erasmus+ KA2 project Re-built (2018-1-RO01-KA203-049214), the Nazarbayev 32 University Research Fund (SOE2017001), and the Natural Science Founda- tion of China (Grant No. 41807223).