Making good meshes

Question

I have been looking at how to make meshes and how to control the size of the boundary elements and the interior elements. A relevent question may be found here. I can see how to do this for regions that are defined by implicit functions but not for regions that are defined by coordinates. I have made some minimum working examples to investigate before tackling my full problem.

For an implicit region I can do

Needs["NDSolve`FEM`"];
r1 = ImplicitRegion[x^2 + y^2 <= 100, {x, y}];
r2 = ImplicitRegion[x^2 + (y - 5)^2 <= 4, {x, y}];
reg = RegionDifference[r1, r2];
bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> 0.5];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 1];
mesh["Wireframe"]

This is a poor mesh if a uniform grid is required. To finish this off I made this dynamic to see how things adjust.

DynamicModule[{bL = 0.5, nn, bmesh, mesh, r = 10.2, mr, mcm, cL, 
  x0 = -10, y0 = 7},
 bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> bL];
 mesh = ToElementMesh[bmesh];
 Column[{
   Row[{"Boundary Edge Length = ", Slider[Dynamic[bL, {(bL = # ) &,
        (bL = # ; 
          bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> bL, 
            AccuracyGoal -> 1];

          mesh = ToElementMesh[bmesh, 
            MaxCellMeasure -> mcm]) &}], {0.05, 5}, 
      Appearance -> "Labeled"],
     "   Cell Area = ", Slider[Dynamic[mcm, {(mcm = # ) &,
        (mcm = # ; 
          bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> bL, 
            AccuracyGoal -> 1];

          mesh = ToElementMesh[bmesh, 
            MaxCellMeasure -> mcm]) &}], {0.05, 5}, 
      Appearance -> "Labeled"]}],
   Row[{"Approximate Cell Edge Length = ", 
     Dynamic[cL = Sqrt[mcm 4/Sqrt[3]]]}],
   Dynamic@
    Show[Graphics[{{Orange, 
        Triangle[{{x0, y0}, {x0 + cL, y0}, {x0 + cL/2, 
           y0 + cL Sqrt[3]/2}}], Line[{{x0, y0}, {x0, y0 + bL}}]}, 
       Circle[{0, 0}, r], 
       Table[Point[r {Cos[n 2 \[Pi]/nn], Sin[n 2 \[Pi]/nn]}], {n, 
         nn = Round[2 \[Pi] r/bL]}]}, ImageSize -> 10 72, 
      Frame -> True],
     mr = 
      MeshRegion[mesh, 
       MeshCellStyle -> {0 -> {PointSize[0.002], Red}, 1 -> Black, 
         2 -> White}]
     ]
   }]
 ]

The orange triangle and line are drawn to the size of the maximum cell measure and the length of the boundary element. I have also drawn a circle that is divided with the boundary edge length. It seems that the boundary element length cannot be larger than (approximately) the cell edge length. Two minor questions. 1) is the boundary edge length adjusted by the maximum cell measure to be an appropriate size? 2) I note there is a node at mid cell edges. This could be used for second order shape functions but I thought only linear interpolation was used. Is this correct?

My main question is about boundary meshes when coordinates are supplied. Here are some coordinates and a module for making a boundary mesh. I have again made a dynamic to look at changing the boundary cell measure and the element cell measure. This time the boundary cell measure does nothing,

Lx = 20; Ly = 30;
coords = {{0, 0}, {Lx, 0}, {Lx, Ly}, {0, Ly}, {Lx/4, Ly/2}, {Lx/2, Ly/
    2}, {Lx/2, (3 Ly)/4}, {Lx/4, 3 Ly/4}};
ClearAll[makeBM];
makeBM[bL_] := ToBoundaryMesh[
   "Coordinates" -> coords,
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        1}}], LineElement[{{5, 6}, {6, 7}, {7, 8}, {8, 5}}]},
   "RegionHoles" -> {{(3 Lx)/8, 5 Ly/8}},
   "MaxBoundaryCellMeasure" -> bL
   ];

DynamicModule[{bL = 0.5, nn, bmesh, mesh, r = 10.2, mr, mcm, cL, 
  x0 = -5, y0 = 25},
 bmesh = makeBM[bL];
 mesh = ToElementMesh[bmesh, MaxCellMeasure -> mcm];
 Column[{
   Row[{"Boundary Edge Length = ", Slider[Dynamic[bL, {(bL = # ) &,
        (bL = # ; bmesh = makeBM[bL];

          mesh = ToElementMesh[bmesh, 
            MaxCellMeasure -> {"Area" -> mcm}]) &}],
      {0.05, 5}, Appearance -> "Labeled"],
     "   Cell Area = ", Slider[Dynamic[mcm, {(mcm = # ) &,
        (mcm = # ; bmesh = makeBM[bL];

          mesh = ToElementMesh[bmesh, 
            MaxCellMeasure -> {"Area" -> mcm}]) &}], {0.05, 5}, 
      Appearance -> "Labeled"]}],
   Row[{"Approximate Cell Edge Length = ", 
     Dynamic[cL = Sqrt[mcm 4/Sqrt[3]]]}],
   Dynamic@
    Show[Graphics[{{Orange, 
        Triangle[{{x0, y0}, {x0 + cL, y0}, {x0 + cL/2, 
           y0 + cL Sqrt[3]/2}}], Line[{{x0, y0}, {x0, y0 + bL}}]}}, 
      ImageSize -> 10 72, Frame -> True],
     mr = 
      MeshRegion[mesh, 
       MeshCellStyle -> {0 -> {PointSize[0.002], Red}, 1 -> Black, 
         2 -> White}]
     ]
   }]
 ]

The boundary length is ignored. The main questions: How do you set the boundary mesh cell size here? Does this carry over to 3D meshes?

Edit Many thanks for your answer which does help to clarify. The suggestion of just using ToElementMesh does sound the way forward. However, I am still unclear about the role of MaxBoundaryCellMeasure when making a boundary mesh. The value seems to be changed when ToElementMesh sets a MaxCellMeasure. Here is an example. The boundary cell measure is set to be small with a length of 0.05. The MaxCellMeasure is set large with a value of 5 area units. The results is a large mesh on the boundary not a small mesh (for an equilateral triangle with area 5 the edge length is about 3.4 which is about the edge length I see. So in this example, unlike the first example boundary edge length has been changed.

Needs["NDSolve`FEM`"];
Lx = 20; Ly = 30;
coords = {{0, 0}, {Lx, 0}, {Lx, Ly}, {0, Ly}, {Lx/4, Ly/2}, {Lx/2, Ly/
    2}, {Lx/2, (3 Ly)/4}, {Lx/4, 3 Ly/4}};
bmesh = ToBoundaryMesh[
   "Coordinates" -> coords,
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        1}}], LineElement[{{5, 6}, {6, 7}, {7, 8}, {8, 5}}]},
   "RegionHoles" -> {{(3 Lx)/8, 5 Ly/8}},
   "MaxBoundaryCellMeasure" -> 0.05
   ];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5]["Wireframe"]

Hence from my viewpoint MaxBoundaryCellMeasure is not overwritten for implicit regions but is ignored for regions defined by coordinates. Is this correct?

Answer 1

Most importantly I'd like to recommend to take a look at the Element Mesh Generation tutorial. That tutorial deal with mesh generation for numerical applications like the Finite Element Method and, I think, it has answers to a few of your questions.

That said, let's look at the first example you provide and your claim that this is a poor mesh.

Needs["NDSolve`FEM`"];
r1 = ImplicitRegion[x^2 + y^2 <= 100, {x, y}];
r2 = ImplicitRegion[x^2 + (y - 5)^2 <= 4, {x, y}];
reg = RegionDifference[r1, r2];
bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> 0.5];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 1];
mesh["Wireframe"]

Now, an element mesh has documented methods and we look some of them:

Max[mesh["MeshElementMeasure"]]
0.999129

Looks as requested. Now, there is no method for the boundary element measure but for this case we can use:

Max[Sqrt[Total[(Subtract @@@ ((NDSolve`FEM`GetElementCoordinates[
            mesh["Coordinates"], #] & /@ 
          ElementIncidents[mesh["BoundaryElements"]])[[1, 
         All, {1, 2}]]))^2, {2}]]]
0.449962

Now, we look at the overall quality of the mesh:

Histogram[mesh["Quality"]]

Which, I think, does not look to bad either. To me this behaves as expected. If you are willing to have more elements you could use a higher "MeshQualityGoal":

mesh = ToElementMesh[bmesh, MaxCellMeasure -> 1, 
   "MeshQualityGoal" -> 1];

But I doubt that this will give you a better solution in a numerical usage. Here is why. The workflow you use is not optimal. Let me show you. Here is the exact integral over the region:

exact = Integrate[1, {x, y} \[Element] reg]
96 \[Pi]

Using your mesh:

exact - NIntegrate[1, {x, y} \[Element] mesh]
0.000277410831756697`

Using the symbolic region directly:

mesh2 = ToElementMesh[reg, "MaxBoundaryCellMeasure" -> 0.5, 
   MaxCellMeasure -> 1];
exact - NIntegrate[1, {x, y} \[Element] mesh2]
-0.00006363248667184962`

Using the symbolic region and an AccuracyGoal

mesh2 = ToElementMesh[reg, AccuracyGoal -> 2];
exact - NIntegrate[1, {x, y} \[Element] mesh2]
4.029767410429486`*^-6

The reason for this is explained in the Element Mesh Generation tutorial](https://reference.wolfram.com/language/FEMDocumentation/tutorial/ElementMeshCreation.html) and has to do with how the second order nodes are moved on the curved mesh. So, yes, an ElementMesh can have curved boundary but it needs to have the information of the curve and if you discretize the region with ToBoundaryMesh then ToElementMesh has no way of getting the information of the curvature of the region.

In version 11.1 you can also do something like

nr = ToNumericalRegion[reg];
manualBoundaryMesh = ToBoundaryMesh[reg, AccuracyGoal -> 2];
SetNumericalRegionElementMesh[nr, manualBoundaryMesh];
mesh3 = ToElementMesh[nr];
exact - NIntegrate[1, {x, y} \[Element] mesh3]
4.029767410429486`*^-6

This allows you to connect a manually generated boundary mesh to a region and thus get a good second order approximation of the region.

I am not sure I understand what your are trying to show with the dynamic. But to answer your questions.

1) Specifying a "MaxCellMeasure" will implicitly define an edge length for each mesh elements that will be overwritten if "MaxBoundaryCellMeasure" is given and we are dealing with a symbolic region.

2) The mid side node is used for second order elements and these can be curved. In your example they are not because ToElementMesh which generates the second order mesh has no information about the curvature. Note, that the visualization is linear though.

Your next example is a bit different as that has no curved edges so a linear boundary representation is exact.

3) When you specify a boundary mesh manually you do specify each edge element and currently (V11.1) "MaxBoundaryCellMeasure" does not split/reduce that further - it is ignored. (Maybe something for a future version). So for now, you need to specify the boundary granularity manually - Alternatively you could try to use a MeshRefinementFunction that changes the granularity closes to the boundary in question.

4) In 3D "MaxCellMeasure" corresponds to the volume of the element and the "MaxBoundaryCellMeasure" corresponds to the area of the boundary element.

Have a look at the Element Mesh Generation tutorial. If anything is unclear there let me know and I'll try to clarify things in the tutorial.