The pure traction (Neumann) problem is classical in the theory of linear elasticity but its numerical solution is delicate because its solution is unique in a quotient space, resulting in a functional setting that is not easily ameanable to Galerkin-type approximations. In standard finite element implementations, uniqueness is usually enforced either by adding carefully chosen essential boundary conditions at the discrete level or by imposing global constraints on the deformation. In practice, these fixes are often awkward to set up and can be computationally inefficient. Here we introduce an alternative strategy that avoids these drawbacks. The proposed method starts from a regularized version of the traction problem whose solution is provably unique. We show that this solution converges to the minimal-norm solution of the original Neumann problem. Crucially, the regularized problem can be discretized with conventional finite elements directly—without introducing extra unknowns or special constraints. We also address a second, common difficulty: when the computational mesh approximates the geometry, the discrete loads may fail to satisfy equilibrium, which can make the discrete pure-traction problem ill-posed. To deal with this, we propose a regularized predictor–corrector finite element formulation that systematically accommodates load incompatibilities. We prove that the resulting approximation converges to the solution of the original Neumann problem as both the mesh size and the regularization parameter go to zero. Representative numerical tests demonstrate the robustness and accuracy of the approach in typical mechanics problems governed by pure traction boundary conditions.