For the numerical implementation of elastoplastic models, many researchers prefer an implicit integral strategy due to its unconditional stability. There are two computational difficulties for solving implicit integral equations for the elastoplastic model, i.e., the non-smoothness from inequality constraints and the convergence of the solution of nonlinear equations. The most popular implicit algorithm may be the return-mapping algorithm, where the inequality constraints are addressed by the loading/unloading estimation and the nonlinear equations are solved by the Newton method. However, it is found that the iterations may not converge for the Newton method when the initial value is far from the final solution. The loading/unloading estimation also makes the stress update procedure more cumbersome. A more concise and robust computational paradigm is eagerly desired for thriving elastoplastic models of geomaterials. This chapter mainly focuses on addressing the non-smoothness in elastoplastic computations. More robust numerical algorithms for solving nonlinear equations will be studied in Chap. 8. First, this chapter presents the standard implicit stress integration format for elastoplastic models. Then, the classic return mapping stress update strategy is reviewed, which comprises elastic prediction and plastic correction. In Chap. 7.3, we introduce a new technique for dealing with non-smoothness in elastoplastic problems, namely the unconstrained stress update strategy, where the loading/unloading inequality is equivalently replaced by a smooth function. Subsequently, a comparative analysis of two stress update strategies is conducted through the examination of a 1D elastoplastic problem. In the remaining content, the NOEP model in Sec. 3.2 was implemented in ABAQUS based on the unconstrained stress update algorithm, and the performance of the algorithm was evaluated based on several boundary value problems.