Abstract: The mathematical programming approach of elastoplastic incremental analysis is one of the effective ways to analyze deformation and strength problems in geotechnical engineering and has unique advantages in dealing with the complex problems such as non-smooth yield surface, contact conditions and multi-surface plasticity. To further simplify the computational framework and overcome the volumetric locking, a novel mixed constant stress-smoothed strain three-node triangle element with a two-level mesh repartitioning scheme is proposed to discretize the generalized Hellinger-Reissner (GHR) variational principle, the boundary value problem of elastoplastic problem can be reformulated as a conic programming problem under the constraint of the associated flow rule, and the cohesive-frictional contact condition is treated as a set of conic constraints and introduced into the conic programming problem of the elastoplastic incremental analysis. Then, an efficient primal-dual interior point algorithm is used to solve it. Finally, the proposed method is applied to two classical geotechnical engineering problems. The results show that the new method is superior to the traditional mixed six-node triangular element in terms of the computational efficiency, convergence and accuracy.