Continuing this topic, once the Dini-Neumann problem is solved:
$$ \begin{cases} \Delta f = 0 & \text{on} \; \Omega \\ \frac{\text{d}}{\text{d} n} f = y\,n_x - x\,n_y & \text{on} \; \partial\Omega \end{cases} $$
with $n = (n_x,\,n_y)$ exterior normal vector to $\partial\Omega$, the goal is the calculation of:
$$ \color{blue}{J \equiv \iint_{\Omega} \left(f_x\,y - f_y\,x\right)\text{d}x\,\text{d}y}\,. $$
On the other hand, we also have:
$$ J = \iint_{\Omega} \left[\left(f\,y\right)_x - \left(f\,x\right)_y\right]\text{d}x\,\text{d}y $$
which by the divergence theorem is equivalent to:
$$ J = \int_{\partial\Omega} f\left(y\,n_x - x\,n_y\right)\text{d}s $$
and for the boundary condition:
$$ J = \int_{\partial\Omega} f\left(\frac{\text{d}}{\text{d}n}f\right)\text{d}s = \int_{\partial\Omega} \left(f\,f_x\,n_x + f\,f_y\,n_y\right)\text{d}s\,. $$
This done, again by the divergence theorem:
$$ J = \iint_{\Omega} \left[(f\,f_x)_x + (f\,f_y)_y\right]\text{d}x\,\text{d}y = \iint_{\Omega} \left(f_x^2 + f_y^2 + f\,\Delta f\right)\text{d}x\,\text{d}y $$
that is, since $f$ is an harmonic function on $\Omega$:
$$ \color{\red}{J = \iint_{\Omega} \left(f_x^2 + f_y^2\right)\text{d}x\,\text{d}y}\,, $$
known in the literature as the Dirichlet integral.
Therefore, after determining the harmonic function:
Needs["NDSolve`FEM`"];
Ω = ToElementMesh[ImplicitRegion[x^4 + y^4 <= 1, {x, y}], "IncludePoints" -> {{0, 0}}];
f = NDSolveValue[{Laplacian[g[x, y], {x, y}] == NeumannValue[y (4 x^3) - x (4 y^3), True],
DirichletCondition[g[x, y] == 0, x == 0 && y == 0]}, g, {x, y} ∈ Ω];
and verified that it's accurate enough:
NIntegrate[f[x, y], {x, y} ∈ Ω]
NIntegrate[Laplacian[f[x, y], {x, y}], {x, y} ∈ Ω]
-0.00000293157
0.00292529
I tried to calculate the Dirichlet integral in its two equivalent versions:
NIntegrate[D[f[x, y], x] y - D[f[x, y], y] x, {x, y} ∈ Ω]
NIntegrate[D[f[x, y], x]^2 + D[f[x, y], y]^2, {x, y} ∈ Ω]
-0.48313
1.84564
but it's evident that something went wrong, but I don't know what! Ideas? Thank you!
EDIT: Thanks to the precious comment of xzczd I understood that:
$(n_x,n_y)$ must necessarily have unitary length;
before the Laplace operator it must absolutely be a minus sign.
In fact, writing:
c[x_, y_] = x^4 + y^4;
{nx, ny} = Normalize[Grad[c[x, y], {x, y}]];
Needs["NDSolve`FEM`"];
Ω = ToElementMesh[ImplicitRegion[c[x, y] <= 1, {x, y}],
"IncludePoints" -> {{0, 0}}];
f = NDSolveValue[{-Laplacian[g[x, y], {x, y}] == NeumannValue[y nx - x ny,
True], DirichletCondition[g[x, y] == 0, x == 0 && y == 0]},
g, {x, y} ∈ Ω];
NIntegrate[f[x, y], {x, y} ∈ Ω]
NIntegrate[Laplacian[f[x, y], {x, y}], {x, y} ∈ Ω]
NIntegrate[D[f[x, y], x] y - D[f[x, y], y] x, {x, y} ∈ Ω]
NIntegrate[D[f[x, y], x]^2 + D[f[x, y], y]^2, {x, y} ∈ Ω]
we get:
0.000000504282
-0.000799624
0.126589
0.126595
which is approximately what I want. The only thing that would remain to understand is whether it were possible to obtain a better approximation of those integrals. Thanks again!
