11
$\begingroup$

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"]

enter image description here

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]

enter image description here

$\endgroup$
0

3 Answers 3

9
$\begingroup$

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]]

enter image description here

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]]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks for your reply, but I'm a little confused. I thought DensityPlot just postprocesses the solution. By doing what you showed, is the same mesh used for solving the problem? $\endgroup$
    – Jiangming
    Mar 18, 2020 at 1:15
  • $\begingroup$ @Jiangming, I am not familiar with how DensityPlot works. For your example case, if you check sol["Grid"] == bmesh["Coordinates"] and sol["ElementMesh"] == bmesh both yield True. So it seems DensityPlot with sol["Grid"] and refines it further. How PlotPoints enters the picture I don't have a clue. $\endgroup$
    – kglr
    Mar 18, 2020 at 2:01
9
$\begingroup$

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"]]

enter image description here

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?

$\endgroup$
8
$\begingroup$

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]

FEM Mesh and Solution

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.