Confusion with Neumann and Dirichlet

Question

If I solve Laplace's equation with Neumann boundary conditions then everything is defined via derivatives. Consequently one needs to fix a point with a specific value to get a solution. However if I fix my value with a Dirichlet condition the solution is distorted. Where am I going wrong?

Edit to question I think I have been asking Mathematica to solve an impossible problem. See my answer below.

Continue with original question

Here is an example

 Needs["NDSolve`FEM`"];
 x2 = 4; y2 = 1;

 reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];

 mesh = ToElementMesh[reg,
   "BoundaryMeshGenerator" -> {"Continuation"},
   MaxCellMeasure -> .002,
   "MaxBoundaryCellMeasure" -> 0.01];

 Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]

Here is a Laplacian with Neumann boundary conditions so that everything is defined through derivatives.

sol = NDSolveValue[{
    Laplacian[u[x, y], {x, y}] == 
     NeumannValue[1, 0 <= y <= y2 && x == 0] + 
      NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2]
    }, u, {x, y} \[Element] mesh];

This gives a message along the expected lines

NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {u}; the result is not unique up to a constant. >>

Now we repeat with a Dirichlet condition that gives a value in the corner

 sol = NDSolveValue[{
        Laplacian[u[x, y], {x, y}] == 
         NeumannValue[1, 0 <= y <= y2 && x == 0] + 
          NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2],
        DirichletCondition[u[x, y] == 0, x == x2 && y == 0]
        }, u, {x, y} \[Element] mesh];
    Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]

The value in the corner is as expected but the whole solution is distorted to reach the value.

If I change the location of the Dirichlet point then the location of the distortion changes as perhaps might be expected.

sol = NDSolveValue[{
    Laplacian[u[x, y], {x, y}] == 
     NeumannValue[1, 0 <= y <= y2 && x == 0] + 
      NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2],
    DirichletCondition[u[x, y] == 0, x == x2 && y == y2]
    }, u, {x, y} \[Element] mesh];
Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]

If I look at the Laplacian of the solution then the point continues to appear

 ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[sol[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}, 
 PlotRange -> All]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]

The values in the above plots should be zero and they are small except where I have my Dirichlet point.

So what do I not understand? Is it something to do with Neumann values being associated with the normal to the boundary and this not being compatible with a defined point? Please enlighten me.

Answer 1

Following a comment from Young it may be that I have set up an impossible problem. As the equations are linear then the principle of superposition applies and we can look at the solutions from each Neumann condition separately. Starting with the second Neumann condition on its own, with the Dirichlet condition we have.

Needs["NDSolve`FEM`"];
x2 = 4; y2 = 1;
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg,
   "BoundaryMeshGenerator" -> {"Continuation"},
   MaxCellMeasure -> .002,
   "MaxBoundaryCellMeasure" -> 0.01];
sol = NDSolveValue[{
    Laplacian[u[x, y], {x, y}] == 
     NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2],
    DirichletCondition[u[x, y] == 0, x == x2 && y == 0]
    }, u, {x, y} \[Element] mesh];
Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]

There is no distortion at {4,0}. Also we can check the solution by looking at the Laplacian

ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[sol[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}, 
 PlotRange -> All]

The Laplacian is good with small values everywhere.

Now we look at the solution with the first Neumann condition and the Dirichlet condition.

sol = NDSolveValue[{
    Laplacian[u[x, y], {x, y}] == 
     NeumannValue[1, 0 <=  y <= y2 && x == 0],
    DirichletCondition[u[x, y] == 0, x == x2 && y == 0]
    }, u, {x, y} \[Element] mesh];
Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]

The distortion is back. However do the boundary conditions prescribe a possible configuration for Laplace's equation? We have a gradient on one boundary and zero gradient on the other boundaries. From a physical fluid mechanics viewpoint this is like having an inlet on one edge but no outlet. Such a condition is not physically possible for an incompressible fluid. We may illustrate the solution using a StreamPlot.

ClearAll[f];
f[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}, "Cartesian"]];
StreamPlot[f[x, y], {x, 0, x2}, {y, 0, y2}, AspectRatio -> Automatic]

There is an inflow on the left and an outflow through the Dirichlet point. This looks like a physical solution but the distortion around the Dirichlet point has become essential to let the fluid out.

In conclusion it looks like I was trying to solve an impossible problem and Mathematica gave me a possible answer to a slightly different problem. I guess it would be difficult to flag up a warning message that the impossible was being asked for.

Answer 2

Update

As a follow-up to my comment, my hypothesis is that the first Neumann boundary is problematic and a variation of the problem could be solved.

Another issue seems to be associated with the meshed region. I lowered the "MeshOrder" to 1 and the Laplacian of the solution becomes perfectly 0.

Needs["NDSolve`FEM`"];
x2 = 4; y2 = 1;
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg, "BoundaryMeshGenerator" -> {"Continuation"},
    MaxCellMeasure -> .002, "MaxBoundaryCellMeasure" -> 0.01, "MeshOrder" -> 1];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 
      NeumannValue[Cos[2 π x/x2], 0 <= x <= x2 && y == y2],
      DirichletCondition[u[x, y] == 0, x == x2 && y == 0]},
      u, {x, y} ∈ mesh];

Plot3D[sol[x, y], {x, y} ∈ mesh, BoxRatios -> {x2, y2, 1}, PlotRange -> All, 
 ColorFunction -> "Rainbow", Mesh -> None]

ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[sol[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}, PlotRange -> All, 
 ColorFunction -> "Rainbow", Mesh -> None]

Answer 3

If you use Dirichlet boundary condition(BC) you are imposing that your solution must be a known value at that particular point. In your case you don't know the value to specify. Your solution will have a discontinuity unless you specify a correct value. You can minimize the problem refining the mesh at the borders. Take a look at the code below:

   Needs["NDSolve`FEM`"];
x2 = 4; y2 = 1;
(*reg=BoundaryDiscretizeRegion[Rectangle[{0,0},{x2,y2}]];*)
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg, "BoundaryMeshGenerator" -> {"Continuation"},
    MaxCellMeasure -> 0.5, "MaxBoundaryCellMeasure" -> 0.001, 
   "MeshOrder" -> 1];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]
op = Laplacian[u[x, y], {x, y}]
(*Subscript[\[CapitalGamma], D]=DirichletCondition[u[x,y]\[Equal]0, x\
\[Equal]0&&0\[LessEqual]y\[LessEqual]y2]*)
Subscript[\[CapitalGamma], D] = 
 DirichletCondition[u[x, y] == 0, x == x2  && y == 0]
(*Subscript[\[CapitalGamma], N]=NeumannValue[-1,y\[Equal]1 &&0\
\[LessEqual]x\[LessEqual]4]*)
Subscript[\[CapitalGamma], N] = 
 NeumannValue[1, x == 0 && 0 <= y <= y2] + 
  NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2]
uif = NDSolveValue[{op == Subscript[\[CapitalGamma], N], 
   Subscript[\[CapitalGamma], D]}, u, {x, y} \[Element] mesh]
ContourPlot[uif[x, y], {x, y} \[Element] mesh, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]
Plot3D[uif[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]
ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[uif[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}, 
 PlotRange -> All]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]