# Generating Meshes with Conforming Boundaries for FEM

I would like to solve a PDE on an annulus with some boundary condition that relates the inner circle to the outer circle. (For example, periodicity along radial lines). If I use a finite element method, it would be convenient if the inner circle and outer circle have conforming boundary meshes. That is, the nodes on the inner and outer circles match. Is there a way to generate an ElementMesh object with this property?

I would prefer triangular elements so that I can apply this to a wider range of problems. More generally, if I have a region with 2n boundaries related by n maps, is there a way to generate a conforming mesh?

• Do you mean the same number of nodes for the inner and outer circles? – Alexei Boulbitch Jul 28 '15 at 11:40
• If F is a map from the inner circle to the outer circle, and x_i are the nodes on the inner circle, I need F[x_i] to be the nodes on the outer circle. This is more than just saying that the number of nodes are the same. – qdice Jul 28 '15 at 12:44
• For your general problem -- be aware that it is easy to introduce an orientation non-preserving map between two boundary components. This will embed a Moebius strip in your region (if, for instance, the boundary condition is equality on the two components). As a consequence, odd modes along paths in the region starting at one associated point and ending in the others will be driven to zero by topology. Are you sure you can reliably make properly oriented maps when this will matter? – Eric Towers Jul 28 '15 at 16:50
• Will the maps between your boundary components have (approximately) unit magnitude of derivative? If not, then you may not be able to get a sufficiently uniform grid when a tiny piece of one boundary is mapped to a long piece of another. – Eric Towers Jul 28 '15 at 17:16
• Can your map between boundary components be "thickened" to include nonintersecting neighbourhoods of each boundary component? (If so, an abstract solution is to transport one neighbourhood along the map to the other component, grid the target and moved neighbourhoods, map back, and repeat for all the pairs. The result leaves only the whole region minus the neighbourhoods to grid compatibly with the "outer" boundaries of the neighbourhoods -- i.e., not compatibly requirements need be enforced. I'm sure I don't want to implement this method.) – Eric Towers Jul 28 '15 at 17:22

Let b1, b2 map the unit circle to the inner, outer boundaries respectively. You might need to handle a list of points with Map (e.g., b1 /@ bdy) depending on b1 and b2. Then here is a simple way:

Needs["NDSolveFEM"]

Clear[b1, b2];
n = 100;
bdy = Table[{Cos[t], Sin[t]}, {t, 0., 2 Pi - 2 Pi/(2 n),
2 Pi/(2 n)}];
hole = {0., 0.};
b1[pts_?MatrixQ | pts_?VectorQ] := 2 pts;  (* to handle list of points or single pts *)
b2[pts_?MatrixQ | pts_?VectorQ] := 3 pts;
bmesh = ToBoundaryMesh[
"Coordinates" -> Join[b1@bdy, b2@bdy],
"BoundaryElements" -> {
LineElement[Partition[Range[2 n], 3, 2, 1][[All, {1, 3, 2}]]],
LineElement[
Partition[Range[2 n + 1, 4 n], 3, 2, 1][[All, {1, 3, 2}]]]
},
"RegionHoles" -> b1@hole
];
emesh = ToElementMesh[bmesh];

Show[
emesh["Wireframe"],
Graphics[{Red, Point[Join[b1@bdy, b2@bdy]]}]
]


I picked n = 100 to show the limitations of this somewhat naive approach. The interior mesh does not match the fineness of the boundary mesh. You can control it with MaxCellMeasure. However, depending on the setting, it might subdivide the boundary, in which case there is no guarantee (as far as I know) that the boundaries will conform. Here is a setting that seems to leave the boundary elements intact:

emesh = ToElementMesh[bmesh,
MaxCellMeasure -> {"Length" -> Sqrt[2] 3 2. Pi/n}];
Show[
emesh["Wireframe"],
Graphics[{Red, Point[Join[b1@bdy, b2@bdy]]}]
]


Maybe that will get you started.

Alternative:

Suppose you have one of the boundaries, and mapping bmap from it to the other boundary, as well as a function hole that identifies the hole in the region. It is better if there is a symbolic description of the region; if there is, then the boundary mesh can be generated accurately. But it the region is from some sort of measurement (such as an image), a discrete representation is fine. Then we can mesh the outer boundary (or whichever is longer) and map the points to the inner boundary. If we mesh the outer boundary with a certain MaxCellMeasure, then the region with the hole is theoretically mesh-able with the same MaxCellMeasure. I can't prove whether it won't result in a modification of the boundary mesh (maybe someone else can comment). It can be checked afterwards though. One can increase the MaxCellMeasure slightly on the second run, to give the mesher some wiggle room.

(* this first part simply constructs the boundary *)
maxCellMu = 0.01;
mesh0 = ToElementMesh[Disk[{0, 0}, 3], MaxCellMeasure -> maxCellMu];
belem = Flatten@ ElementIncidents[mesh0["BoundaryElements"]][[All, All, {1, 3}]];
bpts = mesh0["Coordinates"] ~Part~ belem;
n = Length[bpts]/2;

Clear[bmap, hole];
bmap[pts_?MatrixQ | pts_?VectorQ] := 2/3 pts;
hole[bdy_] := Mean[bdy];    (* works, e.g., for a convex boundary of the hole *)
bmesh = ToBoundaryMesh[
"Coordinates" -> Join[bpts, bmap@bpts],
"BoundaryElements" -> {
LineElement[Partition[Range[2 n], 3, 2, 1][[All, {1, 3, 2}]]],
LineElement[Partition[Range[2 n + 1, 4 n], 3, 2, 1][[All, {1, 3, 2}]]]
},
"RegionHoles" -> hole[bmap@bpts]
];
emesh = ToElementMesh[bmesh, MaxCellMeasure -> (*1.05*) maxCellMu];
emesh["Wireframe"]


Issue:

I discovered that, while the vertices (end nodes) of the boundary LineElements are preserved, the mid-side nodes are moved to the midpoints (mean) of the endpoints. This results in a loss of accuracy. I could not discover a method for preserving them, so the solution I propose below is to move them back after the ElementMesh has been generated. Perhaps someone will find a better way.

adjustBoundary[m_ElementMesh, b_] :=
incidents = ElementIncidents[m["BoundaryElements"]];
quadpoints = m["Coordinates"] ~Part~ Flatten@incidents[[All, All, 3]];
m /. mapping
]



We can visualize the difference. (As explained in FEMDocumentation/tutorial/ElementMeshVisualization, "Wireframe" displays the linear mesh elements as straight lines.)

showmeshpoints[m_, bdy_, opts__?OptionQ] := Show[
Graphics[{Red, PointSize[Medium], Point@bdy}],
emesh["Wireframe"],
Graphics[{Blue, PointSize[Medium], Point[
m["Coordinates"]~Part~
Flatten@ElementIncidents[m["BoundaryElements"]]]}],
opts
]

GraphicsRow[
showmeshpoints[#, Join[bpts, bmap@bpts],
PlotRange -> {{-0.02, 0.16}, {2.94, 3.01}},
AspectRatio -> 1] & /@ {emesh, amesh}
]


The blue points are the boundary points from the respective mesh (emesh left, amesh right). The red points, when not covered by a blue point, reveal a boundary point that has been moved in the mesh from the original boundary point.

• I ran this and checked emesh["Coordinates"] and emesh["BoundaryElements"]. It doesn't seem like these elements conform. – qdice Jul 28 '15 at 12:32
• On second inspection, I think this does in fact conform, but the ordering of the nodes on emesh["BoundaryElements"] is different between the two circles. I think I could still work with this. Thank you. – qdice Jul 28 '15 at 13:13
• @qdice First, yes, I'm not sure the order can be controlled, if you let automatic meshing occur. Second, I think I figured out something else that might have been confusing you about the correspondence of the boundary points. (See "Issue" in update.) Third, I added an alternative that conforms more closely to the specs in your question. – Michael E2 Jul 28 '15 at 15:20
Needs["NDSolveFEM"]

numRings = 10;
numSectors = 40;
points =
Flatten[
Table[
r {Cos[ϕ], Sin[ϕ]},
{r, 1, numRings},
{ϕ, 0, 2 π - 2 π/numSectors, 2 π/numSectors}
],
1
];

triangles =
Flatten[
Table[
{
{
numSectors (i + 1) + j,
numSectors i + Mod[1 + j, numSectors, 1],
numSectors i + j
},
{
numSectors (i + 1) + j,
numSectors (i + 1) + Mod[1 + j, numSectors, 1],
numSectors i + Mod[1 + j, numSectors, 1]
}
},
{i, 0, numRings - 2}, {j, 1, numSectors}
], 2];

mesh = ToElementMesh["Coordinates" -> points,
"MeshElements" -> {TriangleElement[triangles]}
];

mesh["Wireframe"]


points = Flatten[