GPU Matrix-Free Finite Elements Now Handle Moving Quadrature Points Efficiently

Simulations of solids that stretch, crack, or slam into other objects have long posed an awkward problem for finite element practitioners: the very points where the physics is evaluated refuse to stay put. The Material Point Method, a staple of large-deformation solid mechanics, treats material points as integration points that drift freely through a fixed background mesh.

Matrix-free finite element operators are the workhorses of modern supercomputing. Instead of assembling and storing enormous sparse matrices, they recompute the action of an operator on the fly, exploiting structure in the element basis to cut storage to constant cost per degree of freedom and computation to linear cost in the polynomial order. The payoff is higher arithmetic intensity, measured in floating-point operations per byte transferred from memory, and lower memory bandwidth demands per degree of freedom compared with assembled sparse matrices. These techniques, pioneered by the spectral element community and now embedded in libraries such as MFEM and deal.II, achieve very high utilization of modern CPUs and GPUs, which is precisely where the new work aims to deliver.

The difficulty is that classical matrix-free machinery leans on tensor product structure. On hexahedral elements, interpolation and differentiation can be performed as a sequence of one-dimensional contractions, with quadrature points sitting at fixed, consistent locations on every element. Material points in MPM enjoy no such discipline: they move with the material, once per time step, to wherever the physics sends them. Related methods such as CutFEM and other unfitted finite element approaches face the same need to evaluate fields and their spatial derivatives at arbitrarily located quadrature points. The authors treat MPM, following Coombs and colleagues, as a finite element method whose integration points are simply allowed to move independently of the mesh.

Their solution threads the arbitrary points through the operator decomposition of libCEED, a lightweight library providing a purely algebraic interface for linear and nonlinear operators and preconditioners with element-based discretizations. libCEED separates common components, such as reference element configurations, from per-element data like geometric transformations, and selects implementations optimized at run time for CPUs and GPUs, including just-in-time compilation, with interfaces in C, Fortran, Python, Julia, and Rust. The team built their matrix-free MPM operators on this foundation for both CPU and GPU architectures.

The mathematical trick is an intermediate Chebyshev basis. On each reference element, the authors form a matrix that maps nodal solution values to Chebyshev polynomial coefficients; because this matrix is identical for every element in the domain, its product with the standard interpolation matrix can be precomputed and applied in the same number of floating-point operations as ordinary interpolation to fixed quadrature points. From there, evaluating Chebyshev polynomials and their derivatives at each material point, wherever it happens to sit, completes the transfer. The reference coordinates of each material point are precomputed by a pointwise Newton iteration and refreshed once per time step. Because quadrature weights of 1.0 must be used owing to the arbitrary point locations, the integration weights absorb the Jacobian determinant of the mapping instead.

The computational price is quantified precisely. Interpolating values or derivatives at standard quadrature points on tensor product elements costs order q to the fourth power floating-point operations, assuming the polynomial order does not exceed the quadrature order. Doing the same at material points costs order q to the sixth power when the number of material points matches the number of quadrature points, and order q cubed times the number of points per cell when material points outnumber quadrature points. That is a steep surcharge, but the authors note it lands on the same asymptotic order as basis evaluation on elements without tensor product structure. Their scheme also matches the asymptotic performance reported independently by Bergbauer and colleagues, whose similar derivation differs mainly in running on CPUs only and in counting fused multiply-add operations as single operations.

GPU implementation adds its own choreography. Computing Chebyshev coefficients requires information from the whole element, but evaluating polynomial values and interpolated solutions is independent at each material point, so material points can be distributed across threads while the tensor basis operation runs per element. When material points per cell are non-uniform, as is common in MPM, cells are padded with virtual material points whose calculations are computed and then discarded, tying per-cell cost to the maximum point count in the partition. The CPU implementations skip this padding entirely.

Benchmarks from the Center for Efficient Exascale Discretizations put the numbers on the table, run on four AMD Instinct MI300A accelerated processing units, each pairing 24 CPU cores with 228 GPU compute units, 128 gigabytes of unified memory, a peak of 61.3 trillion double-precision operations per second, and 5.3 terabytes per second of memory bandwidth. On the BP2 mass-matrix benchmark, at-points evaluation demanded an order of magnitude more floating-point operations than tensor product evaluation, and linear and quadratic elements lost roughly one order of magnitude in efficiency, with cubic elements degrading further. On the BP4 Poisson benchmark, at-points throughput kept climbing with problem size, peaking near 100 trillion operations per second for quadratic elements, and linear and quadratic at-points operators outperformed non-tensor basis evaluation, a result the authors flag as promising for future multigrid methods under hp-refinement.

The most telling test is a large-deformation static elasticity problem using the Ogden constitutive model, solved on a unit cube under uniaxial compression. Where the Poisson benchmark spends only 45 floating-point operations per quadrature point on the physics, the Ogden linearization spends 1675, and the balance of work shifts dramatically: tensor product operators need 297 operations per quadrature point for the basis transfers while at-points operators need 5443, so the total operation counts converge. Throughput for tensor and non-tensor operators rose sharply as the heavier physics kept the GPU busy, while at-points throughput dipped slightly, yet at-points evaluation retained its efficiency advantage for quadratic elements, which the authors highlight as particularly relevant for finite-strain material point methods, where quadratic elements help prevent artificial locking.

The authors are candid about limits and next steps. Cubic element performance suffers primarily from the enriched quadrature choice used in the benchmarks, and the more common choice of equal quadrature and polynomial order avoids that degradation, while higher quadrature orders remain an optimization target. Production problems typically carry many material points per cell because arbitrarily located points integrate poorly, and performance holds when material points increase at low quadrature order. Robust preconditioners are still needed before at-points operators enter production workflows, though libCEED and PETSc already allow the underlying Ratel library to run on unstructured meshes. With scripts, inputs, and consolidated data publicly archived, the work opens a path to matrix-free MPM simulations of large deformations, history-dependent materials, and complex contact and fracture interactions on exascale-class hardware.

Subject of Research: Matrix-free finite element evaluation at arbitrary quadrature point locations for material point method simulations on GPU hardware

Article Title: Matrix-free finite element methods with arbitrary quadrature point locations

Article References: Atkins, Z. R., Brown, J., Shakeri, R., & Thompson, J. L. (2026). Matrix-free finite element methods with arbitrary quadrature point locations. Computational Science and Engineering, 3(1), Article 3. https://doi.org/10.1007/s44207-026-00010-1

Image Credits: AI Generated

DOI: 10.1007/s44207-026-00010-1

Keywords: finite element methods, material point method, matrix-free operators, GPU computing, libCEED, quadrature points, Chebyshev polynomials, solid mechanics, high-performance computing, CEED benchmarks, Ogden model, tensor product elements