Conservative Taylor least squares reconstruction with application to material point methods

Within the standard material point method (MPM), the spatial errors are partially caused by the direct mapping of material‐point data to the background grid. In order to reduce these errors, we introduced a novel technique that combines the least squares method with the Taylor basis functions, called the Taylor least squares (TLS), to reconstruct functions from scattered data while preserving their integrals. The TLS technique locally approximates quantities of interest such as stress and density, and when used with a suitable quadrature rule, it conserves the total mass and linear momentum after transferring the material‐point information to the grid. The integration of the technique into MPM, dual domain MPM, and B‐spline MPM significantly improves the results of these methods. For the considered examples, the TLS function reconstruction technique resembles the approximation properties of highly accurate spline reconstruction while preserving the physical properties of the standard algorithm.


INTRODUCTION
The material point method (MPM) 1,2 is a continuum-based numerical tool to simulate problems that involve large deformations. In MPM, the material is represented by a set of Lagrangian particles (material points) that move through a fixed Eulerian background grid. The material points carry the physical properties of the continuum such as the mass, strain, and stress. At the beginning of each time step, this information is projected from the particles to the degrees of freedom (DOFs) of the background grid, where, similarly to the finite element method 3 (FEM), the discretized governing equations are assembled and subsequently solved. The obtained information are then mapped back to update the material points. This conjunction of Lagrangian and Eulerian approaches makes MPM well suited for challenging problems involving large strains. Despite its impressive performance for many engineering problems, 4-6 standard MPM still suffers from several shortcomings. For example, when large deformations are considered, the particular case of Shepard interpolation 7 used by the method to project the scattered material-point data to the DOFs can introduce significant numerical inaccuracies. 8,9 In addition, when material points travel from one cell to another, they generate unphysical oscillations in the forces, frequently referred to as grid-crossing error. 10 This is due to the use of piecewise-linear basis functions, whose gradients are discontinuous on element boundaries. Finally, MPM contains FEM-type errors originating from mass-lumping and interpolation, as well as time-stepping errors. 11 Much research has been conducted to overcome these shortcomings. On the one hand, methods such as the generalized interpolation material point (GIMP) method, 10 the convected particle domain interpolation (CPDI) method, 12 the dual domain material point method (DDMPM), 13 and the B-spline material point method (BSMPM) 14,15 have been designed to overcome the grid-crossing error. On the other hand, Sulsky and Gong 8,9 have shown that the shortcomings resulting from the mapping of particle information can be decreased by reconstructing functions with a higher-order technique and evaluating these functions at the desired positions. Furthermore, Tielen et al 15 have pointed out that the use of global cubic-spline function reconstruction technique substantially improves the performance of BSMPM. Although modified mapping generally improves the accuracy of the solution, the standard reconstruction techniques such as spline interpolation might lead to the loss of physical properties of the material point methods. In fact, while MPM projection conserves the mass and linear momentum of the system, most standard mapping techniques do not guarantee that. Therefore, more advanced procedures are required for function reconstruction.
In this paper, we propose a novel Taylor least squares (TLS) reconstruction technique, which is based on the least squares 16 approximation constructed from a set of Taylor basis functions. 17 This technique reconstructs quantities of interest such as stress and density locally, within each active cell, and evaluates them at the integration points. After that, a Gauss quadrature is applied to determine the internal forces and velocities. If a sufficient number of Gauss points is defined within each element, the proposed mapping technique preserves the physical qualities of the standard MPM by conserving the mass and linear momentum of the system.
While the TLS reconstruction can be used autonomously within the MPM algorithm, it can also be combined with methods like CPDI or BSMPM to further reduce the spatial errors. Since previous studies have indicated that DDMPM and BSMPM are viable alternatives not only for MPM but also for GIMP and CPDI, 13,18,19 we have applied the TLS reconstruction technique within MPM, DDMPM, and BSMPM. The methods are tested on two one-dimensional benchmarks describing the deformation of a vibrating bar and column compaction under slowly increasing body force, respectively. The conservation property of the TLS reconstruction technique is verified by computing the total mass and momentum before and after projecting the particle information to the background grid. The accuracy of the material point methods with TLS approximation is investigated either qualitatively or based on the spatial errors and convergence rates. In addition, the obtained results are compared with those computed with the cubic-spline reconstruction technique. 15 We observe that, when material points do not cross element boundaries in the vibrating bar problem, the TLS reconstruction technique has little influence on MPM and DDMPM but significantly improves the convergence of BSMPM. When particles start to cross element boundaries in both examples, the effect of the advanced reconstruction technique becomes clearly evident regardless of the material point method. In fact, the TLS technique smoothens the solutions of MPM, DDMPM, and BSMPM and ensures a close agreement with the reference solutions. Moreover, the TLS approximation outperforms the global spline interpolation in terms of the mass and momentum conservation. Therefore, the TLS reconstruction leads to a conservative projection of the particle data and increases the accuracy of the MPMs. This paper is structured as follows. Section 2 introduces the governing equations capturing the one-dimensional deformation of one-phase continua. Section 3 summarizes the MPM algorithm for numerical discretization and the solution of these equations and outlines the modifications required to obtain DDMPM. Section 4 focuses on the use of B-splines in MPM. Section 5 describes the main concepts of the TLS reconstruction and proves that its use with a sufficiently accurate Gauss quadrature leads to a conservative data projection. Section 6 discusses the numerical results obtained for vibrating bar and column compaction problems. Finally, Section 7 provides the conclusions.

GOVERNING EQUATIONS
In this paper, it is assumed that the considered one-phase continuum occupies a domain Ω 0 ⊆ ℝ at time t 0 and a domain Ω ⊆ ℝ at time t > t 0 . The material does not undergo irreversible deformations. One-dimensional deformations of such a continuum can be described for small strains by a closed coupled system of partial differential equations for velocity and stress in the following manner 20 : The coordinate x refers to the spatial direction, whereas t denotes the time. Moreover, is the density, v is the velocity, is the stress tensor, g is the gravitational acceleration, and E is Young's modulus.
Equation (1) is obtained from the momentum balance equation, whereas Equation (2) describes the linear elastic constitutive relation for small strains. For large deformations, Equation (2) needs to be replaced by 21 This system can be extended by a relation between the velocity v and displacement u given by In order to obtain a unique solution to the above system of equations, initial and boundary conditions have to be prescribed. Therefore, the following initial conditions are required: In addition, two types of boundary conditions are considered.

LOW-ORDER MATERIAL POINT METHODS
The momentum balance equation of a continuum (Equation (1)) is solved by MPM 8,9 in its weak formulation: where v denotes a test function. We introduce vector notation for quantities defined at the DOFs to simplify the description. If the domain discretization generates N n DOFs, we denote the velocity vector by v(t) = [v 1 (t) v 2 (t) … v N n (t)] T . The basis function vector is given by where i represents a basis function associated with the ith DOF. Thus, the velocity is approximated according to the finite-element approach as follows: The displacement u and virtual velocity v are discretized in a similar way.
The semidiscretized weak formulation can then be written as Denoting the acceleration vector a = dv dt , the consistent mass matrix M C = ∫ Ω T dΩ, the external force vector f ext = (x, t) | | Ω − ∫ Ω g dΩ, and the internal force vector f int = ∫ Ω (x, t)∇ dΩ, we obtain It should be noted that, in the considered version of MPM, the consistent mass matrix in Equation (5) is typically replaced by the lumped mass matrix M. The row-sum lumped mass matrix is obtained by summing the off-diagonal entries of each row and adding them to the diagonal entry: Alternatively, for bases satisfying the partition of unity property, the lumping can be performed variationally (ie, M = ∫ Ω dΩ.). Thus, the mass corresponding to the DOFs is stored in the diagonal m of M.
Finally, the stress is updated using the gradient of the basis function vector ∇ (x) and the constitutive relation. For example, for small deformations, Equation (2) yields In MPM, the material points move within the grid and serve as integration points. If the continuum is discretized by N p material points, then for an arbitrary function f, the integral is computed as where V p is the particle volume.
In the considered version of MPM, the Euler-Cromer time-stepping scheme 22 is used for temporal discretization. This method is semi-implicit and applies the forward Euler method to advance the velocity in time and the backward Euler method to advance the displacement. The scheme is energy conservative and conditionally stable.

Algorithm
When the deformation of a one-phase continuum is studied, each particle carries a certain volume V p , density p , position x p , velocity v p , and stress p . All these values are time-dependent, whereas the material-point mass m p is constant throughout the simulation. The superscript t denotes the time level, and Δt is the time step size. At t = 0 s, the material points are initialized (see the work of Al-Kafaji 23 for further details). Assuming that all particle properties are known at time t, the computation for time t + Δt proceeds as follows.
First, the data from the material points is projected to the DOFs. For example, the diagonal of the lumped mass matrix and the internal force vector are computed as The accelerations at the DOFs are then obtained by combining the internal forces with any external forces, namely, Meanwhile, for each DOF, the update-stress-last (USL) scheme 1 calculates the velocities directly from the accelerations, the modified-update-stress-last (MUSL) scheme 2 computes the particle velocities first and then maps the information to the DOFs. Since USL can generate ill-conditioned mass matrices when piecewise-linear basis functions are used, we consider the MUSL scheme. Therefore, in the next step, the velocity of particle p at time t + Δt is determined as v t+Δt The velocities at the DOFs are subsequently obtained from This allows for computing the incremental displacement vector as After these steps, the remaining part of the particle properties is updated where u p and Δ p are the particle displacement and incremental strain, respectively. The particle stress at time t + Δt is computed from t and Δ t+Δt p using a constitutive law. Since only one-dimensional elastic deformations are considered, it follows from Equations (2), (3), and (6) that Finally, the volume and density of each particle are computed from

Piecewise-linear basis functions
The standard MPM uses piecewise-linear (P1) basis functions illustrated in Figure 1. On the one hand, P1 basis functions have several advantages.

P.2.
Each i has compact support. For i = 2, 3, … , N n − 1, i is supported by element i − 1 and i. The compact supports of basis functions 1 and N n are contained in element 1 and N n − 1, respectively. P.3. Each i is nonnegative over its entire support. P.4. They allow for row-sum mass lumping due to P. 3. P.5. They can be implemented in a straightforward manner.
On the other hand, the gradients of the piecewise-linear basis functions, just as the gradients of all other C 0 -continuous basis functions, are discontinuous on the element boundaries, which can lead to unphysical oscillations in the internal forces when material points cross those boundaries. The gradient of a P1 function is shown as ∇ in Figure 1.

Dual domain material point method
Here, DDMPM 13 uses piecewise-linear basis functions but replaces their gradients with smoother ones. For the construction of the new gradients, the method introduces a weight function and the gradient∇ i defined as where V j is the volume associated with node j and given by V = ∫ Ω dΩ. The weight function is required to be zero on the cell boundaries but is not uniquely specified. 13 For the one-dimensional problems considered in this paper, the following expression is adopted: In addition, DDMPM substitutes the gradients of the piecewise-linear basis functions by The gradients used in this method are depicted in Figure 1. From Section 3.1, it follows that DDMPM modifies the computation of the internal forces and strains.

B-SPLINE MATERIAL POINT METHOD
Here, BSMPM 14,15 replaces the piecewise-linear basis functions from Section 3 by higher-order B-splines, which guarantee at least C 0 -continuity of the gradients. The use of higher-order B-splines reduces not only the grid-crossing error but also the interpolation and time-stepping errors. 11 Similarly to piecewise-linear basis functions, B-splines satisfy properties P. 1.-P. 4. but are less straightforward to implement. Therefore, this paper adopts the BSMPM version proposed by Tielen et al 15 that reduces the implementation complexity by defining B-splines recursively following the Cox-de Boor formula. 24 B-spline basis functions are constructed based on a knot vector, a sequence of ordered nondecreasing points in ℝ called knots. A knot vector is denoted by Ξ = { 1 , 2 , … , n + l + 1 } with n and l being the number of basis functions and the polynomial order, respectively. If the knots are distributed equidistantly, they are said to be uniform. Otherwise, the knots are nonuniform. When more than one knot is positioned at the same location, the knots are called repeated. In an open knot vector, the first and last knots are repeated l + 1 times, which ensures that the resulting basis functions are interpolatory at the domain boundaries. A nonempty knot interval [ i i + 1 ) is referred to as a knot span. For an open uniform knot vector, the number of spans is equal to n − l.
The Cox-de Boor formula defines B-spline basis functions recursively, starting with piecewise constants (no repeated knots, ie, l = 0) For l > 0, the basis functions are given by whereΩ is the parametric domain. The gradients of the B-spline basis functions can be computed from 24 B-spline basis functions satisfy the following properties.
P.1. They form a partition of unity:  4. They allow for row-sum mass lumping because of P. 3.

FUNCTION RECONSTRUCTION
An important property of MPM, DDMPM, and BSMPM is that they preserve the total mass  and linear momentum  of the system. Equation (7), Equation (9), and the partition of unity property of the basis functions yield The time index is dropped to simplify the notation. Moreover, MPM can be viewed as a version of FEM, where the particles provide data for the background grid. In fact, the scattered material-point information is projected to the grid by reconstructing a function of interest and evaluating it at the nodes. Sulsky and Gong 8,9 point out that MPM reconstructs functions using a particular case of Shepard interpolation 7 that can introduce significant errors. They improve the accuracy of the method by adopting higher-order reconstruction techniques. In the modified versions, the velocity, density, and stress fields are reconstructed from the particle data and evaluated either at the nodes or element centres. The MPM-integration is then replaced by a one-point quadrature rule. Moreover, Tielen et al 15,25 introduce an alternative reconstruction technique for BSMPM. They reconstruct functions using cubic-spline interpolation and integrate with a two-point Gauss rule on the half of each nonzero interval. However, it should be noted that while the adjusted mapping techniques can increase the accuracy of the solution, they do not necessarily conserve the total mass and momentum of the system.
In this section, we introduce a novel reconstruction technique, called Taylor least squares (TLS), which maintains the conservative properties of the standard algorithm. Its fundamental concepts are described in Sections 5.1 to 5.3. The use of the TLS technique for mapping the particle data within the material point methods is outlined in Section 5.4. A mathematical analysis of the conservation properties of the technique is provided in Section 5.5.

Least squares approximation
Given a set of N p distinct data points, {x p } N p p=1 and the data values of these points, The least squares 16 approximation at a point x ∈ ℝ is the value w * ∈ P that minimizes, among all w ∈ P, the least squares error Using the basis function vector, (x) = [ 1 (x) 2 (x) · · · n b (x)] T and the vector of unknown coefficients, a = [a 1 a 2 · · · a n b ] T , the least squares approximation can be written as In order to compute the coefficient vector, E a i is set to zero for i = 1, 2, … , n, leading to the normal equations Therefore, we obtain the following expression: Defining the least squares solution is given by What remains to be specified is the basis for P.

Taylor basis functions
A viable choice that leads to an overall conservative reconstruction scheme are the local Taylor basis functions. 17 To define these basis functions, we introduce the concept of the volume average of a function f over the cell e where |Ω e | is the volume of cell e. In one dimension, Ω e = [x min , x max ] with x max > x min , and |Ω e | = x max − x min .
The Taylor basis functions are then given by Here, x c = x max +x min 2 is the cell centroid x c of the cell e, and Δx = . An important quality of the Taylor basis that will ensure the conservation property of the reconstruction technique is

Taylor least squares reconstruction
The TLS approach uses local Taylor basis functions for the least squares approximation of a function f Suppose that ∫ Ω e (x) dΩ e = c with c ∈ ℝ should be conserved by the reconstruction. Then, using Equation (12), we obtain Therefore, the reconstruction procedure preserves the integral quantity if we set It should be noted that Equation (14) can be enforced explicitly. We illustrate this property by reconstructing (x) = sin(x) + 2 on [0, 4 ]. In this case, the integral is equal to 8 . The domain is divided into four elements of size and contains 11 data points. Two data points are located at the boundaries of the first element, (ie, 0 and ). In [2 , 3 ], the data points are distributed uniformly in the interior of the domain. The remaining data points have random positions creating different types of data distribution within an element. It should be noted that the data point at the right boundary of the first element is used for both the first and second elements. The TLS approximation is obtained using three Taylor basis functions. We compare its performance with that of the cubic-spline reconstruction in terms of the root-mean-square (RMS) error for function f and the relative error for the integral of f. The RMS error is computed using 100 Gauss points per element, whereas for the numerical integration, the reconstructed function is evaluated only at two Gauss points within each element. Figure 3 visualizes the data point distribution and the cubic-spline and TLS reconstructions of f for 10 Gauss points per element, whereas Table 1 provides the corresponding errors. Table 1 shows that, for the considered example, the TLS technique outperforms the cubic-spline reconstruction when the conservation and accuracy properties are considered. In fact, the TLS approach preserves the integral up to machine precision. However, Figure 3 shows that the performance of the TLS technique depends on the distribution of the data points within each element.
In some rare cases, data distribution can locally decrease the quality of the TLS approximation but has little influence on the cubic-spline interpolation. An example is provided in Figure 4, where [ , 2 ] contains only two data points located at and 5∕3 . This particular data distribution leads to a linear dependence between the columns of matrix D from   Equation (11) and hence distorts the TLS approximation within this interval. The condition number of D can be used to detect the data distributions that decrease the accuracy of the TLS technique. In addition, the quality of the TLS technique in such situations can easily be improved. First of all, the singularity of D can be prevented by reducing the number of basis functions used for the reconstruction on [ , 2 ]. Although this strategy will preserve the conservative properties of the TLS technique and can be implemented in a straightforward manner, it will lower the accuracy of the method. Therefore, we suggest an alternative approach that uses the information from the neighboring intervals to maintain the high quality of the reconstruction scheme. This approach evaluates the TLS approximation of f associated with [2 , 3 ] at 2 and adds the obtained value to the set of data points, upon which the TLS approximation is based on [ , 2 ]. However, this so-called virtual data point at 2 is excluded from the computation of a 1 in Equation (14). Thus, virtual points do not influence the conservative properties of the technique. Figure 5 illustrates the improved approximation.

Mapping of particle data
When the TLS reconstruction is considered as part of a material point method, particles serve as data points. To conserve the integral of a certain quantity within each element, the coefficient of the first basis function is specified according to Equation (14). The remaining coefficients are calculated from Equation (10) excluding 1 , thereby leaving the integral value unchanged. When the conservation of the reconstructed quantity is not required, a standard least-square approach is followed. This implies that all coefficients are treated as unknowns and their values are obtained from Equation (10).
A TLS reconstruction is applied to replace the MPM integration in Equations (8) and (9) by an exact method such as an elementwise Gauss quadrature. We obtain the approximations with a quadratic TLS reconstruction. This implies that only the first three Taylor basis functions are used (ie, n b = 3). In this case, a two-point Gauss rule within each element leads to an exact integration. The internal forces at the DOFs are computed as follows: 1. Apply a quadratic TLS approach to reconstruct the stress field from the particle data within each active element without specifying the coefficient of the first Taylor basis function: where s i is the coefficient corresponding to Taylor basis function i. Outside of element e,̂e is zero. The global approximation of the stress function,̂, is then equal tô 2. Integrate the stress approximation using a two-point Gauss quadrature where N g is the total number of Gauss points, x g is the global position of a Gauss point, and g is its weight.
The material-point velocities are mapped to the DOFs in the following manner.
1. Apply a quadratic TLS approach to reconstruct the density field and the product of density and velocity from the particle data within each active element, while preserving the mass and momentum of the element: 2. Integrate the approximations using a two-point Gauss quadrature to obtain the momentum vector p and and the consistent mass matrix M C It should be noted that M C may be replaced by a lumped mass matrix without loss of the conservation property of the algorithm.

Compute the velocity vector
From Section 5.3, it follows that the optimal performance of the TLS technique requires at least three particles in each element at the beginning of the simulation. Since this particle distribution is not preserved under large deformations, virtual data points may be used to improve the accuracy of the approximation within elements that contain only one or two material points.

Conservation of mass and momentum
As mentioned in Section 5.4, the TLS technique reconstructs the density field and the product of density and velocity inside each element in the following way: According to Equations (13) and (14), this preserves the mass  e and momentum  e of element e. As a result, the total mass and momentum of the system are conserved after the TLS reconstruction The mass-and momentum-conservation property of the mapping obtained using TLS reconstruction and Gauss quadrature can be shown as well.
Since the total mass is equal to the sum of the entries in the mass matrix from Equation (15), it can be written as The last equality is derived using the partition of unity property of piecewise-linear and B-spline basis functions. In the remaining part of the proof, we assume that n b ≤ integration points per element (or knot span), so that the Gauss quadrature is exact. Therefore, the following holds: The last two steps emerge from the conservation of mass per element and Equation (16).
For the linear momentum, we also assume that n b ≤ 2N g N e and that there are integration points per element (or knot span). Following the above steps, the total momentum after the mapping can be written as Therefore, we have shown that if the Gauss quadrature is performed using a sufficient number of integration points per element, the mass and momentum balance is satisfied not only by the TLS function reconstruction but also by its combination with the Gauss quadrature.

NUMERICAL RESULTS
We study the conservation property of the material point methods by calculating the maximum relative errors in the total mass and momentum over all time steps before and after the computation of the velocity at the DOFs. For MPM, DDMPM, and BSMPM, the errors in the mass and momentum are bounded by 1 · 10 −15 for the vibrating bar benchmark, and 1 · 10 −13 for the column compaction example. Therefore, the TLS results are only compared to those obtained with the cubic-spline reconstruction.

Vibrating bar
This example describes the vibration of a one-phase bar with fixed ends. The motion triggered by an initial velocity that varies along the bar is captured by Equations (1)-(4) with g = 0 (ie, the gravitational force is neglected) and the following initial and boundary conditions: For small strains, the analytical solution in terms of displacement, velocity, and stress is given by The material-point solutions are considered at the particle positions. Table 2 provides exemplary parameter values for the vibrating bar benchmark under small deformations. For the spatial convergence analysis, we minimize the contribution of temporal errors by using an artificially small time-step size and a short simulation time. In fact, the time-step size and total simulation time of, respectively, 1 · 10 −7 s and 1.9 · 10 −6 s eliminate the contribution of the temporal errors and are significantly increased for the grid-crossing study presented later in this section. Furthermore, the number of elements (knot spans) is varied from 5 to 40, while the initial number of particles per cell (PPC) is fixed to 12. These settings ensure that grid crossing does not occur, and the maximal observed strain is equal to 5.3 · 10 −7 m.
The results in terms of spatial errors are shown in Figure 6. As expected, MPM with piecewise-linear basis functions demonstrates second-order convergence in both the displacement and velocity. Since the stress is not discretized directly but is computed from the displacement by taking the derivative, its convergence rate is one. The application of the TLS reconstruction technique, as well as the cubic-spline interpolation, has almost no influence on the stress but decreases the displacement error by a factor of 1.7. For DDMPM, the application of the reconstruction techniques tends to reduce not only the error in the displacement but also in the stress. However, velocity results are similar to those obtained for MPM. The use of quadratic B-splines as basis functions leads to a significant decrease in the error and a higher convergence order for the velocity but causes problems at the boundaries of the domain for both stress and displacement. The absolute error in the stress over the domain is shown in Figure 7. The large values of the error at the boundaries prevent the reduction of the RMS error and worsen the convergence properties of the method. However, the use of BSMPM with a function reconstruction techniques eliminates the boundary issues. An example is provided in Figure 8. Consequently, third-order convergence is obtained for all the considered quantities. It should also be noted that the integration of the TLS or spline reconstruction in BSMPM produces more accurate results than the other considered methods. Table 3 compares the relative error in the mass and momentum made using the TLS and cubic-spline reconstructions. The results are provided for MPM, DDMPM, and BSMPM applied to the vibrating bar problem discretized by 40 elements (knot spans) and 12 PPC. They demonstrate that while the cubic-spline interpolation tends to accurately conserve the mass, it produces errors of order 10 −6 for the linear momentum. The errors produced by the TLS approach consistently   Max. initial velocity v 0 0.80 m/s 2 remain close to machine precision and hence are orders of magnitude smaller than those generated by the cubic-spline interpolation.
For large-deformation simulations, the parameters from Table 4 are used. The time-step size and the simulation time are increased to 1 · 10 −4 s and 0.1 s, respectively. The domain is discretized using 20 elements (knot spans) with initially 8 PPC. The maximal strain that is reached is 0.056 m. Since the analytical solution is not available when the vibrating bar experiences large deformations, the numerical results are compared to the solution obtained with the updated Lagrangian finite element method (ULFEM). 26 In the standard MPM simulation, material points cross the element boundaries more than 100 times, leading to significant inaccuracies in the results. Although grid crossing influences the computation of the displacement and velocity, its most evident consequences are in the stress distribution. DDMPM and BSMPM reduce the grid-crossing error, but their results still significantly deviate from the solution provided by ULFEM. This is shown in Figure 9. Figure 9 also illustrates that the application of the TLS approximation or the cubic-spline reconstruction leads to good agreement of the MPM, DDMPM, and BSMPM solutions with that of ULFEM. The maximal reduction of the relative error in L 2 -norm made by standard MPM is achieved when the reconstruction techniques are combined with BSMPM. More precisely, the integration of the spline interpolation or the TLS reconstruction in BSMPM decreases the MPM error by a factor of 13.5 and 9.9, respectively. On the other hand, the conjunction of the reconstruction techniques with DDMPM leads to highly accurate results as well. The spline and TLS reconstruction reduces the MPM error by a factor of 8.9 and 7.8, respectively. As expected, the TLS approximation conserves the total mass and linear momentum significantly more accurately than the spline interpolation, regardless of the material-point method. The conservative properties of the reconstruction techniques are provided in Table 5.

Column compaction
This example considers the response of a column of material to a compressive body force. The body force increases linearly with time during the simulation. A detailed description of the problem that includes the analytical solution is provided by Bardenhagen and Kober. 10     Parameter values used for the simulation are listed in Table 6. The time-step size and total simulation time are 1 · 10 −4 s and 0.5 s, respectively. Moreover, the number of elements (knot spans) is equal to 50, while the initial number of particles per cell is set to 20. When the standard MPM is used, material-points cross element boundaries 180 times, and the maximal strain is equal to 0.27 m.
The obtained results are shown in Figure 10. The influence of grid-crossing error on the solution is clearly evident for the standard MPM. While DDMPM only slightly smoothens the oscillations, BSMPM significantly improves the results. In fact, the use of B-spline basis functions reduces the relative error in L 2 -norm by a factor of 10.
Both the TLS reconstruction technique and cubic-spline reconstruction interpolation provide considerably more accurate results than MPM and DDMPM. They also further improve the accuracy of BSMPM, but their influence is less distinct in this case. It should be noted that the error produced with the cubic-spline reconstruction is slightly lower than that achieved with the TLS approach. However, the TLS technique outperforms the former method by accurately conserving the mass and linear momentum. As was mentioned previously in this section, the maximum relative errors of MPM, DDMPM, and BSMPM are bounded by 1 · 10 −13 . While the TLS approach has little influence on these errors, the spline reconstruction technique increases them to 2.04 · 10 −6 ( Table 7).
Finally, for DDMPM, the TLS mapping is compared with the conservative integration technique proposed by Dhakal and Zhang 18 that is based on subpoints. In Table 8, we present the relative error and computational time for the DDMPM with a varying number of subpoints and compare it with the DDMPM-TLS result from Figure 10.
From the Table, it follows that, for a low number of subpoints, the technique of Dhakal and Zhang effectively decreases the relative error of 4.4555 · 10 −2 produced by DDMPM and is less computationally intensive than the TLS technique. However, its performance may be limited by computational time and memory when a higher accuracy has to be reached. In fact, the method approaches the accuracy of the TLS reconstruction technique when 128 subpoints per particle are used, leading to a significant (factor of 14.4) increase in the computational time.

SUMMARY AND CONCLUSIONS
In this paper, we have introduced a novel TLS reconstruction approach and applied it to MPMs. The proposed technique combines the least squares approximation with local Taylor basis functions to accurately reconstruct the quantities of interest (eg, stress and density fields) from scattered particle data within each element. We have shown that, in contrast to standard reconstruction techniques, the TLS approximation conserves the mass and linear momentum of the system after the material-point data is mapped to the integration points. More importantly, in conjunction with a sufficiently accurate numerical quadrature method, the technique preserves the total mass and momentum after the information is projected to the degrees of freedom of the grid. This implies that the TLS reconstruction maintains the physical properties of the standard MPM.
The TLS technique was applied to MPM, DDMPM, and BSMPM within the vibrating bar and column compaction problems. For the vibrating bar problem in the absence of grid crossing, the reconstruction technique had little influence on MPM and DDMPM but was able to improve the convergence properties of BSMPM. When material points started to cross cell boundaries in both considered examples, the TLS approximation smoothened the solutions of MPM and related methods and brought them closer to the solution computed analytically or by ULFEM.
Furthermore, the novel TLS reconstruction techniques was compared with the global spline reconstruction technique previously used within BSMPM. In general, the differences in spatial accuracy of the techniques were negligible, but the error in total mass and linear momentum was consistently much lower for the TLS reconstruction. For the column compaction problem, it was also shown that the TLS technique integrated into DDMPM reaches a higher accuracy than the conservative DDMPM method with subpoints while requiring lower computational costs. Therefore, this study has demonstrated that the combination of MPM and related methods with the TLS technique leads to a higher accuracy of the material point methods and preserves their fundamental physical properties.