Distinguish between inner and outer boundary in MeshRegion

Question

For example, this is a 2D MeshRegion with one (or more) holes in it:

mr = DiscretizeRegion@RegionDifference[
   Rectangle[{-1, -1}, {1, 1}],
   Disk[{0, 0}, 0.5]
]

I would like to identify outer boundary (MeshCells) and inner boundary of a hole. I can get indices of all boundary MeshCells, but cannot distinguish between them. How can I do that?

HighlightMesh[mr, MeshCellIndex[mr, {1, "Boundary"}]] 

Ultimately I will convert this MeshRegion to ElementMesh and assign different integer markers to LineElements on the boundary.

Answer 1

You can use ConnectedMeshComponents as follows:

cmc = ConnectedMeshComponents[DiscretizeGraphics[MeshPrimitives[mr, {1, "Boundary"}]]]; 

Show[mr,  Epilog -> Thread[{Thickness[.01], {Red, Blue},  MeshPrimitives[#, 1] & /@ cmc}]]

Answer 2

You can convert your MeshRegion to a boundary element mesh with auto generated markers by:

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh[RegionBoundary[mr]];
bmesh["Wireframe"["MeshElementStyle" -> "BoundaryGrouping"]]

The fact that the rectangle has also been split is desirable for FEM. You can look at the markers generated.

bmesh["BoundaryElementMarkerUnion"]
{1, 2, 3, 4, 5}

Here is a mesh that has an outer boundary, an inner boundary and a material boundary:

Needs["NDSolve`FEM`"]
\[CapitalOmega] = 
  ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= (1/
        2)^2 && (x + 1/2)^2 + (y + 1/2)^2 >= (1/2)^2, {{x, -2, 
     2}, {y, -2, 2}}];
(mesh = ToElementMesh[\[CapitalOmega], 
    "RegionHoles" -> {{-1/2, -1/2}}, 
    "RegionMarker" -> {{{1/2, 1/2}, 2}, {{1/2, -1/2}, 1}}, 
    "MaxBoundaryCellMeasure" -> 0.125])["Wireframe"]

You can then extract the boundary element marker union; these markers are auto generated.

groups = mesh["BoundaryElementMarkerUnion"]
{1, 2, 3, 4, 5, 6}

Some colors for the markers:

temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp

mesh["Wireframe"["MeshElement" -> "BoundaryElements", 
  "MeshElementStyle" -> (Directive[#] & /@ colors)]]

To further work with this you could use MeshElementSplitByMarker:

mesbm = Flatten[
   MeshElementSplitByMarker[
    MeshOrderAlteration[mesh, 1]["BoundaryElements"]]];

and extract one of the components:

Graphics[GraphicsComplex[mesh["Coordinates"], 
  Line[ElementIncidents[mesbm[[6]]]]]]

This selection works by looking at normals and grouping them with some threshold. There are a few examples in the ElementMesh ref page in the Scope section under "BoundaryNormals".

Another approach is to use the "VertexBoundaryConnectivity" in conjunction with StronglyConnectedComponents.

bmesh = ToBoundaryMesh[mesh];
belems = bmesh["BoundaryElements"];

bincidents = Join @@ ElementIncidents[belems];
S = bmesh["VertexBoundaryConnectivity"];
S = Transpose[S].S;
components = SparseArray`StronglyConnectedComponents[S];
pls = Map[Line[Developer`ToPackedArray /@ (Part[bincidents, #])] &, 
  components]
Graphics[GraphicsComplex[bmesh["Coordinates"], pls]]