Suppose $$\mathcal{L} =\mathcal{L}(x,y,u,u_x,u_y) = \frac{1}{2} \lVert \nabla u \rVert^2$$
and I want to find $u$ such that the functional
$$ E(u)=\int_{\Omega} \mathcal{L}dxdy $$
is minimized, computing the variational derivative leads me to
$$ \frac{\delta E}{\delta u} = \int_{\partial \Omega}(u_y - u_x)h d\gamma -\int_{\Omega} \Delta u h dxdy $$
therefore the $u$ function I'm looking for is given by
$$ \left\{ \begin{array}{ll} \left(\frac{\partial}{\partial y} - \frac{\partial}{\partial x}\right) u = 0 & u (x,y) \in \partial \Omega \\ \Delta u = 0 & (x,y) \in \text{Int}(\Omega) \end{array} \right. $$
Though the question might be silly, I wonder how the operator
$$ A =\frac{\partial}{\partial y} - \frac{\partial}{\partial x} $$
might be discretized, my attempt would be using be applying the Sobel masks in pixels at $\partial \Omega$ and adding the results, would this be correct? Or is that operator something know that I'm missing? I thought at the beginning it was the divergence operator, but I was wrong since $u$ is a real valued function.
Update (Full derivation of the integral)
We have $$ \frac{E(u + \alpha h) - E(u)}{\alpha} = \frac{1}{\alpha}\int_{\Omega} \mathcal{L}(x,y,u+\alpha h,u_x+\alpha h_x, u_y+\alpha h_y) -\mathcal{L}(x,y,u,u_x, u_y)dxdy = \\ \int_{\Omega} \frac{\partial \mathcal{L}}{\partial u}h + \frac{\partial \mathcal{L}}{\partial u_x}h_x + \frac{\partial \mathcal{L}}{\partial u_y}h_y dxdy = \int_{\Omega} \frac{\partial \mathcal{L}}{\partial u}h dxdy + \int_{\Omega}\frac{\partial \mathcal{L}}{\partial u_x}h_x + \frac{\partial \mathcal{L}}{\partial u_y}h_y dxdy. $$
Consider the integral
$$ \int_{\Omega}\frac{\partial \mathcal{L}}{\partial u_x}h_x + \frac{\partial \mathcal{L}}{\partial u_y}h_y dxdy $$
observing
$$ \frac{\partial \mathcal{L}}{\partial u_x}h_x = \frac{d}{dx}\left(\frac{\partial \mathcal{L}}{\partial u_x}h\right) - \frac{d}{dx} \frac{\partial \mathcal{L}}{\partial u_x} h $$
and likewise for $\frac{\partial \mathcal{L}}{\partial u_y}h_y$ we have $$ \int_{\Omega}\frac{\partial \mathcal{L}}{\partial u_x}h_x + \frac{\partial \mathcal{L}}{\partial u_y}h_y dxdy = \int_{\partial \Omega} \left(\frac{\partial \mathcal{L}}{\partial u_y} - \frac{\partial \mathcal{L}}{\partial u_x}\right)h d\gamma - \int_{\Omega} \left( \frac{d}{dx} \frac{\partial \mathcal{L}}{\partial u_x} + \frac{d}{dy} \frac{\partial \mathcal{L}}{\partial u_y}\right)h dxdy $$
Taking the limit for $\alpha \to 0$ yields
$$ \frac{\delta E}{\delta u} = \int_{\Omega} \left(\frac{\partial \mathcal{L}}{\partial u} - \frac{d}{dx} \frac{\partial \mathcal{L}}{\partial u_x} - \frac{d}{dy} \frac{\partial \mathcal{L}}{\partial u_y}\right)h dxdy + \int_{\partial \Omega} \left(\frac{\partial \mathcal{L}}{\partial u_y} - \frac{\partial \mathcal{L}}{\partial u_x}\right)h d\gamma $$
And explicit computation of the partial derivatives of $\mathcal{L}$ yields to expression above.
The variation should be: $$ \delta\int_\Omega \frac12|\boldsymbol{\nabla} u|^2 = -\int_\Omega\Delta u\ \delta u + \int_{\partial\Omega}\boldsymbol{\nabla} u\cdot\boldsymbol n\ \delta u $$
So the second equation is $$ \left.\boldsymbol{n}\cdot\boldsymbol{\nabla} u = 0 \right|_{\partial\Omega}, $$ where $\boldsymbol n$ is a normal vector to $\partial\Omega$. Or in other words that gradient flow $\boldsymbol{\nabla} u$ should be perpendicular to the boundary $\partial\Omega$.
Discretization of this depends on many things. If boundary is simple (rectangular), you can just set one of the sobel values to zero. If the boundary is complex, you might try to detect pixels of boundaries, calculate the tangent there and write linear equations on sobel values.
Edit
Consider your integral $$ I = \int_\Omega\left(\frac{\partial\mathcal L}{\partial u_x}h_x+\frac{\partial\mathcal L}{\partial u_y}h_y\right)dx\,dy = \int_\Omega (F\,h_x+G\,h_y)\,dx\,dy =\\ \int_\Omega\left(\frac\partial{\partial x}(F h)-F_xh+\frac\partial{\partial y}(G h)-G_yh\right)dx\,dy = \\ -\int_\Omega h\left(F_x+G_y\right)dx\,dy + \int_\Omega\left(\frac\partial{\partial x}(F h)+\frac\partial{\partial y}(G h) \right)dx\,dy. $$
Green's theorem for second integral gives $$ \int_{\partial\Omega}-Ghdx+Fhdy = \int_{\partial\Omega}h(F\,dn_x + G\,dn_y). $$ where $\boldsymbol{dn} = (dy, -dx)$