<p>The Mohr–Coulomb criterion represents a hexagonal pyramidal surface in the three-dimensional principal stress/strain spaces. There exist numerical singularities on the edges and at the apex of the pyramidal surface. An ideal technique for removing such singularities should be simple enough and can minimize the deviation. An early approach seems to meet the requirements but was found to be not smooth and poor in numerical astringency. Therefore, a simple piecewise C2 continuous smoothing technique was presented to remove the singularities posed by the vertexes of the Mohr–Coulomb locus in the deviatoric plane. Meanwhile, the linear Mohr–Coulomb envelope in the meridian plane was still preserved. The singularity at the apex was also removed. A subroutine that still uses the original Mohr–Coulomb criterion to determine the yielding of geomaterials but employs the piecewise C2 continuous smooth yield and potential surfaces to perform the relevant derivation operations was successfully developed in Abaqus using the fully implicit backward Euler integral regression algorithm. Its reliability was extensively validated by numerically modeling the triaxial compression and uniaxial tensile tests, as well as the stability analysis of a typical soil slope. The study highlighted the merit of this technique to be popularized in relevant numerical software.</p>