Problems using DensityPlot when solving Laplace's equation

Question

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]

Answer 1

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]

Answer 2

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.