Abstract: Due to the large number of elements required in the calculation of solid finite element models, large computational resource is consumed in the element state determination process and the factorization of the global tangent stiffness matrix with large dimensions, thus resulting in low efficiency. In this paper, linear tetrahedron and hexahedron isoparametric elements are established based on the inelasticity-separated finite element method. Six direct integration points are considered for hexahedron elements as the nonlinear strain interpolation points instead of eight gauss integration points, of which the computational accuracy is stable and the efficiency is improved. In addition, the main solving process of the governing equation is only the back substitution of the initial stiffness matrix and the matrix-vector multiplication by using the Woodbury formula and combined approximation approach. Therefore, the efficiency is significantly improved. Finally, the computational efficiency of the proposed method based on the time complexity theory indicates that, with the increase of the number of nodal degrees of freedom, the computational efficiency of the proposed method is significantly improved as compared with the traditional variable stiffness method. The numerical examples verify the correctness of the proposed solid element model and the high efficiency of the proposed algorithm.