Piecewise imposes internal boundaries in NDSolve - is this expected?

Question

In the following code I used True as the predicate for DirichletCondition and found that the boundary condition was applied not only at the boundary of the specified region but also at the internal boundary implied by the Piecewise potential.

f[x_, y_] := Piecewise[{{-1, x^2 + y^2 < 0.25}}, 0]

sol = NDSolveValue[{
    Laplacian[u[x, y], {x, y}] == f[x, y],
    DirichletCondition[u[x, y] == 0, True]
    }, u, {x, y} ∈ Disk[]];

Plot3D[sol[x, y], {x, y} ∈ Disk[], PlotPoints -> 50]

I can confirm that there is an internal boundary by looking at the boundary of the FEM mesh:

NDSolve`FEM`ToBoundaryMesh[sol["ElementMesh"]]["Wireframe"]

The result I wanted can be obtained by explicitly giving just the outer boundary to DirichletCondition:

sol = NDSolveValue[{
    Laplacian[u[x, y], {x, y}] == f[x, y],
    DirichletCondition[u[x, y] == 0, {x, y} ∈ Circle[]]
    }, u, {x, y} ∈ Disk[]];

Plot3D[sol[x, y], {x, y} ∈ Disk[], PlotPoints -> 50]

My question is whether the behaviour using True is correct, i.e. expected and documented? My expectation was that the only the boundary of the explicitly given Disk[] region would be used (though I can understand why the internal FEM code creates the internal boundary).

Answer 1

Let us look at what happens a bit and start with your first question (title):

"Piecewise imposes internal boundaries in NDSolve - is this expected"

Yes, that is expected if the function introduces a discontinuity in one of the coefficients and the spatial discretization method is the finite element method (FEM). The FEM will produce better results if the elements do not cross the discontinuity. NDSolve FEM functionality has to assume that the customer using NDSolve does not know anything about the FEM in fact she/he should not need to know that NDSolve uses the FEM. Having an automatic feature detection is a major advantage of using NDSolve and this type of functionality will certainly be extended upon.

That functionality like this exists is documented in the section Partial Differential Equations with Variable Coefficients.

One can switch that feature detection off by hiding the internal of the Piecewise behind a ?NumericQ. I'd advise against this. This will return a result that is less accurate then with the discontinuity detection on.

Now, from the details section of the documentation of DirichletCondition

"Dirichlet conditions are enforced at each point in the discretization of [PartialD][CapitalOmega] where pred is True."

What does this mean? During the conversion of the region what ever is on the boundary is considered for deploying the Dirichlet condition if True is requested. And here is the crux: On one hand the FEM gives better results with the discontinuity represented in the mesh and that has as a consequence that the discontinuity line is also part of the boundary - it is two different materials you are considering they do have boundaries also in nature.

I most certainly understand that you did not expect this. But you also expect the best possible result; so we have a conflict of interests so to speak. I am not sure how this can be resolved. A few things that could be considered: Add an example to the possible issues section of the documentation of DirichletCondition. Add a warning message that True might not be be what you'd expect. Let me know what you think.

So, I understand it's not expected but I would not go as far as saying it's a bug. You already put forward a good way to deal with this by giving an explicit predicate:

sol = NDSolveValue[{
    Laplacian[u[x, y], {x, y}] == f[x, y],
    DirichletCondition[u[x, y] == 0, {x, y} ∈ Circle[]]
    }, u, {x, y} ∈ Disk[]];

Plot3D[sol[x, y], {x, y} ∈ Disk[], PlotPoints -> 50]

Here is another example that (in Version 10.3) does not work out of the box but should and illustrates that automatic feature detection is indeed wanted:

NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 0, 
  DirichletCondition[u[x, y] == 1, x == 0 && y == 0]}, u, {x, 
   y} \[Element] Disk[]]

NDSolveValue::bcnop: "No places were found on the boundary where x==0&&y==0 was True, so DirichletCondition[u==1,x==0&&y==0] will effectively be ignored"

Because the automatic feature detection does not (yet) work for boundary conditions a message is given that the predicate of the Dirichlet condition inside the region can not be satisfied. Now, helping NDSolve a bit with:

if = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1, 
   DirichletCondition[u[x, y] == 1, x == 0 && y == 0]}, 
  u, {x, y} \[Element] Disk[], 
  Method -> {"FiniteElement", 
    "MeshOptions" -> {"IncludePoints" -> {{0, 0}}}}]

Will generate an interpolating function with the expected property:

if[0, 0]
1.0

If we look at the PointElements in the mesh we see that now the center coordinate is part of the mesh:

Needs["NDSolve`FEM`"]
ToBoundaryMesh[if["ElementMesh"]][
 "Wireframe"["MeshElement" -> "PointElements"]]

And a plot looks reasonable:

Plot3D[if[x, y], {x, y} \[Element] if["ElementMesh"], 
 PlotRange -> All]

This should have worked without the need to specify the additional point manually but NDSolve should have included it automatically. The point I like to make is that one does want NDSolve to make changes to the discretized region.

Answer 2

As I noted in a comment, redefining f in a way that does not change its value,

f[x_, y_] := Piecewise[{{-1, 0 <= x^2 + y^2 < 0.25}}, 0]

yields the desired solution, shown in the second plot in the question. The corresponding Wiremesh does not exhibit an internal boundary.

Addendum

In fact, all but one of the definitions I tried give the right answer:

f[x_, y_] := Piecewise[{{-1, (x^2 + y^2 <= 0.25) && (x | y) ∈ Reals}}, 0]

f[x_, y_] := Piecewise[{{-1, Abs[x^2 + y^2] < 0.25}}, 0]

f[x_, y_] := Piecewise[{{-1, x^2 + y^2 < 0.25 && -1 <= x && -1 <= y}}, 0]

f[x_, y_] := Piecewise[{{-1, -10^10 <= x^2 + y^2 <= 0.25}}, 0]

One does not, however.

f[x_, y_] := Piecewise[{{-1, -Infinity < x^2 + y^2 <= 0.25}}, 0]
sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == f[x, y],
    DirichletCondition[u[x, y] == 0, True]}, u, {x, y} ∈ Disk[]];
Plot3D[sol[x, y], {x, y} ∈ Disk[], PlotPoints -> 50, PlotRange -> All]

Although the shape seems correct, the peak value, sol[0, 0] = 0.106813, is too small by nearly 1/3.

Very strange!

Answer 3

Version 10.3.0:

If you just deny Mathematica's access to the interior of f, you get the expected result:

f[x_?NumericQ, y_] := Piecewise[{{-1, x^2 + y^2 < 0.25}}, 0]

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == f[x, y], 
DirichletCondition[u[x, y] == 0, True]}, u, {x, y} \[Element] Disk[]];

Plot3D[sol[x, y], {x, y} \[Element] Disk[], PlotPoints -> 50]

So, for me it's a bug. And even more since, when blocking access to bbgodfrey definition with the -Infinity, I also get the right answer:

f[x_?NumericQ, y_] := Piecewise[{{-1, -Infinity < x^2 + y^2 <= 0.25}}, 0]

---

I recently also had a case where blocking access to the interior of the function also solved the problem:

function3[t_] := If[t < 5, 5, 
  Interpolation[{{5, 5}, {10, 10}}, InterpolationOrder -> 1][t]]

NDSolveValue[{total[0] == 0, total'[t] == function3[t]}, total, {t, 0, 10}][10]
(*62.5*)

InterpolatingFunction::dmval: Input value {0} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

Right answer, but, in my opinion, the warning message shouldn't be there

By denying the access, with _?NumericQ, the warning goes away.

Wolfram Research support answered:

It is an expected behavior for numerical methods of finding solutions to differential equations to iterate in a vicinity of each point, e.g. 0 /-[Epsilon] or 5 /-[Epsilon] where [Epsilon] is a small positive real number. In order to avoid receiving the warning messages, one has to extend the definition of the functions to accommodate for /-[Epsilon] or limit the definition of the function to only evaluate for numerical values.

Which I get it, at the same time I don't...

Answer 4

The behavior of NDSolve described in the question concerns only the versions of Mathematica ealier than 11.1

Compare

Mathematica 11.1 :

Mathematica 11.2 :

This information was discovered in the Wolfram Technology Cofrerence 2017, Paritosh Mokhasi's presentation on Finite Elements , time=7:02 (or 5:25, I don't get reproducbile timing !). The slider title is : "Dirichlet conditions : Multiple materials"