Problems using DensityPlot when solving Laplace's equation

I am trying to solve Laplace's equation in 2 dimensions with potencial boundary conditions at the edges of an external square and an internal circle. Everything seems to work fine until the point where I try to obtain a DensityPlot to visualize the the obtained potential at all points. When I execute the density plot line, I am returned with a diagram that displays the solution of the equation for only a fraction of the region where I solved the equation. How can I fix this ? The code is as follows:

α =
RegionDifference[Rectangle[{-50, -50}, {50, 50}], Disk[{0, 0}, 5]]
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100, x^2 + y^2 == 25],
u[x, -50] == u[x, 50] == u[-50, y] == u[50, y] == 0},
u, {x, y} ∈ α]

DensityPlot[sol[x, y], {x, y} ∈ α, Mesh -> None,
ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic] • I filed this as a bug. Nov 9 '18 at 8:26

It's a bad day for Regions, I guess...

Until NDSolve will dicretize it correctly, you can discretize the region yourself:

α = DiscretizeRegion[
RegionDifference[Rectangle[{-50, -50}, {50, 50}],
Disk[{0, 0}, 5]],
MaxCellMeasure -> {1 -> 1}
];
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100, x^2 + y^2 == 25],
u[x, -50] == u[x, 50] == u[-50, y] == u[50, y] == 0},
u, {x, y} ∈ α];
DensityPlot[sol[x, y], {x, y} ∈ α, Mesh -> None,
ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic] • Using DiscretizeRegion is not ideal for FEM based computations. See my answer for details. Nov 9 '18 at 8:23
• Thanks a lot for your continuous edits of my lazily written posts! Nov 9 '18 at 8:29
• Dito. That's my way of showing that I read them! ;) Nov 9 '18 at 8:30

Using DiscretizeRegion is not the best solution. This is because DiscretizeRegion will only provide a linear approximation for a region. A better approximation is to use ToElementMesh which can generate a second order accurate approximation to the region:

α = RegionDifference[Rectangle[{-50, -50}, {50, 50}], Disk[{0, 0}, 5]]
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100, x^2 + y^2 == 25],
u[x, -50] == u[x, 50] == u[-50, y] == u[50, y] == 0},
u, {x, y} ∈ α]
sol["ElementMesh"]["Wireframe"] DensityPlot understands ElementMesh:

DensityPlot[sol[x, y], {x, y} ∈ sol["ElementMesh"],
Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic] Now to look at the approximation quality difference:

N[exact = Integrate[1, {x, y} ∈ α]]
9921.460183660256

exact - Area[DiscretizeRegion[α]]
-2.610161994172813

exact - Total[sol["ElementMesh"]["MeshElementMeasure"], 2]
-0.004173680044914363`

Clearly the use of an ElementMesh is superior here as it much better approximates the curved boundaries.

This is explained in more detail in the section Comparing ElementMesh and MeshRegion in the Element Mesh Generation tutorial.