Matrix-Free Higher-Order Finite Element Methods for Hyperelasticity

TL;DR

This work develops a matrix-free finite element solver for finite-strain hyperelasticity, leveraging an $hp$-multigrid preconditioner to reduce memory traffic and accelerate high-order discretizations in biomechanics. It presents numerically stable reformulations for compressible, nearly incompressible, and fiber-reinforced constitutive models, and analyzes alternative formulations in material vs spatial configurations, including strategies to precompute and store integration-point data. The proposed $hp$-multigrid with matrix-free operator evaluation demonstrates robust performance improvements over matrix-based AMG, especially at higher polynomial degrees, and shows practical speed-ups in patient-specific biomechanics on realistic geometries. The work also provides detailed guidance on storage choices (scalars vs tensors) and numerical techniques (stable $J$, $J^{-2/3}$, $\ln J$, and fast $\exp$) that are critical for exploiting modern hardware, suggesting significant impact for large-scale biomedical simulations and beyond. These contributions offer a viable pathway to faster, memory-efficient simulations of anisotropic, nearly incompressible tissues with realistic constitutive laws.

Abstract

This work presents a matrix-free finite element solver for finite-strain elasticity adopting an $hp$-multigrid preconditioner. Compared to classical algorithms relying on a global sparse matrix, matrix-free solution strategies significantly reduce memory traffic by repeated evaluation of the finite element integrals. Following this approach in the context of finite-strain elasticity, the precise statement of the final weak form is crucial for performance, and it is not clear a priori whether to choose problem formulations in the material or spatial domain. With a focus on hyperelastic solids in biomechanics, the arithmetic costs to evaluate the material law at each quadrature point might favor an evaluation strategy where some quantities are precomputed in each Newton iteration and reused in the Krylov solver for the linearized problem. Hence, we discuss storage strategies to balance the compute load against memory access in compressible and incompressible neo-Hookean models and an anisotropic tissue model. Additionally, numerical stability becomes increasingly important using lower/mixed-precision ingredients and approximate preconditioners to better utilize modern hardware architectures. Application of the presented method to a patient-specific geometry of an iliac bifurcation shows significant speed-ups, especially for higher polynomial degrees, when compared to alternative approaches with matrix-based geometric or black-box algebraic multigrid preconditioners.

Figures (7)

• Figure 1: Visualization of a multigrid V-cycle (left) and an exemplary $hp$-multigrid hierarchy (right), with degrees of freedom indicated by circles.
• Figure 2: Relative errors in the stress and directional derivative (D/Du stress) for sampled ${\mathrm{Grad}\,}\text{\boldmath${u}$}$ comparing standard and stable formulations in the material (${\Omega_0}$) and spatial (${\Omega_t}$) configurations. The nearly incompressible neo-Hookean and fiber models yield similar results, as the iNH model is one part of the fiber model.
• Figure 3: Relative errors in the stress and directional derivative (D/Du stress) for sampled ${\mathrm{Grad}\,}\text{\boldmath${u}$}$ comparing standard and stable formulations in the material (${\Omega_0}$) and spatial (${\Omega_t}$) configurations with fast evaluation strategies.
• Figure 4: Throughput (top row) and memory transfer (bottom row) using various material models considering integration over the reference configuration (${\Omega_0}$) or the spatial configuration (${\Omega_t}$) and precomputing i) no data, ii) scalar values, or iii) tensorial quantities. Throughput for residual evaluation integrating over ${\Omega_t}$ yields values of 5--8$\times10^7$ DoF/s.
• Figure 5: Discretizations of the iliac bifurcation using $l=0,1,2$ refinement levels and mapping degree $p=1,2$.