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.
f[x_, y_] := Piecewise[{{-1, 0 <= x^2 + y^2 < 0.25}}, 0]
solves the problem! $\endgroup$