# NDSolveValue Quality Issue

I'm having trouble with the resolution of the NDSolveValue solution.

For example this code from here gives the following on my machine (Windows, Version 10.4.1): (Zoomed in to show detail)

Ω = RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]];
sol = NDSolveValue[{
D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100.,
x == 40 && 40 <= y <= 60 ||
x == 60 && 40 <= y <= 60 ||
40 <= x <= 60 && y == 40 ||
40 <= x <= 60 && y == 60],
u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0
}, u, {x, y} ∈ Ω]

DensityPlot[sol[x, y], {x, y} ∈ Ω,  Mesh -> None, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends->Automatic]


This issue persists even when I use:

Method -> {"FiniteElement",
"MeshOptions" -> {"BoundaryMeshGenerator" -> "Continuation"}}


I've tried reducing the MaxCellMeasure to 0.01, increasing AccuracyGoal, PrecisionGoal, and WorkingPrecision, set InterpolationOrder to All. None of these changes seem to produce better results.

Alternatively, you can use

Needs["NDSolveFEM"]
\[CapitalOmega] =
ToElementMesh[
RegionDifference[Rectangle[{0, 0}, {100, 100}],
Rectangle[{40, 40}, {60, 60}]], MaxCellMeasure -> {"Area" -> 1}]


This will generate a second order mesh and you'll need far fewer elements to get an accurate solution then with a first order mesh (as DiscretizeRegion generates)

Ideally you would be able to give the MaxCellMeasure to NDSolve as a FEM option directly, but that does not seem to work. I'll need to investigate that a bit.

Update

This is actually not a problem of NDSolve but of the plotting function. If you discretize the region omega then the plotting function (DensityPlot) will also see the refined mesh. If the you refine the mesh in NDSolve and then extract that mesh and pass it to DensityPlot in place of the undiscretized region omega you'll get the result you expect.

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100.,
x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= 60 ||
40 <= x <= 60 && y == 40 || 40 <= x <= 60 && y == 60],
u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0},
u, {x, y} \[Element] \[CapitalOmega],
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 1}}]

DensityPlot[sol[x, y], {x, y} \[Element] sol["ElementMesh"],
Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic]


Note that here I used \[Element] sol["ElementMesh"]

• When I was struggling with the FEA options in NDSolve, it didn't seem like they were effective. I even tried the 2nd Order mesh inside NDSolve and it had some minor effect but without the ability to change MaxCellMeasure it was almost useless. It also seems like the mesh generation is significantly faster outside of NDSolve. – Young Jun 27 '16 at 3:22
• @Young, I'll investigate that. – user21 Jun 27 '16 at 3:23
• @Young, thank you for taking the time to reformulate the comment question and put that into a full question. – user21 Jun 27 '16 at 3:27

Using DiscretizeRegion[] with RegionDifference[] prior to NDSolveValue was the key.

Ω =
DiscretizeRegion[
RegionDifference[Rectangle[{0, 0}, {100, 100}],
Rectangle[{40, 40}, {60, 60}]],
Method -> "Continuation", AccuracyGoal -> 7, PrecisionGoal -> 7,
MaxCellMeasure -> {"Area" -> 0.1}];

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100.,
x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= 60 ||
40 <= x <= 60 && y == 40 || 40 <= x <= 60 && y == 60],
u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0},
u, {x, y} ∈ Ω, InterpolationOrder -> All]

DensityPlot[sol[x, y], {x, y} ∈ Ω,
Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic] 