# Prelude
Mathematica 11.3, Windows.
Let us say I have a rectangular region with 10 little holes inside. Their coordinates are given by `positionList`. The resulting region `mr` is made through `RegionDifference` and `RegionUnion` of disks and rectangles:
Needs["NDSolve`FEM`"]
radius = 0.1;
Louter = 2;
OuterRegion = Rectangle[{-L/2, -L/2}, {L/2, L/2}] /. L -> Louter;
positionList = {{0.706`, -0.14`}, {0.389`, 0.593`}, {-0.278`,
0.429`}, {0.254`,
0.844`}, {-0.46`, -0.367`}, {0.737`, -0.759`}, {-0.07`, -0.664`}, \
{-0.469`, 0.626`}, {-0.755`, -0.509`}, {-0.455`, -0.015`}};
diskRegions = Disk[#, radius] & /@ positionList;
mr = RegionDifference[OuterRegion, RegionUnion@diskRegions];
LP = ListPlot[positionList -> Table[k, {k, 1, Length@positionList}]];
RP = Show[mr // RegionPlot, LP, ImageSize -> 400]
[![enter image description here][1]][1]
The labels of the disks show the order in which they appear in `positionList`. Ultimately I want a FEM mesh where I can apply specific boundary conditions on disk 1, disk 2, disk 3 etc. First, I make a boundary mesh:
bmesh = ToBoundaryMesh[RegionBoundary@mr,
"RegionHoles" -> positionList, "BoundaryGroupingThreshold" -> 0,
AccuracyGoal -> 1];
Show[bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue]],
ImageSize -> 400]
[![enter image description here][2]][2]
# Problem
As you can see, the **Region Markers are scrambled**! In the boundary mesh, the hole corresponding to disk 2 (in the blue plot above) actually has marker 14. Hole 3 has marker 13, and so on.
**I want the markers of the boundary mesh to be consistent with the disk numbers, here going from 1 to 10. The markers on the walls should have numbers 11-14.**
# A manual solution
I have solved this problem with a `BoundaryMarkerFunction`, writing a condition manually for each disk and each of the walls:
boundaryMarkerFunction =
Compile[{{boundaryElementCoords, _Real,
3}, {pointMarkres, _Integer, 2}},
Module[{pt1 = #[[1]], pt2 = #[[2]]},
Which[
(pt1[[1]] - positionList[[1]][[1]])^2 + (pt1[[2]] -
positionList[[1]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[1]][[1]])^2 + (pt2[[2]] -
positionList[[1]][[2]])^2 < 1.1 radius^2 , 1,
(pt1[[1]] - positionList[[2]][[1]])^2 + (pt1[[2]] -
positionList[[2]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[2]][[1]])^2 + (pt2[[2]] -
positionList[[2]][[2]])^2 < 1.1 radius^2 , 2,
(pt1[[1]] - positionList[[3]][[1]])^2 + (pt1[[2]] -
positionList[[3]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[3]][[1]])^2 + (pt2[[2]] -
positionList[[3]][[2]])^2 < 1.1 radius^2 , 3,
(pt1[[1]] - positionList[[4]][[1]])^2 + (pt1[[2]] -
positionList[[4]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[4]][[1]])^2 + (pt2[[2]] -
positionList[[4]][[2]])^2 < 1.1 radius^2 , 4,
(pt1[[1]] - positionList[[5]][[1]])^2 + (pt1[[2]] -
positionList[[5]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[5]][[1]])^2 + (pt2[[2]] -
positionList[[5]][[2]])^2 < 1.1 radius^2 , 5,
(pt1[[1]] - positionList[[6]][[1]])^2 + (pt1[[2]] -
positionList[[6]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[6]][[1]])^2 + (pt2[[2]] -
positionList[[6]][[2]])^2 < 1.1 radius^2 , 6,
(pt1[[1]] - positionList[[7]][[1]])^2 + (pt1[[2]] -
positionList[[7]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[7]][[1]])^2 + (pt2[[2]] -
positionList[[7]][[2]])^2 < 1.1 radius^2 , 7,
(pt1[[1]] - positionList[[8]][[1]])^2 + (pt1[[2]] -
positionList[[8]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[8]][[1]])^2 + (pt2[[2]] -
positionList[[8]][[2]])^2 < 1.1 radius^2 , 8,
(pt1[[1]] - positionList[[9]][[1]])^2 + (pt1[[2]] -
positionList[[9]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[9]][[1]])^2 + (pt2[[2]] -
positionList[[9]][[2]])^2 < 1.1 radius^2 , 9,
(pt1[[1]] - positionList[[10]][[1]])^2 + (pt1[[2]] -
positionList[[10]][[2]])^2 <
1.1 radius^2 && (pt2[[1]] -
positionList[[10]][[1]])^2 + (pt2[[2]] -
positionList[[10]][[2]])^2 < 1.1 radius^2 , 10,
pt1[[1]] < -0.99 && pt2[[1]] < -0.99, 11, (* left wall *)
pt1[[2]] < -0.99 && pt2[[2]] < -0.99, 12, (* bottom wall *)
pt1[[1]] > 0.99 && pt2[[1]] > 0.99, 13, (* right wall *)
pt1[[2]] > 0.99 && pt2[[2]] > 0.99, 14, (* top wall *)
True, 4 ]] & /@ boundaryElementCoords];
bmesh = ToBoundaryMesh[RegionBoundary@mr,
"RegionHoles" -> positionList, "BoundaryGroupingThreshold" -> 0,
AccuracyGoal -> 2,
"BoundaryMarkerFunction" -> boundaryMarkerFunction];
{Show[bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue]],
ImageSize -> 400], RP}
[![enter image description here][3]][3]
[![enter image description here][4]][4]
# An automatic solution?
Could this be automated, so I can vary the number of disks? I would have happily replaced the `Which` construct by some `Table`, but I am not sure what kind of arguments this mysterious `boundaryMarkerFunction` takes. For instance, I tried
boundaryMarkerFunction[bmesh["BoundaryElements"],
bmesh["Coordinates"]]
which fails because `boundaryElementCoords` is supposed to be 3D.
**Please help me find an automatic solution to label the boundary elements in my desired order.**
# Just for fun
This allows us to use `ElementMarker` with `NDSolve`, specifying precise individual boundary conditions on the disks. For instance, here is the heat equation with a different `DirichletCondition` on each disk:
bmesh = ToBoundaryMesh[RegionBoundary@mr,
"RegionHoles" -> positionList, "BoundaryGroupingThreshold" -> 0,
AccuracyGoal -> 4,
"BoundaryMarkerFunction" -> boundaryMarkerFunction];
mesh = ToElementMesh[bmesh];
boundaryMarkerList = Table[k, {k, 1, Length@positionList}];
op = - Laplacian[u[x, y], {x, y}] + 0.1 u[x, y];
BCedges = {DirichletCondition[u[x, y] == 0, ElementMarker == 11],
DirichletCondition[u[x, y] == 1, ElementMarker == 13]};
BCcircles =
Table[DirichletCondition[u[x, y] == RandomReal[{0, 1}],
ElementMarker == k], {k, 1, Length@positionList}];
BC = Join[BCedges, BCcircles];
ufun = NDSolveValue[{op == 0, BC}, u, {x, y} \[Element] mesh]
{
DensityPlot[ufun[x, y], {x, y} \[Element] mesh,
ColorFunction -> "Temperature", AspectRatio -> Automatic,
PlotPoints -> 100, PlotRange -> All, ImageSize -> 400],
Show[mesh["Wireframe"], ImageSize -> 400]
}
[![enter image description here][5]][5]
[![enter image description here][6]][6]
[1]: https://i.sstatic.net/QSpyc.png
[2]: https://i.sstatic.net/rujEq.png
[3]: https://i.sstatic.net/YRhyo.png
[4]: https://i.sstatic.net/hgJDO.png
[5]: https://i.sstatic.net/b6t8g.png
[6]: https://i.sstatic.net/GXAYY.png