Why the mesh generated for NDSolve is different from the final plotted mesh?

Question

When I try the following code to solve the Laplace's equation, the mesh generated by bmesh = ToElementMesh[bmesh] is different from the mesh shown in DensityPlot. Why?

<< NDSolve`FEM`
region = RegionDifference[Rectangle[{0, 0}, {100, 100}],Disk[{50, 50}, 10]];
bmesh = ToBoundaryMesh[region];
bmesh = ToElementMesh[bmesh];
bmesh["Wireframe"]

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
   DirichletCondition[u[x, y] == 100, Sqrt[(x - 50)^2 + (y - 50)^2] == 10],
   u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0}, u, {x, y} \[Element] bmesh]
DensityPlot[sol[x, y], {x, y} \[Element] bmesh, Mesh -> All, ColorFunction -> "Rainbow", PlotRange -> All, PlotLegends -> Automatic]

Answer 1

If you add the options MaxRecursion -> 0 (to prevent recursive subdivisions) and PlotPoints -> 10 (10 was a lucky guess:) to DensityPlot the mesh lines match the wireframe of bmesh.

Compare DensityPlot output with ColorFunction -> (White&) (middle plot) with bmesh["Wireframe"] (right plot):

Row[{dp1 = DensityPlot[sol[x, y], {x, y} ∈ bmesh, 
    PlotPoints -> 10, MaxRecursion -> 0, 
    Mesh -> All, MeshStyle -> Black, PlotRangeClipping -> False, 
    PlotRangePadding -> 0, Frame -> False, ColorFunction -> "Rainbow",
    PlotRange -> All, ImageSize -> 1 -> 3], 
  dp2 = DensityPlot[sol[x, y], {x, y} ∈ bmesh, 
    PlotPoints -> {10, 10}, MaxRecursion -> 0,
    Mesh -> All, PlotRangeClipping -> False, MeshStyle -> Black, 
    PlotRangePadding -> 0, Frame -> False, ColorFunction -> (White &),
    PlotRange -> All, ImageSize -> 1 -> 3], 
  Show[bmesh["Wireframe"], ImageSize -> 1 -> 3]}, Spacer[3]]

Overlay dp2 with bmesh["Wireframe"] to see the matching:

Legended[Show[bmesh["Wireframe"] /. p_Polygon :> 
 {EdgeForm[{JoinForm["Round"], AbsoluteThickness[5], Red}], FaceForm[], p}, 
  dp2 /. Polygon[a_, ___] :> {Thick, Black, Line /@ (Append[#, First@#] & /@ a)},  
    ImageSize -> 1 -> 5], 
 Placed[LineLegend[{Red, Black}, 
    {"bmesh[\"Wireframe\"]", "DensityPlot mesh"}], Right]]

Answer 2

A few additions to kglr's and Tim's answers.

First, I'd like to point out that the way you generate the mesh is not optimal. Consider your setup:

<< NDSolve`FEM`
region = RegionDifference[Rectangle[{0, 0}, {100, 100}], 
   Disk[{50, 50}, 10]];
bmesh = ToBoundaryMesh[region];
mesh = ToElementMesh[bmesh];
(*mesh["Wireframe"]*)

Now, we compute the symbolic area of the region:

area = Area[region] // N
9685.84

and we compare that to the area of the mesh:

area - Total[mesh["MeshElementMeasure"], 2]
-5.14482

So we are off by a bit. We can do better by using the symbolic region as an input to ToElementMesh directly:

region = RegionDifference[Rectangle[{0, 0}, {100, 100}], 
   Disk[{50, 50}, 10]];
mesh = ToElementMesh[region];
Area[region] - Total[mesh["MeshElementMeasure"], 2]
-0.00636722

Which is substantially better. The reason is the following. When you give a symbolic region to ToElementMesh the mid side nodes of the second order mesh can be moved to the proper position given by the symbolic representation of the region. When you use a boundary element mesh then that information is lost, since there an edge is a straight line. For more information on this and how to still use a boundary mesh to get a high quality approximation to a region please see the ElementMesh generation tutorial.

If you still need a bounday mesh you could get it by using:

bmesh = ToBoundaryMesh[mesh];

Now, concerning your question about using the FEM mesh for visualization. This is not necessarily a good idea. The requirements for a good visualization mesh are different to that for a good FEM mesh. For example, a good visualization mesh needs to resolve areas with steep gradients. While this might also be a desirable feature for FEM this is not available right now. But you can tweak the plotting functions to accept the FEM mesh as the visualization mesh.

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
    DirichletCondition[u[x, y] == 100, 
     Sqrt[(x - 50)^2 + (y - 50)^2] == 10], 
    u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0}, 
   u, {x, y} \[Element] mesh];
dp = DensityPlot[sol[x, y], {x, y} \[Element] mesh, Mesh -> All, 
   ColorFunction -> "Rainbow", PlotRange -> All, 
   PlotLegends -> Automatic, MaxRecursion -> 0, PlotPoints -> 14];
Show[dp, mesh["Wireframe"]]

The MaxRecursion->0 prevents adaptive mesh refinement for the visualization functions. The PlotPoints -> 14 limit the number of sample points used in each direction. This will also trigger a refinement in the visualization functions. For my up to 14 plot points did not alter the mesh. So this setting requites some experimenting.

Another alternative is to use:

Show[
 ElementMeshContourPlot[sol, ColorFunction -> "Rainbow", Mesh -> All]
 , mesh["Wireframe"], PlotPoints -> 45, PlotLegends -> Automatic]

ElementMeshContourPlot extracts the FEM mesh from the interpolating function and uses that. However, this also has drawbacks, as options like PlotPoints and PlotLegends do not work and this functionality does not exist for DensityPlot. Would you find something like that useful?

Answer 3

I think your main issue is that the FEM mesh is second order and the plotting mesh is first order. You could change the element order of your FEM mesh to first order, but that will affect accuracy.

You can see that the FEM mesh has 6 coordinates per triangle by the following:

<< NDSolve`FEM`
region = RegionDifference[Rectangle[{0, 0}, {100, 100}], 
   Disk[{50, 50}, 10]];
bmesh = ToBoundaryMesh[region];
mesh = ToElementMesh[bmesh];
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
    DirichletCondition[u[x, y] == 100, 
     Sqrt[(x - 50)^2 + (y - 50)^2] == 10], 
    u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0}, 
   u, {x, y} \[Element] mesh];
dp = DensityPlot[sol[x, y], {x, y} \[Element] mesh, PlotPoints -> All,
    MaxRecursion -> 4, ColorFunction -> "Rainbow", PlotRange -> All, 
   PlotLegends -> Automatic];
lp = ListPlot[mesh["Coordinates"], PlotMarkers -> Automatic, 
   PlotStyle -> Red];
Show[dp, mesh["Wireframe"], lp, ImageSize -> Large]