Comparison and unification of material-point and optimal transportation meshfree methods

Both the material-point method (MPM) and optimal transportation meshfree (OTM) method have been developed to efficiently solve partial differential equations that are based on the conservation laws from continuum mechanics. However, the methods are derived in a different fashion and have been studied independently of one another. In this paper, we provide a direct step-by-step comparison of the MPM and OTM algorithms. Based on this comparison, we derive the conditions, under which the two approaches can be related to each other, thereby bridging the gap between the MPM and OTM communities. In addition, we introduce a novel unified approach that combines the design principles from B-spline MPM and the OTM method. The proposed approach does not contain user-defined parameters and can decrease the costs of the standard OTM method. Moreover, it allows for the use of a consistent mass matrix without stability issues that are typically encountered in MPM computations.


Introduction
The material-point method (MPM) [48,51] is a numerical technique suited to model large deformations in continuum mechanics. MPM has been successfully applied in the numerical simulation of complex engineering problems [3,50,59,61]. It originates from the fluid-dynamics-oriented particlein-cell (PIC) method [25] and is mainly based on a weak formulation of the momentum-balance equation. To avoid mesh distortion when the material undergoes large deformations, MPM combines a fixed Eulerian background grid and a set of Lagrangian point masses called material points. The MPM grid contains no permanent information about the continuum and is used only to assemble and solve the discretized governing equations. The material points travel freely through the background grid, while carrying the information about the continuum. The optimal transportation meshfree (OTM) method [27] has been developed to simulate general solid and fluid flows and applied to a wide range of problems [19,28,29,37]. The OTM method uses the concepts from optimal transportation theory (an overview is provided by Villani [54]) to translate the mass and momentum balance equations into a minimization problem for the total action of the solid over a time interval. The OTM method employs two sets of points: nodal points and material points. Nodal points carry information about the positions, while material points represent the continuum. This updated Lagrangian method is typically used with local maximum-entropy (maxent) basis functions [4,46] that are fully defined by the nodal set and the domain of analysis.
Many similarities can be found between MPM and the OTM method. For instance, they both employ the idea of material points that represent the continuum, but are not used directly to compute the solution of the governing equations. Moreover, an alternative derivation of the OTM scheme has been provided by Weißenfels and Wriggers [56], where the method is obtained from the weak form of the equation of motion. Despite this, MPM method and the OTM method have been evolved and studied independently from each other. An in-depth analysis and direct comparison of the methods provide a better understanding of their relation, with potential improvements of MPM based on the present knowledge of the OTM method and vice versa.
This paper consists of two parts: The first part offers new insights into the relation between MPM and the OTM method by drawing a detailed comparison of their algorithms. Based on this comparison, it identifies the conditions, under which the two approaches can be related to each other, and highlights their fundamental differences. The second part of the study presents a novel unified approach that combines the design principles from B-spline MPM (BSMPM) [43] and the OTM method. The choice of BSMPM is motivated by previous studies. First of all, a number of studies have indicated that BSMPM is a viable alternative not only for MPM, but also for its more advanced versions such as dual domain material-point (DDMP) [60] and convected particledomain interpolation (CPDI) [40] methods [20,45,53,57]. Secondly, Cyron et al. [13] pointed out the strong similarities between B-spline and maxent basis functions. For example, both of them are non-negative, smooth and possess the partition of unity property. Although in the present paper, the unified approach is applied to relatively simple examples, all derivations are presented in their general form enabling its straightforward extension to more complex problems.
This paper is structured as follows: Section 2 presents the governing equations that serve as a basis for the development of both MPM and OTM. To be self-contained and allow for a direct comparison of the methods, Sects. 3 and 4 describe, respectively, the MPM and OTM schemes and provide their computational algorithms. Section 5 gives a step-by-step comparison of the algorithms and identifies the conditions, under which they can be related, as well as their main differences. Aside from the computational algorithms, the MPM and OTM methods can be classified by means of the basis functions typically used within them. For this reason, Sect. 6 describes the piecewise-linear Lagrange and maximum entropy basis functions. Since our unified approach and BSMPM adopt the B-spline basis functions, the section also presents the B-spline basis functions. After that, Sect. 7 introduces the unified approach that combines the advantages of BSMPM and the OTM method. Section 8 reports the results that show the difference between the MPM and OTM methods numerically. In addition, it demonstrates the results obtained with our unified approach and gives a comparison with BSMPM. Finally, Sect. 9 summarizes the main conclusions of this study.

Governing equations
To describe the governing equations for the MPM and OTM schemes, we start by introducing some notation from contin-uum mechanics. It is assumed that the considered one-phase continuum occupies the domain Ω 0 ⊆ R 3 at the initial time t 0 and domain Ω t ⊆ R 3 at any later time t > t 0 . The initial position of the material is denoted by In abstract formulation, the deformation mapping is defined as where T is the final time. For a fixed time t, the mapping Ω 0 → Ω t can be considered as the 'push forward' operator, which needs to be bijective, so that the current domain can also be 'pulled back' to the initial one via ϕ −1 (x, t) : Ω t → Ω 0 . Likewise, it links the initial and current positions: The displacement, velocity, and acceleration vectors are denoted by respectively. The displacement at time t is defined as the difference between the current and initial positions: The velocity and acceleration are defined by means of the material time derivative.
T , the material time derivative can be written as Since the convective effects can be neglected in Lagrangian computations (e.g., Donea et al. [16]), the velocity and acceleration are obtained from Furthermore, the deformation gradient tensor F is defined as where I is the identity matrix. For a continuum under purely mechanical loading, the motion is governed by the conservation laws given by the following system of coupled partial differential equations (PDEs): where ρ is the mass density, σ is the Cauchy stress tensor, and b is a body force. Equations (8) and (9) describe the mass and momentum balance in the updated Lagrangian form, respectively. To specify the material, the above system needs to be supplemented with a suitable constitutive equation. In this paper, only linear elastic and neo-Hookean (hyperelastic) constitutive relations are considered. A multidimensional linear elastic model can be described as where λ and μ are Lamé's first parameter and the shear modulus, respectively. Denoting the determinant of F by J , the neo-Hookean material model that can be used to describe a nonlinear stress-strain material behavior is given by The details on the derivation of the conservation laws for a solid continuum and information on various constitutive models can be found in the works of Malvern [32] and Spencer [42]. It should be noted that while the choice of the constitutive model should not affect the approaches, it may lead to more complex algorithms for more elaborate nonlinear stress-strain relations.
3 Material-point method MPM originates from the particle-in-cell (PIC) [25] and fluid-implicit particle (FLIP) [10] methods. After its introduction, several related methods have been proposed. For example, Bardenhagen and Kober [6] introduced a generalization of MPM, the generalized interpolation material-point (GIMP) method, where particles are represented by particle characteristic functions. In addition to the related methods, several versions of the MPM algorithm have been suggested. First, the original update stress last (USL) [48] version of MPM was transformed into the modified update stress last (MUSL) [51]. The main difference between the USL and MUSL algorithms is that USL computes the nodal velocities directly from the nodal accelerations, while MUSL obtains them from the material-point velocities, significantly improving the stability of the scheme. Next to USL and MUSL, the update stress first (USF) [5] and update stress averaged (USAVG) [36] MPM algorithms have been introduced. Both USF and USAVG attempt to improve the energy conservation within USL and MUSL. However, for the comparison of MPM and OTM algorithms, we focus on the most basic version of the MPM scheme, USL. Finally, a version of MPM, the so-called moving-mesh MPM, has been proposed where the background grid does not remain fixed as the simulation progresses in time. Moving-mesh MPM has been successfully applied to model complex problems such as the biological mechanics of cells [21] and the texture evolution in polycrystalline nickel [62]. MPM simplifies the momentum equation, Eq. (9), by assuming that the mass density is constant in time and solves the simplified equation on the fixed background grid adopting a variational formulation. The weak form of the momentum equation is given by where ω denotes an element of the test space that consists of all functions, which are sufficiently smooth and zero on the part of the boundary where essential boundary conditions are prescribed, Γ t = ∂Ω t is the boundary of the domain Ω t , and n is the outward unit normal vector. By adopting the standard approach in the finite element method (FEM) [11], Eq. (12) can be written in the following matrix-vector form for each direction x k : whereā k = [a k,1 , a k,2 , . . . , a k,N ] T is the unidirectional vector of the unknown acceleration coefficients, 1 M = M i j ∈ R N ×N is the consistent mass matrix, andf k = [ f k,1 , f k,2 , . . . , f k,N ] T is the unidirectional force vector. Here, N being the total number of nodes on the background grid. The entries of the mass matrix and force vector are given by respectively. Here, τ (x, t) is the prescribed traction at the boundary. In general, MPM can be implemented combining Eq. (13) for each k into one linear system. The implementation procedure of a multi-dimensional MPM can be found, for instance, in the thesis of Kafaji [2]. In MPM computations, the consistent mass matrix is typically replaced by the row-sum lumped mass matrix M L , which can be obtained by summing the off-diagonal entries of M in each row, adding them to the diagonal entry, and subsequently setting the off-diagonal entries to zero. If the basis functions maintain a partition of unity within the domain, that is the diagonal entry of the lumped mass matrix can be expressed as follows: Furthermore, in MPM material points represent the material and carry all physical information about it (e.g., the mass, strain, and stress). While most material-point properties vary in time, the mass is time independent. This assures that the conservation of mass, Eq. (8), is satisfied. Furthermore, throughout a simulation, the integrals in Eqs. (14) and (15) are approximated by projecting the material-point information onto the background grid. Let the continuum be discretized by M material points, then the integral of an arbitrary vectorvalued function g(x) is approximated by where V p and x p represent the material-point volume and position, respectively. The information obtained by solving Eq. (13) is then mapped back to update the material points. Compute mass matrix M s , linear momentum vectorq s k , and force vectorf s k : Compute nodal accelerations: Compute nodal velocities and incremental nodal displacements: Update material-point coordinates:

9
Update material-point deformation gradients: 10 Update material-point volumes: Update material-point densities: In the considered version of MPM, the Euler-Cromer time-stepping scheme [12] is used for temporal discretization. The scheme applies the forward Euler method to advance the velocity and employs the backward Euler method for the displacement. The choice of the time integration technique is discussed in detail by Wallstedt and Guilkey [55]. A schematic representation of the MPM solution strategy that is followed within each of the N time steps is provided in Fig. 1. The complete form of the original MPM algorithm is presented in Algorithm 1, where it is assumed that only a constant body force is acting externally. In the algorithm, the notation used for the acceleration and force vectors,ā k andf k , is also adopted for global nodal vectors (e.g., the nodal velocity vector in the direction

Optimal transportation meshfree method
The OTM method is a meshless updated Lagrangian method that is based on the concepts from optimal transportation theory. In contrast to MPM, the OTM method explicitly includes the dependence of the mass density on time. For an arbitrary time interval t 0 , T , it assumes that the density at time t 0 and T is prescribed: where x 0 and x are given in Eq. (1). Benamou and Brenier [8] note that the mass and momentum balance equations together with Eqs. (32) and (33) can be translated into a minimization problem for the action of the solid over the time interval. The action over t 0 , T is given by where K (ρ, v) is the kinetic energy, which is equal to Although Eqs. (34) and (35) are expressed in terms of mass density and velocity, in the OTM framework the flow is described by means of the deformation mapping ϕ that is defined in Eq. (1). The deformation mapping is related to velocity and density in the following way [27]: Benamou and Brenier [8] also demonstrate that the deformation mapping that minimizes the action in Eq. (34) is given in terms of McCann's displacement interpolation [33]: Here, ϕ x 0 , T is the optimal transportation map of ρ 0 into ρ T with respect to the cost function [27]: in which γ x 0 , t is a generic mapping of mass density.
To generate a numerical scheme, Eq. (34) is discretized in space and time. The time interval t 0 , T is divided into sub-intervals t s , t s+1 with s = 0, 1, . . . , N − 1, to which the above theory developed for t 0 , T is still applicable. The OTM method approximates a flow using the concept of the free energy of the solid U [29]: where f is the local free-energy density per unit volume. Furthermore, the method employs the Wasserstein distance d W between mass densities at two consecutive time instances [27,28]: For elastic materials in unforced systems, the semidiscrete action sum is equal to [28] where ϕ s is the deformation mapping at time t s . For the spatial discretization of Eq. (42), the OTM method employs two sets of points: nodal points and material points (see Fig. 2). Nodal points carry position information, while material points represent the continuum body and arise from the spatial approximation of the mass densities by point masses: where δ x − x s p is the Dirac delta distribution centered at x s p . Material points are convected by the deformation: in which ϕ s→s+1 h (x) is the incremental deformation map. Fedeli et al. [19] explain that it can be described by general linear interpolation schemes of the form: whereby the basis functions are assumed to be consistent. Consistent basis functions satisfy the following conditions [19]: -Linear completeness Moreover, material points carry a fixed mass, serve as integration points, and store all local state data. Compute mass matrix M s , linear momentum vectorq s k , and force vectorf s k : Update nodal coordinates: Update material-point coordinates: Update material-point volumes: Update material-point densities: Update material-point deformation gradients: Set s = s + 1 13 end Algorithm 2: Optimal transportation meshfree method.
Using the definition of the Wasserstein distance and free energy, as well as Eq. (43), the fully discrete form of Eq. (42) is obtained (see Appendix A for more details): The OTM algorithm originates from applying discrete Hamilton's principle [24] to the fully discrete action [22,27]. A basic OTM algorithm for a solid material is summarized in Algorithm 2. For more details on the OTM method, we refer to the work of Li et al. [27] and Habbal [22].

Comparison of algorithms
While the derivations of MPM and the OTM method are fundamentally different, the resulting algorithms have many similarities. This section provides a side-by-side comparison of the computational steps from Algorithm 1 (USL-MPM) and 3 (OTM). Based on this comparison, it summarizes the conditions required to further relate the methods and highlights the principal differences between them.
-In the beginning of the simulation, both algorithms initialize the nodal coordinates and material-point properties. Assuming that the material-point velocity in MPM can be written as A direct substitution of Eq. (57) into the expression for the linear momentum in the MPM algorithm (Eq. 20) leads to the linear momentum formula used in the OTM method (Eq. 49). -Furthermore, steps 5 and 6 of the USL algorithm are implicitly included in step 5 of the OTM algorithm, where the nodal coordinates at time step s + 1 are computed. More precisely, from step 5 in the OTM algorithm, we obtain At the same time, the incremental nodal displacement in MPM can be written as It can be seen that Eq. (59) is equal to Eq. (58) for constant time-step sizes. From this, we conclude that the update of nodal positions is identical for both methods when the time-step size is fixed. -The definition of the incremental transport map (Eq. 45) implies that in step 6 of the OTM algorithm, material-point positions are obtained from while step 7 in the MPM algorithm states that The equality of the above expressions can be shown if linear completeness (see Eq. 47) of the MPM basis functions is imposed. That is, if the following condition is satisfied: Equation (61) can be rewritten as -Moreover, the OTM scheme avoids a direct update of the material-point velocity by adopting Eq. (57), whereas MPM performs the update in step 8. Nevertheless, assuming that it is possible to relate the methods again. Substituting Eq. (64) in step 8 of MPM gives v s+1 Substituting the updated nodal velocity from step 6 of the MPM algorithm yields v s+1 Therefore, the updated material-point velocity in MPM is equal to The above equalities follow from the computation of the incremental nodal displacements in steps 6, and the update of material-point coordinates in step 7 of the MPM algorithm. The final expression in Eq. (68) is identical to the implicit material-point velocity update in the OTM algorithm.
Although an extra assumption is required to establish a connection between the methods, it does not lead to an essential disparity between them. -Step 9 of both schemes can be shown to be identical as well. From step 6 in the OTM algorithm and the OTM definition of material-point velocity presented in Eq. (57), it follows that Therefore, step 9 is the same for the MPM and OTM algorithm. -Next, the update of the material-point volumes is investigated. MPM performs this update in step 10, while the OTM method computes the volume in step 7. To prove the equivalence of those steps, we need to show that This identity can be proved in the following way: In the above, Eq. 70 and the definition of the deformation gradient tensor are used. -The remaining steps are identical for MPM and OTM.
In this section, we have shown that under certain conditions the MPM and OTM algorithms can be related to each other. Namely, assuming a constant time step and the validity of Eqs. (57), (62), and (64) for the MPM scheme, the only difference between the standard algorithms emerges from the update of the basis functions and their gradients.
However, this difference is fundamental. Since the nodes are fixed to their initial positions in the MPM basis-function update (step 3 of the USL algorithm), the method is considered to be a combination of Lagrangian and Eulerian approaches. At the same time, the OTM method is an updated Lagrangian particle method. For the implementation of the methods, this implies that the OTM method only discretizes the initial material domain, while MPM discretizes the complete domain, where the material is allowed to move, as well as the initial material domain. Consequently, the OTM method does not include inactive elements.

Basis functions
Section 5 identifies strong similarities between the OTM and MPM schemes without specifying the basis functions used within the methods. Since MPM typically uses piecewiselinear basis functions, while the OTM method exploits the local maximum-entropy (maxent) basis functions, both piecewise-linear and maxent basis functions are presented in this section. In addition, this section describes the Bspline basis functions frequently combined with the MPM algorithm and adopted in the unified approach provided in Sect. 7. For simplicity, the piecewise-linear and B-spline basis functions are presented in their one-dimensional form. The extension of these basis functions to quadrilateral elements is obtained by means of the tensor product.

Piecewise-linear basis functions
Given a set of nodes {x i } N i=1 , piecewise-linear basis functions are defined as The corresponding derivatives are equal to These nonnegative basis functions possess the partition of unity property and have a compact support. In addition, they can be implemented in a straightforward manner. However, their gradients are discontinuous at the element boundaries. Within MPM, this can lead to unphysical results when material points cross element boundaries [6]. Figure 3 provides an example of quadratic B-spline basis functions and their gradients. Generally, B-spline basis functions are defined in the parametric space, based on a knot vector. A knot vector in one dimension is a sequence of ordered non-decreasing coordinates. It is typically denoted as Ξ = ξ 1 , ξ 2 , . . . , ξ nb+d+1 , where ξ j ∈ R is the jth knot, nb is the total number of basis functions, andd is the polynomial order. The knot vector is uniform when its knots are distributed equidistantly. The knots are called repeated when more than one knot is positioned at the same coordinate in the parametric space. An open knot vector contains the first and last knotsd + 1 times ensuring that the resulting basis functions are interpolatory at the boundaries of the domain.

B-spline basis functions
A non-empty knot interval ξ j ξ j+1 is referred to as a knot span. For an open uniform knot vector, the number of spans is equal to nb −d.
The Cox-de Boor formula [14] defines B-spline basis functions recursively, starting with piecewise constants (no repeated knots, i.e.,d = 0): Ford > 0, the basis functions are given by whereẑ is the parametric domain. The derivatives of the Bspline basis functions can be computed from [14]: It should be noted that in the considered implementation of BSMPM, the parametric and physical domains are the same. Moreover, first-order B-splines basis functions are identical to piecewise-linear basis functions described previously. B-spline basis functions were introduced within MPM by Steffen et al. [43,44]. They form as a potential remedy to the grid-crossing and quadrature errors present in MPM. The concept of BSMPM was further investigated by combining B-spline basis functions with more advanced integration techniques [20,53,57] and adopting the method to unstructured grids [15].

Maximum-entropy basis functions
Maxent basis functions were introduced by Sukumar [46] for the construction of polygonal interpolants. Arroyo and Ortiz [4] presented local maxent basis functions and first used them within a meshfree method. After that, local maxent basis functions have been integrated within several meshfree schemes such as point collocation methods [18]. The schemes that combined meshfree methods with maxent basis function have been applied to a variety of problems including simulations of shear-deformable plates [23] and thin-shell analysis [34].
The construction of maxent basis functions combines elements from probability theory and optimization. In fact, within a convex hull of the nodal set x i (i.e., the smallest convex set that contains all nodes), the set of local maxent basis functions φ i (x) ≥ 0 N i=1 form the solution of the following constrained optimization problem [38]: where w i (x) is a nonnegative weight function or prior estimate of φ i . The solution of this problem is typically found using the method of Lagrange multipliers [7] and can be written as where λ represents the Lagrange multipliers. In practice, the primal problem of maximization is transferred into the dual problem of minimization [46]. Considering the new formulation, λ is expressed as [1,9] λ = arg min ln(Z ). (83) The Lagrange multipliers are typically found via Newton's method. The solution procedure of this method (e.g., [39]) requires an initial guess for λ and both first-and second-order partial derivatives of ln(Z ) with respect to λ. The first derivatives of the local maxent basis functions are given by (e.g., [58]) in which the matrices H and A are computed in the following way: Here, is the dyadic product (i.e., the dyadic product of any two vectors is equal to a b = ab T ).
In the above description, the weight function remained unspecified due to a large number of viable options (e.g., Gaussian prior, cubic or quartic spline). In general, the prior functions are defined by means of the normalized radius of the support domain, r i , for node i: where || · || is the L 2 -norm, and d i is the size of the domain of support of node i, which is a user-defined parameter. In this paper, d i is equal to in which d max is a factor with a typical value between 2.0 and 4.0 (this value is selected by the user), and ζ i is the distance between node i and its nearest neighboring node. Alternative definitions of the normalized radius of the support domain can be found, for instance, in the work of Sukumar and Wright [47] and Yaw et al. [58]. Cubic spline weight functions, which are used in this paper, are given by (89) Figure 4 shows the one-dimensional cubic spline basis functions and their derivatives, as well as the corresponding local maxent basis functions together with the derivatives. Local maxent basis functions possess many desirable properties for meshfree algorithms. First of all, they are entirely defined by the nodal set and the domain of analysis. They are also non-negative, satisfy the partition of unity property, and provide an exact approximation for affine functions [4]. Furthermore, the local maxent basis functions have the so-called weak Kronecker delta property [4]. However, the local maxent basis functions are defined only within a convex hull of the nodal set. If non-convex domains are considered, the basis functions lose the weak Kronecker delta property at the non-convex parts of the domain [4]. In addition, the calculation of Lagrange multipliers is numerically challenging [22] and frequently requires significant computational costs.

Unified approach
The unified approach combines the OTM scheme and BSMPM. Similar to BSMPM, the proposed approach uses Bspline basis functions. However, it computes them based not only on the updated material-point positions, but also on the advected degrees-of-freedom set. B-spline basis functions of any order possess the partition of unity property, but only first-order B-spline basis functions satisfy the linear completeness property. Therefore, according to the definition in Sect. 4, higher-order basis functions are not consistent. Since Algorithm 2 is designed for consistent basis functions, the unified algorithm is mainly based on Algorithm 1 to ensure the compatibility of the unified algorithm with the higherorder B-spline basis functions. Furthermore, the proposed approach employs the consistent mass matrix. The unified algorithm is presented in Algorithm 3.
Both MPM and OTM provide motivation for the proposed unified method. On the one hand, the addition of the advected nodal points to update the basis functions is supposed to stabilize the computation when BSMPM is combined with a consistent mass matrix. Consistent mass matrices frequently cause stability issues in MPM [31], and BSMPM inherits these issues. For this reason, MPM and BSMPM are generally used with a lumped mass matrix. While mass lumping has little influence on the solution quality of lower order methods, its O(h 2 ) approximation of the consistent mass matrix [44] can significantly influence the spatial convergence of higher-order methods such as BSMPM. Moreover, as was mentioned earlier, previous studies demonstrate that moving-mesh MPM can be used successfully for complex simulations.
On the other hand, within the OTM framework, the use of B-spline basis functions is expected to significantly reduce the computational costs. In contrast to maxent basis functions, B-spline basis functions do not require the adoption of iterative methods and have a purely analytical definition. It has been pointed out by Cyron et al. [13] that maxent and higher-order B-spline basis functions have many common properties. For example, they are both smooth, non-negative, have compact support, and satisfy the partition of unity property. Therefore, the use of B-spline basis functions provides a viable alternative to the massively parallel implementation of OTM (pOTM) [30]. In addition, B-spline basis functions do not require the use of search algorithms, frequently added to the standard OTM scheme for stabilization [22,56]. Finally, the adoption of higher-order B-spline basis functions can lead to higher-order spatial convergence. Compute mass matrix M s , linear momentum vectorq s k , and force vectorf s k : Compute accelerations at degrees of freedom: Compute incremental displacement and updated coordinates for degrees of freedom: Update material-point coordinates: Update material-point deformation gradients: Update material-point volumes: Update material-point densities:
We remark that in contrast to the OTM scheme, the proposed unified approach cannot be viewed as a meshless method. However, the study by De Koster et al. [15] shows the potential of BSMPM on arbitrary grids. The extension of the BSMPM to arbitrary grids in combination with an efficient remeshing technique (e.g., the remeshing strategy for large deformations proposed by Erhart et al. [17]) might bring the unified approach closer to the meshless algorithms. In this paper, the examples are restricted to relatively simple problems to study the basic properties of the considered methods. Thus, further research is required to evaluate the effect of mesh distortion on the proposed unified approach.

Numerical results
In this section, three benchmarks are considered to illustrate the performance of the discussed methods. The onedimensional benchmarks describe the vibration of a bar, but have fundamentally different motion triggers and boundary conditions. In the first benchmark, where both ends of the bar are fixed, the domain contains only filled elements allowing for a straightforward implementation and analysis. The second benchmark, where a traction force is acting at one of the boundaries, contains multiple empty cells throughout the simulation, thereby serving as a representative example for the stability analysis. The last benchmark is two-dimensional; it further extends the numerical analysis of the considered algorithms.
The results are provided for the USL version of the MPM scheme (Algorithm 1) and OTM algorithms, as well as the proposed unified approach. For MPM, piecewiselinear, second-order B-spline and maxent basis functions are employed. The OTM algorithm is used only with consistent basis functions (i.e., piecewise-linear and maxent basis functions). For the one-dimensional benchmarks, the factor d max is set to 2.0 to compute the maxent basis functions, while for the two-dimensional problem, its value depends on the considered algorithm.

Bar with fixed ends
This example describes the vibration of a linear elastic bar of length l with fixed ends. The motion is generated by an initial velocity prescribed along the bar: Here, the length of the bar is set to 1 m, Young's modulus is set to 4 kPa, while the initial density and amplitude of the velocity v 0 are equal to 1 kg/m 3 and 0.6 m/s, respectively. The total simulation time is set to 0.001 s, while the time-step size is equal to 10 −5 s. This relatively small time step is required to minimize the contribution of the temporal error to the total one. Moreover, the number of nodes varies between 8 and 512, while the number of material points per element remains equal to 12. The analytical solution in terms of the displacement, velocity, and stress can be found in the work of Wobbes et al. [57]. For the convergence analysis, the root-mean-square (RMS) error in the displacement is computed. RMS error is defined as follows:  Table 1 provides the results in terms of computational time of all considered methods for 32 and 512 nodes. We remark that the maxent computations were performed without a search algorithm. For the coarse discretization, the MPM scheme-based approaches outperform the ones adopting the OTM algorithm. However, when 512 nodes are used, the efficiency of the computation depends on the choice of the basis functions instead of the algorithm. More precisely, replacing maxent by B-spline basis functions reduces the computational time by at least a factor of 2.
In addition, Fig. 6 illustrates the spatial convergence behavior of the considered methods at the end of the simulation. When piecewise-linear basis functions are used, both MPM and the OTM method demonstrate second-order convergence for relatively coarse grids. However, for fine meshes, the methods behave differently. In fact, MPM suffers from grid-crossing errors that result in a loss of the convergence. The OTM-P1 method preserves the secondorder convergence until the final refinement, where only For MPM and the OTM method, the number of nodes is equal to 512, while for the unified approach, only 32 nodes were used  [44,49,55]. Furthermore, Fig. 6 shows that maxent basis functions lead to significantly lower errors than the piecewise-linear basis functions for both MPM and OTM schemes. With maxent basis functions, the convergence of MPM and the OTM method varies between linear and quadratic. It should be noted that the accuracy of MPM and the OTM method with maxent basis functions can be further improved by adapting more advanced implementations [4,35,52,63]. The use of B-spline basis functions leads to similar results for BSMPM and the unified approach. These methods have third-order convergence until the limiting value is reached.

Bar with dynamic traction boundary conditions
This benchmark describes the motion of a neo-Hookean bar with one free end. The bar is fixed at x 0 = 0 and subjected to a traction force at the free end x 0 = l. The forcing function is equal to τ x 0 , t = δ x 0 − l σ x 0 , t . Defining ω = π/l, the stress is given by The initial length of the bar is set to 1 m, the density is equal to 100 kg/m 3 , and Young's modulus is equal to 100 Pa. The length of the computational domain is set to 1.25 m. A more detailed description that includes an analytical solution for displacement is provided by Steffen et al. [43]. To illustrate the stress profile obtained with different methods, the material domain is discretized by 68 nodes, which results in 85 nodes for the complete domain. The material-domain discretization is sufficient for the unified and OTM methods due to their updated Lagrangian nature, whereas MPM requires the discretization of the complete domain. Each active ele- Fig. 7 Comparison of the considered methods to the analytical solution. The material domain is discretized by 68 nodes ment initially contains 4 particles. The computational time is set to 0.4 s, and the time-step size is equal to 10 −4 s. Figure 7 depicts the obtained results. It shows that maxent and B-spline basis functions eliminate the grid-crossing error in MPM. However, MPM-maxent and BSMPM do not follow the analytical solution at the right edge of the bar. Within BSMPM, these inaccuracies can be significantly reduced by increasing the initial number of particles per elements. This suggests that the errors are caused by insufficient accuracy of the numerical integration in MPM. Thus, advanced numerical integration techniques (e.g., Taylor least squares [57]) may improve the BSMPM solution at the boundary. The inaccuracies within MPM-maxent have a different origin. They are most probably caused by the incomplete set of maxent basis functions, which arise from the presence of inactive elements throughout an MPM simulation. Figure 7 also shows that the OTM-P1, OTM-maxent, and unified methods provide significantly more accurate solutions than their MPM equivalents.
Similar to the benchmark discussed in Sect. 8.1, the use of B-spline basis functions instead of maxent basis functions considerably decreases the computational time for both MPM and OTM methods. In fact, the unified approach and BSMPM computations are approximately 10 times faster than the OTM and MPM computations with maxent basis functions (without a search algorithm) with the settings used for Fig. 7.
Furthermore, the unified approach and BSMPM have the lowest RMS error and highest convergence rates compared to the other methods. This is illustrated in Fig. 8. To minimize quadrature and time integration errors, this figure is obtained placing 12 particles per cell at the beginning of the simulation and reducing the computational time to 0.1 s. In general, the Fig. 8 Convergence behavior of the considered methods. The number of nodes, in this case, refers to the number of nodes used to discretize the material domain obtained convergence orders of the considered algorithms are slightly lower than expected. This can be related to the discontinuities in the solution for the stress field.
The main advantage of the unified method over BSMPM arises from its stability properties. When a material point enters an empty element, BSMPM inherits stability issues from MPM. For instance, changing the total number of nodes to 81 results in a termination after 0.3 s. This issue in MPM is discussed in detail by, for example, Kafaji [2] and requires the use of the MUSL algorithm with a lumped mass matrix to circumvent the breakdown.

Plate undergoing axis-aligned displacement
The final benchmark describes a two-dimensional neo-Hookean plate that is fixed at the entire boundary. The plate is assumed to be a unit square (l × l with l = 1 m), and its motion is triggered by a body force. The corresponding displacement is given by Here, A denotes the maximum amplitude of the displacement, which is set to 0.005 m, whereas c is a constant defined as with Young's modulus E being equal to 1 · 10 7 Pa and initial mass density ρ 0 being set to 1·10 7 Pa. Furthermore, Poisson's ratio is equal to 0.3. A detailed description of this benchmark is available [40] including an analytical solution obtained with the method of manufactured solutions [26,55]. The domain is discretized by 33 nodes in each direction, and each element contains initially 16 particles. The computational time is set to 3.5·10 −3 s, and the time-step size equals 10 −4 s. The tensor product of the one-dimensional basis functions described in Sect. 6 is adopted for discretization. For the unified approach, the knot vector in each direction is updated based on the velocity field to approximate the moving nodes. For large deformations, however, the use of unstructured Bsplines is preferred, such as PS-splines [41], as it gives the possibility to update the inner nodes independently of each other. Figure 9 shows the normal stress in x 1 -direction σ 11 along the plate for material points with x 0 2 ≈ 0.07 m. Employing piecewise-linear basis functions within both MPM and the OTM method leads only to a piecewise-constant approximation of the stress field. Due to grid-crossing errors, the stress profile obtained using MPM with piecewise-linear basis functions deviates significantly from the analytical solution. The simulations performed with the maxent basis functions show considerably more accurate stress approximations for both methods. To obtain these results with the MPM algorithm, the user-defined factor d max is set to 3.0, while for the OTM method, d max = 2.0 is taken. The unified approach and BSMPM lead to an even smoother stress profiles which are in close agreement with the analytical solution.

Conclusions
The first part of this paper provides a comparison between the MPM and OTM schemes. While the methods were derived in fundamentally different manners, the resulting algorithms are closely related. In fact, assuming a constant time step, the validity of the backward Euler scheme for material-point displacement in MPM, as well as the linear completeness of the MPM basis functions and their ability to translate nodal velocities into material-point velocities, the only difference between the algorithms emerges from the update of the basis functions. However, this difference is fundamental. Since MPM uses initial nodal positions in the basis-function update, it is viewed as a combination of Lagrangian and Eulerian approaches. At the same time, the OTM method is a fully updated Lagrangian method. Moreover, MPM is typically used with piecewise-linear basis functions, whereas the OTM method generally employs maximum-entropy basis functions.
In the second part of the paper, a unified approach is proposed. This approach combines BSMPM and the OTM method. Similar to BSMPM, the proposed approach uses B-spline basis functions. However, it computes the basis functions based not only on the updated material-point positions, but also on the advected degrees-of-freedom set. The obtained numerical results demonstrate that the proposed method preserves the convergence properties of BSMPM and remains stable when a consistent mass matrix is adopted. Furthermore, the unified approach does not contain user-defined parameters and is faster than computations with maxent basis functions for large size problems.