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}+$ε$^{1/2}$N$^{−1}ln$N$). This is validated by numerical results.