Efficient Residual and Matrix-Free Jacobian Evaluation for Three-Dimensional Tri-Quadratic Hexahedral Finite Elements With Nearly-Incompressible Neo-Hookean

Soft materials such as rubber, elastomers, and soft biological tissues mechanically deform at large strain isochorically for all time, or during their initial transient (when a pore fluid, typically incompressible, does not have time to flow out of the deforming polymer or soft tissue porous skeleton).  Simulating these large isochoric deformations computationally, such as with Finite Element Analysis (FEA), requires higher order (typically quadratic) interpolation functions and/or enhancements through hybrid/mixed methods to maintain stability.  Lower order (linear) finite elements with hybrid/mixed formulation may not perform stably for all mechanical loading scenarios involving large isochoric deformations, whereas quadratic finite elements with or without hybrid/mixed formulation typically perform stably, especially when large bending or folding deformations are being simulated.  For topology-optimization design of soft robotics, for instance, the FEA solid mechanics solver must run efficiently and stably.  Stability is ensured by the higher order finite element formulation (with possible enhancement), but efficiency for higher order FEA remains a concern.  Thus, this paper addresses the efficiency concern from the perspective of computer science algorithms and programming. The proposed efficient algorithm utilizes the parallel computing library PETSc for fast parallel computation, along with the libCEED library for efficient compiler optimized tensor-product-basis computation to demonstrate an efficient nonlinear solution algorithm. For preconditioning, a scalable p-multigrid discontinuous Galerkin method is presented whereby a hierarchy of levels is constructed. Each level represents a different approximation polynomial order $p$ instead of the mesh width $h$ at each iteration of the nonlinear Newton solver.  The effect of different relaxation schemes such as block Jacobi, Chebychev smoothers, and Vanka smoothers in form of block Gauss-Seidel are considered. For a Neo-Hookean hyperelastic problem,  we examine a residual and matrix-free Jacobian formulation of a tri-quadratic hexahedral finite element with enhancement. The approximation orders of $p=4$ and $p=2$ are considered in the p-multigrid preconditioning. The Algebraic MultiGrid (AMG) linear equation solver is then applied to the assembled $Q_1$ (linear) coarse element on the nodes of the quadratic $Q_2$ (quadratic) mesh for the unstructured mesh.  This allows low storage that can be efficiently used to accelerate the convergence to  solution. Efficiency estimates on different architectures such as AVX-2 and AVX-512 based on CPU time are provided as a comparison to similar simulation (and mesh) of isochoric large deformation hyperelasticity as applied to soft materials conducted with the commercially-available FEA software program ABAQUS. The particular problem in consideration is the simulation of an assistive device in the form of finger-bending in 3D.