We consider a linear singularly perturbed reaction-diffusion problem in three dimensions and its numerical solution by a Galerkin finite element method with trilinear elements. The problem is discretised on a Shishkin mesh with $N$ intervals in each coordinate direction. Derivation of an error estimate for such a method is usually based on the (Shishkin) decomposition of the solution into distinct layer components. Our contribution is to provide a careful and detailed analysis of the trilinear interpolants of these components. From this analysis it is shown that, in the usual energy norm the errors converge at a rate of $\mathcal{O}(N^{-2}+\varepsilon^{1/2}N^{-1}\ln N)$. This is validated by numerical results.