Unstructured MLS Material Point Methods

• Unstructured Moving Least Squares Material Point Methods are a hybrid Eulerian-Lagrangian simulation technique providing continuous C1 interpolation on unstructured simplicial grids.
• They modify standard MLS weights with diminishing factors to ensure smooth gradient reconstruction and second-order spatial accuracy even under large deformations.
• They balance computational efficiency and simulation precision, making them ideal for complex geometries in applications like geomechanics and biomedical simulations.

Unstructured Moving Least Squares Material Point Methods (UMLS-MPM) provide a hybrid Eulerian-Lagrangian simulation technique for solid mechanics under significant deformation, designed to overcome the limitations of classical Material Point Methods (MPM) on unstructured grids. Conventional MPM typically utilizes $\mathcal{C}^0$ finite-element basis functions on structured quadrilateral or hexahedral grids, which are insufficient for handling complex geometries due to discontinuous velocity gradients and resulting cell-crossing errors. UMLS-MPM introduces a stable moving least squares (MLS) kernel with diminishing sample weights to achieve continuous ($\mathcal{C}^1$) shape functions and velocity gradients over general unstructured triangulations (2D) and tetrahedral tessellations (3D) (Cao et al., 2023).

1. Motivation and Core Challenges

Standard updated-Lagrangian MPM architectures are effective for simple domains using structured grids, but encounter two principal obstacles on unstructured simplices:

• Cell-Crossing Error: Material points crossing from one simplex to another encounter a discontinuous gradient in the interpolation weight function ($\nabla w$), generating spurious stress oscillations.
• Geometric Complexity: Unstructured triangulations or tetrahedralizations are necessary to match complex boundaries (e.g., soil-rock interfaces, biological tissue), but the absence of gradient continuity impairs accuracy, especially under large deformation (Cao et al., 2023).

Previous attempts to construct $\mathcal{C}^1$-continuous interpolation functions have either been limited to structured backgrounds, have not supported unstructured backgrounds, or have not been generalized beyond 2D triangular meshes. The UMLS-MPM framework addresses these limitations by enforcing analytic continuity in the MLS kernel itself.

2. MLS Kernel Construction with Diminishing Weights

The MLS approximation for a field $u(x)$ uses $n$ grid-vertex samples $\{x_v, u_v\}$ from the 1-ring neighborhood of $x$. The linear basis is $p(\xi) = [1, \xi^{\mathsf{T}}]^\mathsf{T}$, $\xi = x_v - x$, with initial sample weight $\mathcal{C}^1$0 typically chosen as a cubic B-spline in $\mathcal{C}^1$1.

2.1 Standard MLS Formulation

• Moment Matrix: $\mathcal{C}^1$2.
• Right-Hand Vector: $\mathcal{C}^1$3.
• MLS Solution: $\mathcal{C}^1$4, yielding $\mathcal{C}^1$5 and $\mathcal{C}^1$6 via $\mathcal{C}^1$7.

2.2 Diminishing-Weight Modification

To ensure $\mathcal{C}^1$8 continuity as the sampling neighborhood $\mathcal{C}^1$9 evolves (especially when crossing simplex boundaries):

• Weights are modified by a diminishing factor $\nabla w$0:

$\nabla w$1

• For simplices,

$\nabla w$2

where $\nabla w$3 are vertices of the host simplex, $\nabla w$4 is the barycentric coordinate, and $\nabla w$5 encodes adjacency.

• By construction, $\nabla w$6 for $\nabla w$7 vertices and diminishes to zero for vertices entering or leaving the 1-ring as $\nabla w$8 crosses edges/faces, ensuring smoothness.

2.3 MLS Shape Function and Gradient

With these modified weights:

• Shape function:

$\nabla w$9

• Gradient:

$\mathcal{C}^1$0

$\mathcal{C}^1$1 is computed analytically, resulting in continuous gradients across element boundaries.

3. Continuous Gradient Reconstruction and Theoretical Continuity

As a material point traverses a simplex boundary, the neighborhoods ($\mathcal{C}^1$2, $\mathcal{C}^1$3) of added and removed nodes shift, but for each $\mathcal{C}^1$4 in these sets, $\mathcal{C}^1$5, with $\mathcal{C}^1$6 the crossing distance. Perturbation analysis shows that the disturbance in $\mathcal{C}^1$7 and $\mathcal{C}^1$8 is also $\mathcal{C}^1$9. Thus, by standard matrix perturbation theory, both the reconstructed field $u(x)$0 and its gradient $u(x)$1 vary only by $u(x)$2. Consequently, $u(x)$3 and $u(x)$4 across all element interfaces, demonstrating full $u(x)$5 continuity (Cao et al., 2023).

A detailed proof utilizing moment matrix perturbations is provided in Appendix A of the source.

4. Extension to 2D and 3D Simplicial Tessellations

The construction generalizes to both triangles (2D) and tetrahedra (3D):

• 2D (Triangles): $u(x)$6 comprises 3 vertices, and the 1-ring $u(x)$7 typically involves 8–12 vertices. $u(x)$8 is calculated by summing barycentric coordinates of 0-ring vertices adjacent to $u(x)$9.
• 3D (Tetrahedra): $n$0 contains 4 vertices, with the 1-ring extending to up to ~30 vertices. The same adjacency-based definition for $n$1 applies.

The B-spline support radius is chosen as the maximal mesh edge length $n$2, ensuring all neighboring vertices within the 1-ring are sampled at each relevant location.

Computational costs per particle involve assembling and inverting a $n$3 moment matrix, with a bounded neighbor count, making the approach practical for $n$4.

5. Algorithmic Implementation Pipeline

The simulation cycle consists of the following principal steps, applied per time step:

1. Neighborhood Identification: Locate the host simplex for each particle and form $n$5 and $n$6 adjacency sets.
2. MLS Weight Computation: For each neighbor vertex, compute $n$7, barycentric coordinates $n$8, $n$9, combined weights $\{x_v, u_v\}$0, and their gradients.
3. Particle-to-Grid (P2G) Projection: Aggregate mass, momentum, and internal/external forces to the grid using MLS weights.
4. Grid Update: Advance grid velocities via time integration and apply boundary conditions.
5. Grid-to-Particle (G2P) Back-Projection: Update particle velocities, positions, deformation gradients, volumes, and stresses, including plasticity routines as necessary.
6. Grid Reset: Reinitialize grid fields for the subsequent time step.

Essential data structures include the unstructured mesh with vertex adjacency lists and an efficient spatial search mechanism (e.g., bounding-box hierarchy or spatial hashing) to locate containing simplices.

6. Numerical Experiments and Convergence Properties

A suite of benchmarks demonstrates the practical properties and convergence rates of UMLS-MPM:

Test Case	Key Outcome	Metric/Result
1D Vibrating Bar	2nd-order convergence, energy conservation	RMSE, <0.1% energy err
2D Colliding Disks	Momentum, kinetic energy match, no oscillation	Stress/time-history
2D Cantilever (mesh rot.)	Errors <5% up to 45°, 2nd-order convergence	RMSE
3D Slope Failure	Agreement in strain/stress vs B-spline-MPM	Run-out distance/stat.
3D Elastic Sphere Expansion	Smooth contact forces, nonlinear boundary	Force profiles

In all cases, UMLS-MPM mitigates cell-crossing error, maintains energy conservation comparable to high-order MPM variants, and demonstrates second-order spatial accuracy consistent with the theoretical error bounds of a $\{x_v, u_v\}$1 MLS interpolant (Cao et al., 2023).

7. Significance and Practical Considerations

The UMLS-MPM framework:

• Eliminates cell-crossing artifacts: Stress oscillations are suppressed by virtue of enforced $\{x_v, u_v\}$2 MLS kernels, maintaining a smooth velocity gradient across arbitrary simplices.
• Ensures accuracy and convergence: Achieves verified second-order accuracy in spatial discretization for smooth problems.
• Maintains stability and energy conservation: Utilizes symplectic time-integration and continuous gradient reconstruction to provide stable updates and nearly conservative energy behavior without special quadrature rules.
• Balances computational cost: The principal overhead arises from per-particle assembly and inversion of small dense systems and neighborhood searches; this is modest compared to the typically dominant constitutive model workloads in large-deformation MPM.

The approach is compatible with any particle–grid codebase capable of supporting unstructured simplicial backgrounds, furnishing both geometric flexibility and high-order accuracy for complex domains (Cao et al., 2023).