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?