Analytical Mohr–Coulomb return mapping and consistent tangent operator: a finite element implementation in Mathematica

An analytical solution for the projected stresses is presented for the perfectly plastic Mohr–Coulomb criterion within the context of elastoplastic initial-value problem integration. The solution is obtained via the closest-point projection method in a rotated Haigh–Westergaard space (HWR), where a distance function is defined in terms of the trial stresses using an energy norm. The projected stresses are expressed in closed form as functions of the trial stresses. The methodology readily extends to other classical associative plasticity models, such as Drucker–Prager and von Mises. The derivation of the fully consistent tangent operator is also presented in a form applicable to general yield criteria. Its computation requires only the projected principal stresses and the principal directions of the trial stresses. A dedicated regularization is introduced for repeated principal values, typically encountered in edge-return projections for Mohr–Coulomb, without resorting to smoothing or active-set logic. This yields a continuous and fully consistent tangent across faces, edges, and the apex. The framework is demonstrated for the associative Mohr–Coulomb model, showing robustness, accuracy, and straightforward implementation in standard finite-element codes.