# 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. – user21 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. – user21 Nov 9 '18 at 8:23
• Thanks a lot for your continuous edits of my lazily written posts! – user21 Nov 9 '18 at 8:29
• Dito. That's my way of showing that I read them! ;) – Henrik Schumacher 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.