Zero average in solution to Laplace problem with Lagrange multipliers 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?
 A: What follows is an explanation how I would impose the mean value zero condition on the discrete level. Afterwards I discuss how this reflects to the continuous weak form.
Let $\varphi_i$, $i\in\{1,\cdots,N\}$, be the finite element basis functions.
Let $A \in \mathbb{R}^{N \times N}$ be the stiffness matrix of the original problem with the components $A_{ij}=\int_\Omega \nabla \varphi_j\cdot \nabla \varphi_i\,\mathrm{d}x$ and $f \in \mathbb{R}^N$ be the load vector of the original problem with the components $f_j = \int_\Omega f \varphi_i \,\mathrm{d}×$. Let us further define the (column) vector $g  \in \mathbb{R}^N$ with the components $g_j = \int_\Omega \varphi_j\,\mathrm{d}x$. Then the mean value zero is imposed by solving the modified matrix system:
$$
\begin{bmatrix} A & g \\
g^T & 0
\end{bmatrix}
\begin{bmatrix}
u \\
\lambda
\end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix},
$$
where $u \in \mathbb{R}^N$ is the solution vector and $\lambda \in \mathbb{R}$ is the Lagrange multiplier.
Notice that the second equation is simply $g^Tu=0$ (= mean value of $u$ is zero).
The corresponding continuous and constrained minimization problem reads:
$$
\min_{u \in H^1_0(\Omega),\\\int_\Omega u\,\mathrm{d}x=0} \frac12 \int_\Omega |\nabla u|^2\,\mathrm{d}x - \int_\Omega fv\,\mathrm{d}x.$$
Using a Lagrange multiplier we may unconstrain this by searching for the saddle-point of
$$
\frac12 \int_\Omega |\nabla u|^2\,\mathrm{d}x - \int_\Omega fv\,\mathrm{d}x+\lambda \int_\Omega u\,\mathrm{d}x,$$
where $\lambda \in \mathbb{R}$ is the Lagrange multiplier. Computing the derivative of this expression with respect to $u$ and $\lambda$ gives the two equations:
$$
\begin{aligned}
\int_\Omega \nabla u \cdot \nabla v \,\mathrm{d}x+\lambda \int_\Omega v\,\mathrm{d}x&=\int_\Omega fv\,\mathrm{d}x, \quad \forall v \in H^1_0(\Omega),\\
\int_\Omega u \,\mathrm{d}x &= 0.
\end{aligned}$$
What we wrote before is essentially the discretization of this system.
