  Pointwise error analysis of the linear finite element approximation for
$-\Delta u + u = f$ in $\Omega$, $\partial_n u = \tau$ on $\partial\Omega$,
where $\Omega$ is a bounded smooth domain in $\mathbb R^N$, is presented. We
establish $O(h^2|\log h|)$ and $O(h)$ error bounds in the $L^\infty$- and
$W^{1,\infty}$-norms respectively, by adopting the technique of regularized
Green's functions combined with local $H^1$- and $L^2$-estimates in dyadic
annuli. Since the computational domain $\Omega_h$ is only polyhedral, one has
to take into account non-conformity of the approximation caused by the
discrepancy $\Omega_h \neq \Omega$. In particular, the so-called Galerkin
orthogonality relation, utilized three times in the proof, does not exactly
hold and involves domain perturbation terms (or boundary-skin terms), which
need to be addressed carefully. A numerical example is provided to confirm the
theoretical result.
