Matrix-Free Higher-Order Finite Element Methods for Hyperelasticity

TL;DR

This paper introduces a matrix-free finite element solver with an $hp$-multigrid preconditioner for hyperelasticity, achieving significant speed-ups by reducing memory traffic and optimizing evaluation strategies for biomechanics applications.

Contribution

It develops a matrix-free approach for finite-strain elasticity with tailored storage and evaluation strategies, improving computational efficiency on modern hardware.

Findings

Significant speed-ups over traditional matrix-based methods.

Effective strategies for balancing compute load and memory access.

Enhanced performance for higher polynomial degrees in biomechanics models.

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.