MLS-MPM: Moving Least Squares Material Point Method

• MLS-MPM is a hybrid computational method that combines Moving Least Squares reconstruction with MPM for stable and accurate simulation of large-deformation solid mechanics on unstructured meshes.
• It employs continuous MLS shape functions with C1 continuity to effectively reduce cell-crossing errors and achieve second-order spatial convergence even under extreme deformations.
• The method integrates with explicit MPM schemes like FLIP and APIC, conserving linear and angular momentum while efficiently handling complex geometries and adaptive updates.

The Moving Least Squares Material Point Method (MLS-MPM) is a hybrid Eulerian-Lagrangian computational approach for large-deformation solid mechanics. It augments the conventional Material Point Method (MPM) by employing local moving least squares (MLS) reconstruction to devise higher-order, adaptively smooth interpolation kernels for transferring physical quantities between particles and background grid. Recent developments extend MLS-MPM to fully unstructured meshes, generalizing its capabilities to arbitrary simplex tessellations in 2D and 3D, while ensuring continuous gradients and robust convergence—even across nonuniform geometries and extreme deformations (Cao et al., 2023, Su et al., 2021).

1. Mathematical Formulation and Kernel Design

MLS-MPM recasts interpolation and gradient evaluation as localized least-squares problems centered at particle locations. For a particle at position $x$ within an unstructured simplicial mesh $\Omega$, the support of the MLS kernel is defined by the 1-ring vertex neighborhood $\mathcal{N}^1(x)$ (the element vertices and their direct neighbors), mitigating the "kernel degeneration" that occurs when support changes abruptly at cell boundaries.

A function $u(x)$ is reconstructed via a local MLS fit to its samples $u_v$ at vertices $v$ in $\mathcal{N}^1(x)$, using a linear polynomial basis $p(z) = [1, z^\top]^\top$ and a diagonal sample-weight matrix $D(x)$ with entries $D_{vv}(x) = d(x_v - x)$, where $d(\cdot)$ is smooth, compactly supported, and positive (e.g., B-spline profile). The coefficients $a(x)$ minimizing the weighted least squares error,

$$\sum_{v \in \mathcal{N}^1(x)} d(x_v - x) \|u_v - a(x)^\top p(x_v-x)\|^2,$$

yield the MLS reconstruction: $a(x) = M(x)^{-1} b(x),$ where

$$M(x) = \sum_v d(x_v-x) p(x_v-x) p(x_v-x)^\top,\qquad b(x) = \sum_v d(x_v-x) p(x_v-x) u_v.$$

The reconstructed value and gradient at $x$ then follow as

$$[\hat{u}(x);\, \nabla\hat{u}(x)] = M(x)^{-1} \sum_{v\in \mathcal{N}^1(x)} d(x_v-x) p(x_v-x) u_v.$$

To address discontinuous neighborhood transitions on unstructured meshes, the MLS kernel incorporates a diminishing function $\eta_v(x)$, vanishing smoothly as a vertex $v$ exits or enters the support: $$\eta_v(x) = \sum_{n \in \mathcal{N}^0(x)} B_n(x) A_{vn},$$ where $B_n(x)$ are barycentric coordinates of $x$ in its cell, and $A_{vn}$ encodes vertex adjacency. The MLS sample weights become $\bar{w}_v(x) = \eta_v(x) d(x_v - x)$. The MLS shape function and its gradient are given by

$w_v(x) = p(x_v - x)^\top M(x)^{-1} \bar{w}_v(x),$

with the gradient derived by differentiating $w_v(x)$ with respect to $x$ (Cao et al., 2023).

2. Continuity, Cell-Crossing, and Gradient Properties

Standard ($\mathcal{C}^0$) piecewise-linear FEM or MPM basis functions on unstructured meshes yield discontinuous gradients across element faces, inducing "cell-crossing errors" manifested as spurious stress oscillations. The MLS-MPM approach, enhanced with the diminishing factor, achieves $\mathcal{C}^1$ continuity of both the shape functions and their gradients. This is mathematically justified: as a particle transitions from one element to another, added or removed node sets contribute weights that vanish smoothly, ensuring that both $w_v(x)$ and $\nabla w_v(x)$ remain continuous under mesh traversal—a critical requirement for accuracy and stability in high-deformation physics (Cao et al., 2023).

3. Algorithmic Integration: Explicit UMLS-MPM and APIC+MLS

MLS-MPM is integrated within the explicit MPM (FLIP) pipeline, substituting MLS-based functions for classical grid-to-particle and particle-to-grid transfer. The UMLS-MPM algorithm proceeds as follows:

• Particle-to-Grid (P2G): Mass, momentum, and internal/external forces are rasterized using the MLS shape functions $w_{ip} = w_i(x_p)$, with gradients $\nabla w_{ip}$ applied for computing stress divergence.
• Grid Update: Node velocities are advanced via explicit symplectic Euler integration, incorporating computed accelerations.
• Grid-to-Particle (G2P): Updated node velocities and accelerations are interpolated back to particle positions.
• Deformation and Stress Update: The deformation gradient $F_p$, particle volume, and stress tensor $\sigma_p$ are updated using the MLS kernel.
• Reset: Grid nodal quantities are cleared and time advances.

In the context of Arbitrary Updated Lagrangian MPM (A-ULMPM), the MLS kernel naturally interpolates between total-Lagrangian and Eulerian limits. Reference configuration updates are triggered based on a physics-driven criterion involving the deviation of the local Jacobian from unity. Kernel and support recalculations are only performed upon such updates, improving computational efficiency in scenarios with localized large deformations (Su et al., 2021). Velocity rasterization is further augmented by the APIC scheme, conserving angular and linear momentum up to machine precision and leveraging first-order velocity derivatives reconstructed by the MLS operator.

4. Stability, Robustness, and Convergence Properties

Stability derives from the partition of unity ($\sum_i w_{ip} = 1$), non-negativity (if $\eta_v, d \ge 0$ and $M$ is well-conditioned), and the continuity of MLS gradients, which suppresses stress spikes at cell crossings. Well-posedness requires non-degenerate sampling in the 1-ring support: the minimum singular value $\sigma_{\min}(M)$ must remain bounded away from zero. Empirical and theoretical analyses confirm that, under mesh gradations $R \le 1.5$, the mean-square error in particle displacement converges as $\mathrm{RMSE}(h) = O(h^2)$ for coarse grids, with temporal errors $O(\Delta t)$ under the CFL condition (Cao et al., 2023, Su et al., 2021).

5. Numerical Validation and Benchmarks

Comprehensive benchmarks validate the performance and robustness of MLS-MPM and UMLS-MPM:

Test Problem	Key Metrics / Observations	Reference
1D Vibrating Bar	Max RMSE $<0.554\%$ @ $dx=0.5$ m; 2nd-order spatial convergence; $<0.1\%$ energy error	(Cao et al., 2023)
2D Colliding Elastic Disks	$\leq 1\%$ momentum loss vs. high-res B-spline MPM; no spurious oscillations	(Cao et al., 2023)
2D Rotated Cantilever	Tip-displacement error (1.27–4.72%) for varying mesh rotations; maintains 2nd order	(Cao et al., 2023)
3D Slope Failure (Elastoplastic)	Plasticity, mean stress, and run-out histories align with total-Lagrangian MPM	(Cao et al., 2023)
3D Spherical Container Expansion	Conformal handling of detailed curved boundaries, impact, and expansion phases	(Cao et al., 2023)
3D Multi-Material, Fluid/Solid Coupling	Large deformation, angular/lateral momentum conserved, adaptive reference update	(Su et al., 2021)

Within these tests, UMLS-MPM consistently recovers high-order spatial convergence and eliminates cell-crossing artifacts, even for complex 2D and 3D unstructured meshes. Energy and momentum conservation metrics are consistently within strict tolerances; equivalent results on standard benchmarks match or exceed results from B-spline or structured-grid MLS-MPM (Cao et al., 2023, Su et al., 2021).

6. Advantages, Limitations, and Computational Considerations

MLS-MPM and its variants retain several advantages over standard MPM approaches:

• Adaptivity and Robustness: The MLS kernel adapts to deformation and mesh quality, eliminating "numerical fracture" and restoring high-order convergence, even under extreme deformations.
• Generalization to Unstructured Meshes: UMLS-MPM exploits general simplex tessellation in both 2D and 3D, enabling accurate simulation of complex geometries.
• Conservation Properties: Integration with APIC methods enables exact conservation of both linear and angular momentum, up to floating-point precision (Su et al., 2021).
• Computational Profile: Moment matrix inversion per particle is required upon reference reset; otherwise, most kernel operations are localized and efficient for modern accelerator hardware.
• Reference Update Criterion: Physics-based triggers for reference configuration change balance accuracy and computational efficiency, e.g., thresholds $\varepsilon\approx0.01$ for fluids, $0.5$ for soft solids.

Limiting factors include potential ill-conditioning of the MLS moment matrix for very sparse or poorly graded meshes, addressed via regularization or robust solvers, and stringent timestep restrictions for highly stiff materials (Su et al., 2021). Meshes with sufficiently small gradation ratios $R$ and well-defined 1-ring neighborhoods are prerequisites for stability and optimal convergence.

7. Connections, Extensions, and Related Methodologies

MLS-MPM generalizes the standard MPM by adapting shape function support, reproducing both constant and linear fields exactly. It provides a seamless bridge between total Lagrangian and fully Eulerian formulations, unified by the A-ULMPM approach, and is compatible with particle-in-cell schemes such as FLIP or APIC for transfer steps. The method's $\mathcal{C}^1$ kernel construction on unstructured meshes distinguishes it from prior efforts, and the diminishing factor in kernel support offers a general recipe for extending MLS-based approaches to arbitrary tessellations. Applications include elastoplasticity, multi-material interaction, and large-deformation fluid-solid coupling, making MLS-MPM a versatile tool for computational physics and engineering simulation (Cao et al., 2023, Su et al., 2021).

---

References:

• "Unstructured Moving Least Squares Material Point Methods: A Stable Kernel Approach With Continuous Gradient Reconstruction on General Unstructured Tessellations" (Cao et al., 2023)
• "A-ULMPM: An Arbitrary Updated Lagrangian Material Point Method for Efficient Simulation of Solids and Fluids" (Su et al., 2021)