Extension of B-spline Material Point Method for unstructured triangular grids using Powell–Sabin splines

The Material Point Method (MPM) is a numerical technique that combines a ﬁxed Eulerian background grid and Lagrangian point masses to simulate materials which undergo large deformations. Within the original MPM, discontinuous gradients of the piecewise-linear basis functions lead to the so-called grid-crossing errors when particles cross element boundaries. Previous research has shown that B-spline MPM (BSMPM) is a viable alternative not only to MPM, but also to more advanced versions of the method that are designed to reduce the grid-crossing errors. In contrast to many other MPM-related methods, BSMPM has been used exclusively on structured rectangular domains, considerably limiting its range of applicability. In this paper, we present an extension of BSMPM to unstructured triangulations. The proposed approach combines MPM with C 1 -continuous high-order Powell–Sabin spline basis functions. Numerical results demonstrate the potential of these basis functions within MPM in terms of grid-crossing-error elimination and higher-order convergence.


Introduction
The Material Point Method (MPM) has proven to be successful in solving complex engineering problems that involve large deformations, multi-phase interactions and historydependent material behaviour. Over the years, MPM has been applied to a wide range of applications, including modelling of failure phenomena in single-and multi-phase media [1,35,40], crack growth [14,20] and snow and ice dynamics [12,31,34]. Within MPM, a continuum is discretised by defining a set of Lagrangian particles, called material points, which store all relevant material properties. An Eulerian background grid is adopted on which the equations of motion are solved in every time step. The solution on the background grid is then used to update all material point properties such as displacement, velocity and stress. In this way, MPM incorporates both Eulerian and Lagrangian descriptions. Similarly to other combined Eulerian-Lagrangian techniques, MPM attempts to avoid the numerical difficulties arising from non-linear convective terms associated with an Eulerian problem formulation, while at the same time preventing grid distortion, typically encountered within mesh-based Lagrangian formulations (e.g. [10,32]).
Classically, MPM uses piecewise-linear Lagrange basis functions, also known as 'tent' functions. However, the gradients of these basis functions are discontinuous at element boundaries. This leads to the so-called grid-crossing errors [2] when material points cross this discontinuity. Gridcrossing errors can significantly influence the quality of the numerical solution and may eventually lead to a lack of convergence (e.g. [28]). Different methods have been developed to mitigate the effect of grid-crossings. For example, the Generalised Interpolation Material Point (GIMP) [2] and Convected Particle Domain Interpolation (CPDI) [25] methods eliminate grid-crossing inaccuracies by introducing an alternative particle representation. The GIMP method represents material points by particle-characteristic functions and reduces to standard MPM, when the Dirac delta function centred at the material point position is selected as the characteristic function. For multivariate cases, a number of versions of the GIMP method are available such as the uniform GIMP (uGIMP) and contiguous-particle GIMP (cpGIMP) methods. The CPDI method extends GIMP in order to accurately capture shear distortion. Much research has been performed to further improve the accuracy of the CPDI approach and increase its range of applicability [16,22,26]. On the other hand, the Dual Domain Material Point (DDMP) method [41] preserves the original point-mass representation of the material points, but adjusts the gradients of the basis functions to avoid grid-crossing errors. The DDMP method replaces the gradients of the piecewise-linear Lagrange basis functions in standard MPM by smoother ones. The DDMP method with sub-points [7] proposes an alternative manner for numerical integration within the DDMP algorithm.
The B-spline Material Point Method (BSMPM) [28,29] solves the problem of grid-crossing errors completely by replacing piecewise-linear Lagrange basis functions with higher-order B-spline basis functions. B-spline and piecewise-linear Lagrange basis functions possess many common properties. For instance, they both satisfy the partition of unity, are non-negative and have a compact support. The main advantage of higher-order B-spline basis functions over piecewise-linear Lagrange basis function is, however, that they have at least C 0 -continuous gradients which preclude grid-crossing errors from the outset. Moreover, spline basis functions are known to provide higher accuracy per degree of freedom as compared to C 0 -finite elements [15]. On structured rectangular grids, adopting B-spline basis functions within MPM not only eliminates grid-crossing errors but also yields higher-order spatial convergence [3,11,28,29,37]. Previous research also demonstrates that BSMPM is a viable alternative to the GIMP, CPDI and DDMP methods [11,19,30,39]. While the CPDI and DDMP methods can be used on unstructured grids [16,22,41], to the best of our knowledge, BSMPM for unstructured grids does not yet exist. This implies that its applicability to real-world problems is limited compared to the CPDI and DDMP methods.
In this paper, we propose an extension of BSMPM to unstructured triangulations to combine the benefits of Bsplines with the geometric flexibility of triangular grids. The proposed method employs quadratic Powell-Sabin (PS) splines [23]. These splines are piecewise higher-order polynomials defined on a particular refinement of any given triangulation and are typically used in computer-aided geometric design and approximation theory [9,18,24,27]. PS-splines are C 1 -continuous and hence overcome the grid-crossing issue within MPM by design. We would like to remark that although this paper focuses on PS-splines, other options such as refinable C 1 splines [21] can be used to extend MPM to unstructured triangular grids. The proposed PS-MPM approach is analysed based on several benchmark problems.
The paper is organised as follows. In Sect. 2, the governing equations are provided together with the MPM solution strategy. In Sect. 3, the construction of PS-spline basis functions and their application within MPM is described. In Sect. 4, the obtained numerical results are presented. In Sect. 5, conclusions and recommendation are given.

Material Point Method
We will summarise the MPM as introduced by Sulsky et al. [32] to keep this work self-contained. First, the governing equations are presented, afterwards the MPM strategy to solve these equations is introduced.

Governing equations
The deformation of a continuum is modelled using the conservation of linear momentum and a material model. It should be noted that MPM can be implemented with a variety of material models that either use the rate of deformation (i.e. the symmetric part of the velocity gradient) or the deformation gradient. However, for this study it is sufficient to consider the simple linear-elastic and neo-Hookean models that are based on the deformation gradient. Using the Einstein summation convention, the system of equations in a Lagrangian frame of reference for each direction x k is given by where u k is the displacement, t is the time, v k is the velocity, ρ is the density, g k is the body force, D kl is the deformation gradient, δ kl is the Kronecker delta, σ kl is the stress tensor, x 0 is the position in the reference configuration, λ is the Lamé constant, E kl is the strain tensor defined as 1 2 (D kl + D lk ) − δ kl , J is the determinant of the deformation gradient and μ is the shear modulus. Equations 2 and 4 describe the conservation of linear momentum and the material model, respectively.
Initial conditions are required for the displacement, velocity and stress tensor. The boundary of the domain Ω can be divided into a part with a Dirichlet boundary condition for the displacement and a part with a Neumann boundary condition for the traction: where x = [ x 1 x 2 ] T and n is the unit vector normal to the boundary ∂Ω and pointing outwards. We remark that the domain Ω, its boundary ∂Ω = ∂Ω u ∪ ∂Ω τ and the normal unit vector n may all depend on time.
For solving the equations of motion in the Material Point Method, the conservation of linear momentum in Eq. 2 is needed in its weak form. For the weak form, Eq. 2 is first multiplied by a continuous test function φ that vanishes on ∂Ω u and is subsequently integrated over the domain Ω: in which a k = ∂v k ∂t is the acceleration. After applying integration by parts, the Gauss integration theorem and splitting the boundary into Dirichlet and Neumann part, the weak formulation becomes whereby the contribution on ∂Ω u equals zero, as the test function vanishes on this part of the boundary. Equation 8 can be solved using a finite element approach by defining n bf basis functions φ i (i = 1, . . . , n bf ). The acceleration field a k is then discretised as a linear combination of these basis functions:

Discretised equations
in whichâ k, j is the time-dependent j th acceleration coefficient corresponding to basis function φ j . Substituting Eq. 9 into Eq. 8 and expanding the test function in terms of the basis functions φ i leads to which holds for i = 1, . . . , n bf . By exchanging summation and integration, this can be rewritten in matrix-vector form as follows: where M denotes the mass matrix,â k the coefficient vector for the acceleration, while F trac k , F int k and F body k denote, respectively, the traction, internal and body force vector in the x k -direction. When density, stress and body force are known at time t, the coefficient vectorâ k is found from Eq. 12. The properties of the continuum at time t + Δt can then be determined using the acceleration field at time t.

Solution procedure
Within MPM, the continuum is discretised by a set of particles that store all its physical properties. At each time step, the particle information is projected onto a background grid, on which the momentum equation is solved. Particle properties and positions are updated according to this solution, as illustrated in Fig. 1.
To solve Eq. 12 at every time step, the integrals in Eq. 11 have to be evaluated. Since the material properties are only known at the particle positions, these positions are used as integration points for the numerical quadrature, and the particle volumes V are adopted as integration weights. Assume that the total number of particles is equal to n p . Here, subscript p corresponds to particle properties. A superscript t is assigned to particle properties that change over time. The mass matrix and force vectors are then defined by in which m p denotes the particle mass, which is set to remain constant over time, guaranteeing conservation of the total mass. The coefficient vectorâ t k at time t for the x k -direction is determined by solving Unless otherwise stated, Eq. 17 is solved adopting the consistent mass matrix. Next, the reconstructed acceleration field is used to update the particle velocity: The semi-implicit Euler-Cromer scheme [5] is adopted to update the remaining particle properties. First, the velocity field v t+Δt k is discretised as a linear combination of the same basis functions, similar to the acceleration field (Eq. 9): in whichv t+Δt k is the velocity coefficient vector for the velocity field at time t + Δt in the x k -direction. The coefficient vectorv t+Δt k is then determined from a density weighted L 2projection onto the basis functions φ i . This results in the following system of equations: where M t is the same mass matrix as defined in Eq. 13. P t+Δt k denotes the momentum vector with the coefficients given by P t+Δt Particle properties are subsequently updated in correspondence with Eqs. 1-4. First, the deformation gradient and its determinant are obtained: Here, I denotes the identify matrix and ε t+Δt p denotes the symmetric part of the velocity gradient: For linear-elastic materials, the particle stresses are computed as follows: For neo-Hookean materials, the stresses are obtained from The determinant of the deformation gradient is used to update the volume of each particle. Based on this volume, the density is updated in such way that the mass of each particle remains constant: Finally, particle positions and displacements are updated from the velocity field: The described MPM can be numerically implemented by performing the steps in Eqs. 13-31 in each time step in the shown order. Note that all steps can be applied with a variety of basis functions, without any essential difference in the algorithm. In this paper, we investigate the use of Powell-Sabin spline basis functions, which are described in the next section.

Powell-Sabin grid refinement
To construct PS-splines on an arbitrary triangulation, a Powell-Sabin refinement is required, dividing each of the original main-triangles into six sub-triangles as follows (see also Fig. 2).
1. For each triangle t j , choose an interior point Z j (e.g. the incentre), such that if triangles t j and t m have a common edge, the line between Z j and Z m intersects the edge. The intersection point is called Z jm . 2. Connect each Z j to the vertices of its triangle t j with straight lines. 3. Connect each Z j to all edge points Z jm with straight lines. In case t j is a boundary element, connect Z j to an arbitrary point on each boundary edge (e.g. the edge middle).
The area consisting of the main-triangles around V i is referred to as the molecule Ω i of V i (see Fig. 2). Each PSspline is associated with a main vertex V i and will be defined on the molecule Ω i .

Control triangles
For each main vertex V i , the set of PS-points is given by the union of V i and the sub-triangle edge midpoints directly around it, see Fig. 3. A control triangle for V i is then defined as a triangle that contains all its PS-points. Note that this control triangle is not uniquely defined. However, it preferably has a small area to ensure good stability properties of the resulting PS-splines [8].
A sufficiently good control triangle can be constructed by considering only control triangles with two or three edges shared with the convex hull of the PS-points and taking the one with smallest surface. Further details about the implementation of such an algorithm can be found in [6].
Each control triangle defines three basis functions, all associated with vertex V i and the molecule Ω i . Therefore, the basis functions are indexed φ The triplets of the three PS-spline basis functions corresponding to vertex V i are determined from the control triangle and the position of V i . Let the Cartesian coordinates of the control triangle vertices be 4 A control triangle lifted to z = 1 in one of its vertices for determining a triplet of a basis function. The location in the lifted plane above the vertex V i in the middle is marked with a dot; the value and gradient there correspond to the triplet of the associated basis function The three triplets are then determined by solving The control triangle has a direct geometric interpretation for the triplet. Solving for the triplet α

Bernstein-Bézier formulation for quadratic splines
i is constructed as a piecewise quadratic polynomial over each sub-triangle in the PS-refinement in barycentric coordinates. Cartesian coordinates (x, y) can be converted to barycentric coordinates (η 1 , η 2 , η 3 ) with respect to the considered triangle vertices Next, we define 6 quadratic Bernstein polynomials B i, j,k for each sub-triangle in barycentric coordinates, which are non-negative and form a partition of unity. Any desired quadratic polynomial b(η) can be uniquely constructed as a linear combination of the 6 quadratic Bernstein polynomials B i, j,k over a sub-triangle, where b i, j,k are called the Bézier ordinates. The desired quadratic polynomial can thus be fully defined by its Bézier ordinates, which can be schematically represented by associating Bézier ordinate b i, j,k with barycentric coordinates i 2 , j 2 , k 2 , as shown in Fig. 5.

Construction of PS-splines
Next, consider the molecule Ω i around a vertex V i and one of the main-triangles t j in the molecule, as depicted in Fig. 6. The grey main-triangle has been divided into 6 sub-triangles in the PS-refinement. In each of the sub-triangles, the 6 locations with barycentric coordinates i 2 , j 2 , k The Bézier ordinates for basis function φ 1 i , corresponding to the vertex V i considered in Fig. 6. There are three different basis functions, note that for each of these basis functions, α, β and γ should correspond with the triplet (α q , β q , γ q ) corresponding to PS-triangle vertex Q q with β = β(x 2 − x 1 ) + γ (y 2 − y 1 ) and The resulting PS-splines corresponding to the shown control triangle of vertex V i are shown on the associated molecule in Fig. 3.

Projecting on PS-splines and lumping within PS-MPM
In this section, we show by numerical experiments that the projection on PS-splines leads to third-order spatial convergence in the L 2 -error between an analytic function and its projection. Furthermore, we also investigate the effect of performing a projection on PS-splines using a lumped variant of the mass matrix instead of the consistent one. Lumping of the mass matrix is often used in MPM with piecewise-linear basis functions, but we will show that lumping is less suitable for MPM with PS-spline basis function when no special measures are taken. We consider the projection of a two-dimensional sine function onto a basis of PS-splines, using a consistent and a lumped mass matrix, respectively. The results are also com-pared with a projection onto a basis of piecewise-linear basis functions. We project the function f = sin(x/π ) sin(y/π ) on the unit square onto a basis of PS-splines using a standard L 2 -projection: The entries of the mass matrix M and the right-hand side b are evaluated using high-order Gauss integration such that the integration error is insignificant. Finally, we may obtain the projectionf by either solving M c = b consistently, or by using the lumped mass matrix, b i = c i /M ii , in which the lumped mass matrixM has been created by lumping all the mass of each row onto its diagonal element,M ii = n bf j=1 M i j ,M i j = 0 for i = j. A basis of PS-splines, and of piecewise-linear functions were adopted, respectively, both defined on a structured triangulation. Figure 8 shows the L 2 -error in the function value and xderivative for various refinements of the grid. When using a consistent mass matrix, the expected orders of convergence are obtained for piecewise-linear basis functions and PS-spline basis functions, both in the function value and the x-derivative. The use of a lumped mass matrix leads to a firstorder accurate approximation of the function value for both types of basis functions. However, the x-derivative does not converge at all when adopting a lumped mass matrix within a PS-spline basis. As the gradient of the reconstructed field is often used in MPM (e.g. Eq. 24), lumping within PS-MPM will often lead to spatial oscillations as will be shown in more detail in Sect. 4.3.

Numerical results
To validate the proposed PS-MPM, several benchmarks involving large deformations are considered. The first benchmark describes a thin vibrating bar, where the displacement is caused by an initial velocity field. For this benchmark, PS-MPM on a structured triangular grid is compared with a reference solution. A second benchmark considers a unit square undergoing axis-aligned displacement with known analytical solution. The spatial convergence of PS-MPM on an unstructured grid is determined for this benchmark. The third benchmark describes a solid column under self-weight. In this benchmark, we investigate the use of a lumped mass matrix in PS-MPM to stabilise the simulation when part of the original domain becomes empty. We will show that a

Vibrating bar
In this section, a thin linear-elastic vibrating bar is considered with both ends fixed. A one-dimensional UL-FEM [36] solution on a very fine grid serves as reference. The grid used for PS-MPM and the initial particle positions are shown in Fig. 9.
The bar is modelled with density ρ = 25 kg/m 3 , Young's modulus E = 50 Pa, Poisson ratio ν = 0, length L = 1 m, width W = 0.05 m and initial maximum velocity of 0.1 m/s. The chosen parameters result in a maximum normal strain in the x-direction of approximately 7%. At the left and right boundary, homogeneous Dirichlet boundary conditions are imposed for both x-and y-displacement, whereas at the top and bottom boundary, a homogeneous Dirichlet boundary condition is imposed only for the y-displacement, and a free-slip boundary condition for the x-direction. The initial displacement equals zero, but an initial x-velocity profile is set to v x (x 0 , y 0 ) = 0.1 sin(x 0 π/L). The initial y-velocity is equal to zero.
The time step size for the simulation is Δt = 5·10 −3 s. The Courant number for two-dimensional problems is defined as C = 2Δt h √ E/ρ, in which h is the typical element length, and √ E/ρ is the characteristic wave speed. Due to the ambiguity of h for a PS-refinement on an unstructured triangular grid, we estimate the average edge length of the sub-triangles by h ≈ 0.025 m. In this case, the Courant number is C ≈ 0.56 < 1, satisfying the CFL condition [4]. Figure 10 shows the displacement in the middle of the bar over time and the stress profile through the entire bar at the end of the simulation. Although a relatively coarse PS-MPM grid and few particle are used, the method yields accurate results and a smooth stress profile.
In case of small and large deformations, the energy in the system over time is shown in Fig. 11. Results for small deformations have been obtained by setting the initial maximum velocity equal to 0.001 m/s. In the limit of small deformations, the PS-MPM solution is described exactly by a harmonic oscillator, and the simulation with small deformations is indeed in good agreement with the (sampled) exact solution. For large deformations, however, the solution is no longer perfectly periodic. Nonetheless, the energy over time obtained with PS-MPM strongly resembles the UL-FEM reference solution. Only after 7 s, the PS-MPM simulation

Axis-aligned displacement on an unstructured grid
In this section, we consider two-dimensional axis-aligned displacement on the unit square (L × W , with L = W = 1 m) to investigate the grid-crossing error and spatial convergence of PS-MPM on an unstructured grid. An analytical solution for this problem is constructed using the method of manufactured solutions: the analytical solution is assumed a priori, and the corresponding body force is calculated accordingly. This benchmark has been adopted from [25] and [38]. The analytical solution for the displacement in terms of the reference configuration is given by u y = u 0 sin 2π y 0 sin E/ρ 0 π t + π .
Here, ρ 0 = 10 3 kg/m 3 , u 0 = 0.05 m and Young's modulus E = 10 7 Pa. The corresponding body forces [25] are in which the Lamé constant λ, the shear modulus μ and the normal components of the deformation gradient D x x and D yy are defined as D yy = 1 + 2u 0 π cos(2π y 0 ) sin E/ρ 0 π t + π . (46) Here, ν = 0.3 and all solutions are again given with respect to the reference configuration. This benchmark was simulated with standard MPM and PS-MPM, using an unstructured triangular grid with material points initialised uniformly over the domain, as shown in Fig. 12. A time step size of Δt = 2.25 · 10 −4 s was chosen, which results in a Courant number of C ≈ 0.72. For the adopted parameters, the imposed solution in Eqs. 40-41 has period T = 0.02 s.

Grid-crossing error
First, it is shown that the grid-crossing error typical for standard MPM does not occur in PS-MPM, by comparing the normal horizontal stress resulting from these methods (see Fig. 13). The configurations for PS-MPM and standard MPM contain the same number of particles (5120, 4 times as many as in Fig. 12) and a comparable number of basis functions, 289 for standard MPM and 243 for PS-MPM. As expected, the stress field in standard MPM suffers severely from grid-crossing, whereas with PS-MPM a smooth stress field is observed. The same conclusion can also be drawn when investigating individual particle stresses over time, as is shown in Fig. 14. For this figure, the particle positioned at (x 0 , y 0 ) ≈ (0.25 m,0.47 m) has been traced throughout the simulation. The figure depicts the stress and trajectory of the particle. For standard MPM, it is observed that the stress profile is highly oscillatory, resulting in a disturbed trajectory. PS-MPM shows no oscillations, but only a small offset, which was found to disappear upon further refinement of the grid.

Spatial convergence
PS-MPM with quadratic basis functions is expected to show third-order spatial convergence. To determine the spatial convergence, the time averaged root-mean-squared (RMS) error is used: where n t denotes the total number of time steps of the simulation. Under the assumption that the spatial error is much larger than the error produced by the numerical integration and time-stepping scheme, the spatial convergence of standard MPM and PS-MPM is determined by varying the typical element length h. This length was defined as the average edge length for standard MPM and the average sub-triangle edge length for PS-MPM. It has been observed that the timeintegration error is indeed sufficiently small, but the number of particles required for an adequately accurate integration increases rapidly as h decreases. Figure 15 shows the spatial convergence of both standard MPM and PS-MPM for different numbers of particles per element on an unstructured grid. Provided that a sufficient number of particles is adopted, standard MPM shows secondorder spatial convergence, and PS-MPM converges with third order. Both orders of convergence correspond to the order of the projection-error, as shown in Sect. 3.5. Besides higherorder spatial convergence, PS-MPM also results in a smaller RMS-error compared to standard MPM, even for the coarsest grid.
To achieve the optimal convergence order for both standard MPM and PS-MPM, the number of integration particles required increases rapidly, as the integration error otherwise dominates the total error, as shown in Fig. 15. Note that the convergence rate in the number of particles for both standard MPM and PS-MPM is measured to be of first order. However, for a fixed typical element length h and number of particles per element, the use of PS-MPM leads to a lower RMS-error, illustrating the higher accuracy per degree of freedom. The inaccurate integration in MPM due to the use of particles as integration points is known to limit the spatial convergence. Different measures have been proposed to decrease the quadrature error in MPM based on function reconstruction techniques like MLS [13], cubic splines [37] and Taylor least squares [39]. The combination of PS-MPM with function reconstruction techniques to obtain optimal convergence rates with a moderate number of particles is subject to future research to further improve PS-MPM for practical applications.

Column under self-weight
In the previous examples, we considered PS-MPM with a consistent mass matrix. However, lumping of the mass matrix is common practice in many applications of MPM, as it speeds up the simulation and increases the numerical stability. In case elements become (almost) empty, the consistent mass matrix in Eq. 13 has very small entries and becomes ill-conditioned, leading to stability issues, which is a wellknown phenomena in MPM [32]. Using a lumped version of the mass matrix to solve for the acceleration and velocity fields is known to overcome the ill-conditioning [32]. However, as explained in Sect. 3.5, lumping within PS-splines can cause spatial oscillations. We show that this is indeed the case and propose an alternative based to lumping to overcome the ill-conditioning of the mass matrix while mitigating spatial oscillations as much as possible.
First, we demonstrate that the standard lumping technique is poorly suited for PS-MPM, by considering a solid col- Fig. 12 The exact solution in which particles (marked with dots) move back and forth along the marked vectors (left) and the unstructured grid with the initial particle configuration (right) Fig. 13 The interpolated particle stress in the x-direction at t = 0.016s  umn under self-weight, as shown in Fig. 16. The column in this exemplary benchmark is modelled as a linear-elastic material. It is fixed in all directions at the bottom, and completely free at the top. At the left and right boundary, we impose a free-slip boundary condition for movement in the y-direction, but the displacement in the x-direction is fixed to be zero. The column is compressed by gravity due to its self-weight.
The solid column is modelled with density ρ = 1 · 10 3 kg/m 3 , Young's modulus E = 1 · 10 5 Pa, Poisson ratio ν = 0, gravitational acceleration g = −9.81 m/s 2 , height H = 1 m and width W = 0.1 m. The maximum strain when adopting these parameters is approximately 18%. Figure 16 illustrates the grid considered for this benchmark as well as the particle positions at maximal deformation.
Due to the fact that part of the grid becomes empty, the use of a consistent mass matrix leads to stability issues with unlumped PS-MPM. The use of a lumped mass matrix within PS-MPM restores the stability, but causes spatial oscillations in the stress profile, as was also discussed in Sect. 3.5. The oscillations degrade the solution of the displacement and velocity as well, as shown in Fig. 17 (left column). Since these spatial oscillations are typical for lumped PS-MPM, common remedies to solve problems due to lumping the mass matrix, like the momentum formulation method [33] or distribution coefficient method [17], do not reduce the spatial oscillations in lumped PS-MPM. Instead, we propose the use of partial lumping.

Partial lumping to mitigate spatial oscillations
As full lumping within PS-MPM causes spatial oscillations, we will instead consider a partial-lumping approach to overcome the stability issues due to almost empty elements, while at the same time minimising spatial oscillations due to lumping. Partial lumping only lumps those rows responsible for the ill-conditioning. These rows typically correspond to the basis functions in the part of the grid where very few particles are left. By only lumping these few rows, the simulation is stabilised, while introducing a much smaller lumping error compared to full lumping. Note that the choice of rows to lump is crucial to ensure stability. Lumping more rows increases the robustness of the simulation, but decreases the accuracy. The trade-off between stability and accuracy within partially lumped PS-MPM should be reconsidered for each benchmark.
For the column under self-weight, we implemented partial lumping by lumping all rows corresponding to a basis function with at least one empty main-triangle in its molecule. Additional empty elements are added at the top of the column, to ensure that the top-most basis functions are always lumped. Figure 17 (right column) shows the displacement Fig. 18 The equilibrium stress at individual particles in a damped solid column. The numerical solutions were determined with fully lumped and partially lumped PS-MPM at t = 10 s, at which the solution was in equilibrium and velocity of a particle over time, as well as the stress profile at t = 2.5 s obtained with partially lumped PS-MPM. Compared to the results obtained with fully lumped PS-MPM, shown in Fig. 17 (left), the stress profile, the velocity and displacement over time significantly improve. Furthermore, the total energy fluctuates less for partially lumped PS-MPM than for fully lumped PS-MPM and also shows less numerical dissipation. However, for very long simulations, partially lumped PS-MPM requires smaller time steps than fully lumped PS-MPM in order to remain stable.
Finally, we compared the equilibrium solution for the stress profile in the solid column with partial and full lumping. The analytical equilibrium solution is given by [41]: σ yy (y eq ) = σ yy (0) (1 − y eq /H ) + κ(y eq /H ) where y eq refers to the equilibrium position of the particles. Furthermore, σ yy (0) = H ρg and κ = σ yy (0) 2E . To obtain the static solution based on the dynamic calculations, a damping term has been added: where, α = 0.6. The damping force damps the PS-MPM solution at each time step in the opposite direction of the velocity field. The equilibrium solution has then been determined by simulating the system with both fully lumped and partially lumped PS-MPM until a stable solution was reached. Both the analytical and numerical static solutions are shown in Fig. 18. For fully lumped PS-MPM, the spatial oscillations in the stress profile remain visible in the equilibrium solution. However, the stress profile of partially lumped PS-MPM is in close agreement with the exact solution.
Hence, the use of partial lumping within PS-MPM combines the advantages of adopting a consistent and lumped mass matrix, while minimising their drawbacks, and thereby improves the overall performance of the method. Future research will focus on techniques to further mitigate or fully overcome the oscillations caused by lumping, in particular when PS-MPM is applied for practical problems.

Conclusion
In this paper, we presented an alternative for B-spline MPM suited for unstructured triangulations using piecewise quadratic C 1 -continuous Powell-Sabin spline basis functions. The method combines the benefits of smooth, higher-order basis functions with the geometric flexibility of an unstructured triangular grid. PS-MPM yields a mathematically sound approach to eliminate grid-crossing errors, due to the smooth gradients of the basis functions. As a first validation, a vibrating bar was considered, for which the PS-MPM solution yields accurate results on a relatively coarse grid. Numerical simulations obtained for a unit square undergoing axis-aligned displacement have shown higher-order convergence for the particle stresses and displacements.
The use of a lumped mass matrix to increase stability within PS-MPM leads to spatial oscillations in the stress profile. A partial lumping strategy was proposed to combine the advantages of adopting a consistent and lumped mass matrix, which successfully mitigates this issue. Investigation of alternative unstructured spline technologies, in particular, refinable C 1 splines on irregular quadrilateral grids [21] is underway.

Compliance with ethical standards
Conflict of interest Pascal de Koster has received funding support from Dutch research institute Deltares.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.