NDSolveValue Quality Issue

Question

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.

Answer 1

Alternatively, you can use

Needs["NDSolve`FEM`"]
\[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"]

Answer 2

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]