Tutorial/Laplace equation with Neumann boundary condition

We propose here to solve a first very simple academic example with GMSH and GetDP. The problem is the following: \begin{equation} \begin{cases}\label{eq:problemU} \Delta u + u = f & \text{in } \Omega\\ \displaystyle{\frac{\partial u}{\partial \mathbf{n}} = 0} & \text{on }\partial\Omega, \end{cases} \end{equation} where $\Omega$ is the unit square, $\partial\Omega$ its boundary and $f$ the two-variables function defined by $$ \forall x,y\in [0,1]^2,\qquad f(x,y) = (1+2\pi^2)\cos(\pi x)\cos(\pi y) $$ Then, the unique solution of problem \eqref{eq:problemU} is $$ \forall x,y\in[0,1]^2, \qquad u(x,y) = \cos(\pi x)\cos(\pi y). $$ The variationnal formulation of problem \eqref{eq:problemU} reads $$ \left\{\begin{array}{l} \text{Find } u\in H^2(\Omega) \text{ such that, }\\ \displaystyle{\forall v\in H^1(\Omega), \qquad \int_{\Omega} \nabla u\cdot\nabla v \;d\Omega + \int_{\Omega}uv \;d\Omega - \int_{\Omega}fv\;d\Omega = 0}. \end{array}\right. $$