The Galerkin finite element solutionuhu_hof the Poisson equation−Δu=f-\Delta u=funder the Neumann boundary condition in a possibly nonconvex polygonΩ\varOmega, with a graded mesh locally refined at the corners of the domain, is shown to satisfy the following maximum-norm stability:‖uh‖L∞(Ω)≤Cℓh‖u‖L∞(Ω),\begin{align*} \|u_h\|_{L^{\infty }(\varOmega )} \le C\ell _h\|u\|_{L^{\infty }(\varOmega )} , \end{align*}whereℓh=ln(2+1/h)\ell _h = \ln (2+1/h)for piecewise linear elements andℓh=1\ell _h=1for higher-order elements. As a result of the maximum-norm stability, the following best approximation result holds:‖u−uh‖L∞(Ω)≤Cℓh‖u−Ihu‖L∞(Ω),\begin{align*} \|u-u_h\|_{L^{\infty }(\varOmega )} \le C\ell _h\|u-I_hu\|_{L^{\infty }(\varOmega )} , \end{align*}whereIhI_hdenotes the Lagrange interpolation operator onto the finite element space. For a locally quasi-uniform triangulation sufficiently refined at the corners, the above best approximation property implies the following optimal-order error bound in the maximum norm:‖u−uh‖L∞(Ω)≤{Cℓhhk+2−2p‖f‖Wk,p(Ω)amp;if r≥k+1,Cℓhhk+1‖f‖Hk(Ω)amp;if r=k,\begin{align*} \|u-u_h\|_{L^\infty (\varOmega )} \le \begin {cases} C\ell _h h^{k+2-\frac {2}{p}} \|f\|_{W^{k,p}(\varOmega )} &\text {if $r\ge k+1$}, \\ C\ell _h h^{k+1} \|f\|_{H^{k}(\varOmega )} &\text {if $r=k$}, \end{cases} \end{align*}wherer≥1r\ge 1is the degree of finite elements,kkis any nonnegative integer no larger thanrr, andp∈[2,∞)p\in [2,\infty )can be arbitrarily large.