Saving Gradients on Unstructured Grids — Moving Least-Squares Approximation

The Day a Gradient Broke on an Unstructured Grid#

At two in the afternoon, on an unstructured tet mesh, an LSQ (least-squares) gradient flipped negative in a single cell. Neighboring cells read 0.3 and 0.5, but the middle one read -4.1. I walked back through the mesh: that cell's seven neighbors were at wildly different distances. One nearby point was being weighted the same as six faraway ones in the normal equations.

The fix had already existed since 1981, under a different name in a different field.

Move h, the polynomial order, or the sample count N and the derivative RMSE moves with them.

Shrink $h$ to 0.025 and at some evaluation points there are too few active nodes — $\mathbf{A}$ becomes nearly singular and RMSE explodes. Raise $h$ to 0.20 and the amplitude is averaged out, blunting the derivative. There is a narrow sweet spot between the two.

1981 — A Loan from Surface Fitting#

Lancaster and Salkauskas made a promise in a surface reconstruction paper: "at every evaluation point, refit a polynomial using only data near that point." They called it moving least-squares, MLS. It was a tool for building smooth surfaces in computer graphics.

Twenty years later Belytschko's EFG (Element-Free Galerkin) borrowed the same equation as the backbone of meshless structural analysis. Ten more years and unstructured FVM took it for gradient reconstruction. SPH's consistency-corrected derivative lives in the same place. The shared starting point — points instead of grids — calls for the same equation in every field.

A Weight Function Turns Distance into Influence#

The starting equation is simple. Near an evaluation point $\mathbf{x}^*$, fit the field with a polynomial basis.

$\mathbf{p}$ is the polynomial basis vector (linear 2D: $[1, \xi, \eta]$; quadratic 2D: $[1, \xi, \eta, \xi^2/2, \xi\eta, \eta^2/2]$), and $\mathbf{a}$ is the coefficient vector.

The crucial part: $\mathbf{a}$ is a function of the evaluation point $\mathbf{x}^*$. Move the point and the coefficients are re-solved. That is the "moving" in moving least-squares.

The coefficients minimize the weighted sum of squared residuals:

The weight function $w$ does the heavy lifting. Points close to $\mathbf{x}^*$ get large weights; far points get nearly zero. Gaussian $w = \exp(-(r/h)^2)$, cubic spline, fourth-order polynomial — they all share one trait. Beyond a distance $h$ they cut off to zero. That $h$ is the kernel radius (or support radius).

Normal Equations and the Shape Function#

Setting $\partial J / \partial \mathbf{a} = 0$ produces the normal equations:

with

$\mathbf{u}$ holds the nodal values, $w_i = w(\|\mathbf{x}^* - \mathbf{x}_i\|)$. Substituting back:

$\Phi_i(\mathbf{x}^*)$ is the shape function. It carries a small but annoying property: $\Phi_i(\mathbf{x}_j) \neq \delta_{ij}$. Interpolation is not exact at the nodes. That bites EFG when imposing Dirichlet boundary conditions. For CFD gradient reconstruction it does not matter — we want a smooth derivative, not exact nodal interpolation.

Python — Differentiating sin from 25 Points#

Choose the basis in Taylor form, $p_k(\xi) = \xi^k/k!$, and the coefficients become derivatives directly — $a_k \approx u^{(k)}(\mathbf{x}^*)$.

import numpy as np
from math import factorial
 
def mls_taylor_1d(xs, us, x_star, h, order=2):
    """Weighted LSQ with Taylor basis at evaluation point x_star.
    Returns a[k] approximating u^(k)(x_star).
    """
    m = order + 1
    A = np.zeros((m, m))
    B = np.zeros(m)
    for x_i, u_i in zip(xs, us):
        dx = x_i - x_star
        w = np.exp(-(dx / h) ** 2)
        if w < 1e-10:
            continue
        p = np.array([(dx ** k) / factorial(k) for k in range(m)])
        A += w * np.outer(p, p)
        B += w * p * u_i
    return np.linalg.solve(A, B)
 
# Validation — recover the derivative of sin(2 pi x)
np.random.seed(7)
xs = np.sort(np.random.uniform(0, 1, 25))
us = np.sin(2 * np.pi * xs)
 
x_eval = np.linspace(0.05, 0.95, 200)
du_hat = np.array([mls_taylor_1d(xs, us, xq, h=0.07, order=2)[1] for xq in x_eval])
du_true = 2 * np.pi * np.cos(2 * np.pi * x_eval)
 
rmse = np.sqrt(np.mean((du_hat - du_true) ** 2))
print(f"derivative RMSE: {rmse:.3e}")
# derivative RMSE: 2.1e-01   (h=0.07, order=2)

Twenty-five random points, derivative RMSE around 0.21. The amplitude is $2\pi \approx 6.28$, so the relative error is about 3%. Bumping order to 3 and $h$ to 0.10 drops it further. But push $h$ too high and a different problem appears.

For the Gaussian $w=\exp(-(r/h)^{2})$, the weight falls to one half where $(r/h)^{2}=\ln 2$. At $h=0.07$ that radius is about 0.058. Mean spacing is about 0.042. The snippet drops a point once $w<10^{-10}$, about $4.8h$ out: 0.34 at $h=0.07$, and 0.12 at $h=0.025$.

Order 2 means three unknowns. The $h=0.025$ window is about 0.24 across, only five or six points at a spacing of 0.042. Seed 7 leaves gaps under that average, and $\mathbf{A}$ gets ill-conditioned there. Divide the derivative RMSE 0.21 by the amplitude $2\pi\approx 6.28$ and the ratio is 0.033, about 3%.

Set $h$ to 1.5 times the mean spacing and $0.042\times 1.5\approx 0.063$, almost the snippet's 0.07. Push $h$ to 3 times, about 0.13, and the half-weight radius passes two neighbors, dulling the higher derivative.

How unweighted least squares splits on a skewed cell is written out in Green–Gauss versus weighted least squares on a skewed mesh. The −4.1 that afternoon was a near point and a far point sharing one coefficient in the normal equations.

Kernel Radius h — What Each Extreme Destroys#

In gradient reconstruction, $h$ trades two things against each other.

• Small $h$: the condition number of $\mathbf{A}$ blows up. One or two nearby points carry all the weight, and you try to fit a degree-$m$ polynomial with fewer than $m$ active points. The result becomes noise-sensitive.
• Large $h$: information from distant points mixes in and averages out. The zeroth moment survives, but higher derivatives go dull.

A practical rule: pick $h$ at 1.5–3× the local average node spacing, so the support contains 1.5–2× $m$ active nodes. On unstructured grids you must adapt $h$ per cell — a global $h$ almost always fails.

Scale up to 2D. Twenty-eight scattered nodes carry $u = \sin(2\pi x)\cos(2\pi y)$. Drag the yellow query point and watch the weights shift.

Drop the query into the empty upper-left region and the active count falls below $m$ — $\mathbf{A}$ goes singular. Drag it back into a dense cluster and the error returns to small. This is MLS's structural weakness on unstructured grids, and the reason an adaptive $h$ is non-negotiable in practice.

Echoes in SPH, EFG, and FVM#

MLS is alive in three places today under different names.

• SPH consistency correction: the standard SPH derivative loses even first-order accuracy on non-uniform particle distributions. Corrected SPH and reproducing-kernel SPH inject a weighted least-squares correction to restore first-order consistency — essentially first-order MLS.
• EFG / MLPG meshless structural analysis: the shape function is literally MLS. Non-exact nodal interpolation makes Dirichlet boundaries awkward; Lagrange multipliers or penalty methods are the standard fix.
• Unstructured FVM gradient reconstruction: Compact Least-Squares Reconstruction solves the normal equation over cell centers to extract gradients and even curvature. OpenFOAM's leastSquares is a first-order polynomial MLS variant.

Three fields, one set of normal equations. Only the weight function and the basis dimension change.

Back to the Broken Gradient#

Back to the opening cell. That LSQ gradient flipped to -4.1 because two close neighbors entered the normal equations with the same weight as five distant ones. Add a distance weight and the value moved to 0.4. The divergence vanished.

Locality through a weight function, normal equations re-solved at every point, and the non-separability of the shape function — three pieces MLS lent us from graphics. Today every meshless integration and half of all unstructured integrations stand on that scaffolding.

The cell whose gradient jumps to −4.1 when every neighbor shares one weight#

Two near neighbors and five far ones entered with the same coefficient, and that cell's gradient jumped to −4.1. A distance weight brought it back to 0.4. Cut $h$ far below the node spacing and $\mathbf{A}$ leans singular instead. The sign can look fine while the magnitude shakes again.

Once the gradient is back, the next call is how much of it a cell is allowed to keep. The stretch where a second-order reconstruction builds fake extrema on an unstructured mesh is set down in where a Barth–Jespersen limiter cuts. Hand the slope back with weights and then skip the limiter, and that cell grows an extremum again.