Let us consider the problem: find $u:\Omega \rightarrow \mathbb{R}$ such that
\begin{equation} \begin{cases} \Delta u = -f \quad &\text{on } \Omega,\\ \nabla u \cdot \textbf{n} = 0 \quad &\text{on } \partial \Omega, \end{cases} \end{equation}
where $\Omega$ is the domain of the problem, $\partial \Omega$ is the boundary and $f:\Omega \rightarrow \mathbb{R}$. Let us also consider the condition $\int_{\Omega} f d \textbf x= 0$, which is a necessary compatibility condition.
This problem is not well-posed. In fact, if $\tilde u$ is solution of the problem, also $\tilde u + C$ for any $C \in \mathbb{R}$ is. By adding the additional constraint of null average on the domain, i.e. \begin{equation} \dfrac{1}{|\Omega|}\int_{\Omega} u d \textbf x= 0, \end{equation} the problem has a unique solution.
Here is the question. Let us suppose that I want to solve the problem using FEM. I know that the imposition of the null average on the solution of the problem can be imposed by using Lagrange multipliers, but I do not understand if I should introduce additional terms in the weak formulation to do so. How can I use Lagrange multipliers to impose the constraint of zero average to this problem?