---
title: Unstructured MLS Material Point Methods
url: https://www.emergentmind.com/topics/unstructured-moving-least-squares-material-point-methods
type: topic
---

# Unstructured MLS Material Point Methods

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) [2312.10338].

## 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 [2312.10338].

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 $d_v(x)$ typically chosen as a cubic B-spline in $\|x - x_v\|$. 

### 2.1 Standard MLS Formulation

- **Moment Matrix:** $A(x) = \sum_{v \in N(x)} d_v(x)\, p(x_v - x) p(x_v - x)^\mathsf{T} \in \mathbb{R}^{(d+1) \times (d+1)}$.
- **Right-Hand Vector:** $b(x) = \sum_{v \in N(x)} d_v(x)\, u_v\, p(x_v - x) \in \mathbb{R}^{d+1}$.
- **MLS Solution:** $A(x)a(x) = b(x)$, yielding $u(x)$ and $\nabla u(x)$ via $[a_0(x); a_1(x)] = a(x)$.

### 2.2 Diminishing-Weight Modification

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

- Weights are modified by a diminishing factor $\eta_v(x)$:
  $$
  w_v(x) = \eta_v(x)\, d_v(x)
  $$
- For simplices,
  $$
  \eta_v(x) = \sum_{n \in V^0(x)} B_n(x) A_{v,n}
  $$
  where $V^0(x)$ are vertices of the host simplex, $B_n(x)$ is the barycentric coordinate, and $A_{v,n}$ encodes adjacency.
- By construction, $\eta_v(x) = 1$ for $V^0(x)$ vertices and diminishes to zero for vertices entering or leaving the 1-ring as $x$ crosses edges/faces, ensuring smoothness.

### 2.3 MLS Shape Function and Gradient

With these modified weights:

- Shape function:
  $$
  \phi_v(x) = p(0)^\mathsf{T} A(x)^{-1} B_v(x)
  $$
- Gradient:
  $$
  \nabla \phi_v(x) = -[\partial_x p(x_v - x)]^\mathsf{T} A(x)^{-1} B_v(x)
  - p(0)^\mathsf{T} A(x)^{-1} [\partial_x A(x)] A(x)^{-1} B_v(x)
  $$

$\nabla \phi_v(x)$ 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 ($N_a$, $N_r$) of added and removed nodes shift, but for each $v$ in these sets, $\eta_v(x) = \mathcal{O}(\epsilon)$, with $\epsilon$ the crossing distance. Perturbation analysis shows that the disturbance in $A(x)$ and $b(x)$ is also $\mathcal{O}(\epsilon)$. Thus, by standard matrix perturbation theory, both the reconstructed field $\hat{u}(x)$ and its gradient $\nabla \hat{u}(x)$ vary only by $\mathcal{O}(\epsilon)$. Consequently, $\phi_v(x) \in \mathcal{C}^0$ and $\nabla \phi_v(x) \in \mathcal{C}^0$ across all element interfaces, demonstrating full $\mathcal{C}^1$ continuity [2312.10338].

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):** $V^0(x)$ comprises 3 vertices, and the 1-ring $V^1(x)$ typically involves 8–12 vertices. $\eta_v(x)$ is calculated by summing barycentric coordinates of 0-ring vertices adjacent to $v$.
- **3D (Tetrahedra):** $V^0(x)$ contains 4 vertices, with the 1-ring extending to up to ~30 vertices. The same adjacency-based definition for $\eta_v(x)$ applies.

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

Computational costs per particle involve assembling and inverting a $(d+1) \times (d+1)$ moment matrix, with a bounded neighbor count, making the approach practical for $d=2,3$.

## 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 $V^0(p)$ and $V^1(p)$ adjacency sets.
2. **MLS Weight Computation:** For each neighbor vertex, compute $d_v$, barycentric coordinates $B_n$, $\eta_v$, combined weights $w_{vp}$, 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 $\mathcal{C}^1$ MLS interpolant [2312.10338].

## 7. Significance and Practical Considerations

The UMLS-MPM framework:

- **Eliminates cell-crossing artifacts:** Stress oscillations are suppressed by virtue of enforced $\mathcal{C}^1$ 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 [2312.10338].

Source: https://www.emergentmind.com/topics/unstructured-moving-least-squares-material-point-methods