Search

The paper investigates matrix-free high-order implementation of finite element discretization with p-multigrid preconditioning for the compressible Neo-Hookean hyperelasticity problem at finite strain on unstructured 3D meshes in parallel. We consider two formulations for the matrix-free action of the Jacobian in Neo-Hookean hyperelasticity: (i) working in the reference configuration to define the second Piola-Kirchhoff tensor as a function of the Green-Lagrange strain S(E) (or equivalently, the right Cauchy-Green tensor C = I+2E), and (ii) working in the current configuration to define the Kirchhoff stress in terms of the left Cauchy-Green tensor τ(b). The proposed efficient algorithm utilizes the Portable, Extensible Toolkit for Scientific Computation (PETSc), along with the libCEED library for efficient compiler optimized tensor-product-basis computation to demonstrate an efficient nonlinear solution algorithm. We utilize p-multigrid preconditioning on the high-order problem with algebraic multigrid (AMG) on the assembled linear Q1 coarse grid operator. In contrast to classical geometric multigrid, also known as h-multigrid, each level in p-multigrid is related to a different approximation polynomial order p, instead of the element size h. A Chebyshev polynomial smoother is used on each multigrid level. AMG is then applied to the assembled Q1 (trilinear hexahedral elements), which allows low storage that can be efficiently used to accelerate convergence to a solution. For the compressible Neo-Hookean hyperelastic constitutive model we exploit the stored energy density function to compute the stored elastic energy density of the Neo-Hookean material as it relates to the deformation gradient. Based on our formulation, we consider four different algorithms each with different storage strategies. Algorithms 1 and 3 are implemented in the reference and current configurations respectively and store ∇Xξ, det(∇ξX), and ∇Xu. Algorithm 2 in the reference configuration stores, ∇Xξ, det(∇ξX), ∇Xu, C−1, and λ log (J). Algorithm 4, in the current configuration, stores det(∇ξX), ∇xξ, τ, and μ – λ log(J). x refers to the current coordinates, X to the reference coordinates, and ξ to the natural coordinates. We perform 3D bending simulations of a tube composed of aluminum (modulus of elasticity E = 69 GPa, Poisson’s ratio v = 0.3) using unstructured meshes and polynomials of order p = 1 through p = 4 under mesh refinement. We explore accuracy-time-cost tradeoffs for the prediction of strain energy across the range of polynomial degrees and Jacobian representations. In all cases, Algorithm 4 using the current configuration formulation outperforms the other three algorithms and requires less storage. Similar simulations for large deformation compressible Neo-Hookean hyperelasticity as applied to the same aluminum material are conducted with ABAQUS, a commercial finite element software package which is a state-of-the-art engineering software package for finite element simulations involving large deformation. The best results from the proposed implementations and ABAQUS are compared in the case of p = 2 on an Intel system with @2.4 GHz and 128 GB RAM. Algorithm 4 outperforms ABAQUS for polynomial degree p = 2.

We examine a residual and matrix-free Jacobian formulation of compressible and nearly incompressible (v → 0.5) displacement-only linear isotropic elasticity with high-order hexahedral finite elements. A matrix-free p-multigrid method is combined with algebraic multigrid on the assembled sparse coarse grid matrix to provide an effective preconditioner. The software is verified with the method of manufactured solutions. We explore convergence to a predetermined L2 error of 10-4, 10-5 and 10-6 for the compressible case and 10-4, 10-5 for the nearly-incompressible cases, as the Poisson's ratio approaches 0.5, based upon grid resolution and polynomial order. We compare our results against results obtained from C3D20H mixed/hybrid element available in the commercial finite element software ABAQUS that is quadratic in displacement and linear in pressure. We determine, for the same problem size, that our matrix-free approach for displacement-only implementation is faster and more efficient for quadratic elements compared to the C3D20H element from ABAQUS that is specially designed to handle nearly-incompressible and incompressible elasticity problems. However, as we approach the near incompressibility limit, the number of Conjugate Gradient iterations required to achieve the desired solution increases significantly. &nbsp

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 such as water, 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 the Finite Element Method (FEM), 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 FEM 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 FEM remains a challenge. Thus, this paper addresses efficiency from the perspective of computer science algorithms and programming. The proposed efficient algorithm utilizes the Portable, Extensible Toolkit for Scientific Computation (PETSc), 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 method is presented whereby a hierarchy of levels is constructed. In contrast to classical geometric multigrid, also known as h-multigrid, each level in p-multigrid is related to a different approximation polynomial order, p, instead of the element size, h. A Chebyshev polynomial smoother is used on each multigrid level. Algebraic MultiGrid (AMG) is then applied to the assembled Q1 (linear) coarse mesh on the nodes of the quadratic Q2 (quadratic) mesh. This allows low storage that can be efficiently used to accelerate the convergence to solution. For a Neo-Hookean hyperelastic problem, we examine a residual and matrix-free Jacobian formulation of a tri-quadratic hexahedral finite element with enhancement. Efficiency estimates on AVX-2 architecture 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 FEM software program ABAQUS. The particular problem in consideration is the simulation of an assistive device in the form of finger-bending in 3D.