Mathematical modeling of real-world phenomena often leads to formulations of problems in infinite dimensions containing variational inequalities. In this talk, we will focus on so-called elasto-plastic problems which require the pointwise bounding of the gradient of the displacement, i.e., the stress on a body at each point in the presence of a given force. We will then derive the first order optimality conditions. To use the semismooth Newton method to solve the problem numerically, we require a regularized penalization of our problem. Numerical path-following strategies will be developed and analytical results can be numerically verified.