We present the analysis of an [Formula presented] version Finite Element Method for the approximation of the solution to singularly perturbed reaction–diffusion problems posed in smooth domains [Formula presented]. The method uses piecewise polynomials of degree [Formula presented] in each variable, defined on an exponentially graded mesh, optimally constructed for the approximation of exponential layers. We establish robust, optimal convergence rates in a variety of norms and illustrate our theoretical findings through numerical computations.