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.