I am trying to solve the Helmholtz PDE $-\nabla^2u + u = f \quad \text{on} \; \Omega$ with homogeneous Neumann boundary condition $\partial_n u = 0 \; \text{on} \; \partial \Omega$.
I noted that an unique solution does exist given the PDE contains the "u" term as itself. I proceeded to discretize the domain using linear quad elements and further use Galerkin method to arrive at the stiffness expression $[A]\{u\} = \{F\}$, where the elemental stiffness matrix was $[A_e]_{ij} = \int \limits_{\Omega_e} (\phi_i\phi_j + \nabla \phi_i \nabla \phi_j) d\Omega_e$.
However, I find that the matrix are singular and I am unable to decipher why they are so. Given that the BC is homogeneous Neumann everywhere, I suppose no special consideration be given to include them ( Am I wrong to assume so ? ).
Thanks for the help!