Thank you for your reply. I am sorry, I have sent you the wrong variational problem. Here is the variational functional which I think should implement the constraint:
F = F_{\rm Poisson} + F_{\rm constraint}
where F_{\rm Poisson} is the functional related to the Poisson equation that you wrote:
F_{\rm Poisson}[u, v] = \int_\Omega (\frac{1}{2}[ \nabla (u+v) ]\cdot[ \nabla (u+v)] + f (u+v)) ~ \mathrm{d}x.
If I vary F_{\rm Poisson} with respect to u and v I obtain
\delta F_{\rm Poisson} = \int_{\Omega} \{[ \nabla (u+v) ]\cdot[ \nabla (\nu_u+\nu_v)] ] + f (\nu_u + \nu_v)\} = \int_{\Omega} \{[ -\nabla^2 (u+v) +f] (\nu_u + \nu_v)\}, where \nu_u = \delta u, \nu_v = \delta v are the test functions related to u and v, respectively.
\delta F_{\rm Poisson} =0 is the equation which I obtain if I consider the original PDE
\nabla^2 (u+v)-f=0, multiplied it by \nu_u + \nu_v, and integrated.
Following the example in your post, the functional F_{\rm constraint} which enforces the constraint u^2-v^2=g could be (?)
F_{\rm constraint}[u,v] = \frac{\alpha}{2} \int_{\Omega} dx (u^2 -v^2-g)^2.
I don’t know why in your post you wrote the integral over the boundary only. Here F_{\rm constraint} is an integral over \Omega, dx, because I want to enforce the constraint in \Omega. If I vary F_{\rm constraint} I get
\delta F_{\rm constraint} = \alpha \int_{\Omega} dx (u^2 -v^2-g) (2 u \nu_u - 2 v \nu_v).
Would I obtain the solution of the PDE above with the constraint above if I solved the variational problem \delta F_{\rm Poisson} + \delta F_{\rm constraint} = 0?