# Interior of discretized region

I have some domain of interest:

r = 0.3;
dom = ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= r^2, {{x, 0, 1}, {y, 0, 1}}];


Now I can extract mesh of the whole domain or just boundary:

Needs["NDSolveFEM"]
Boundary = First@ToBoundaryMesh[dom];
Domain = First@ToElementMesh[dom];


If I ListPlot those two, Boundary is what I expected: All points that lie on the boundary of dom. However, can I somehow extract only the interior points of dom? That is, something like "Domain - Boundary"? Is there a function that does it in a nice way? I'm really interested only in pure positions of points, that's why I'm taking First@ of the whole dataset.

Thanks.

It can be done like this:

Needs["NDSolveFEM"]
mesh = ToElementMesh[dom];

boundaryIncidents = Union@Flatten@ElementIncidents[mesh["BoundaryElements"]];
allIncidents = Union@Flatten@ElementIncidents[mesh["MeshElements"]];
interiorIncidents = Complement[allIncidents, boundaryIncidents];

allCoordinates = mesh["Coordinates"];
boundaryCoordinates = allCoordinates[[boundaryIncidents]];
interiorCoordinates = allCoordinates[[interiorIncidents]];

Show[
mesh["Wireframe"],
Graphics[{
PointSize[Large],
Red, Point[boundaryCoordinates],
Blue, Point[interiorCoordinates]
}]
]


You can think of incidents as the indices of coordinates in the list mesh["Coordinates"]. The idea is the same as that of GraphicsComplex if you have ever used that. mesh["BoundaryElements"] gives you all the line elements that make up the boundary described using incidents. mesh["MeshElements"] gives you all the triangles that make up the mesh described using the same incidents. All we have to do to get the coordinates that are not on the boundary is take the incidents of the full mesh and remove the incidents that are used in the boundary elements.

Here is another approach:

mesh = ToElementMesh[dom, "MeshOrder" -> 2];
bmesh = ToBoundaryMesh[mesh, "MeshOrder" -> 2];
boundary = bmesh["Coordinates"];
interior = Complement[mesh["Coordinates"], boundary];


I specify "MeshOrder" just to make sure it will be the same in the mesh and the boundary mesh. I noticed that by default the boundary mesh gets mesh order 1 and the full mesh gets mesh order 2.

Another option, provided by user21:

mesh = ToElementMesh[dom, "MeshOrder" -> 2];
coords = mesh["Coordinates"];
boundaryIDs = Flatten[ElementIncidents[mesh["PointElements"]]];
interiorIDs = Complement[Range[Length[coords]], boundaryIDs];

• Thank you, the second solution is actually very short and elegant! Aug 26, 2017 at 12:29
• One more thing I would like to know...how can I write down the boundary normal vector facing towards the interior? I want to have it like list of {nx,ny}, but after you make those lists there is no way to reproduce it. Doesn't ElementMesh have something like that? Aug 26, 2017 at 21:01
• @user16320 Have you seen the BoundaryNormals property? It is documented in the ElementMesh documentation page. Aug 26, 2017 at 22:21
• @user16320 Yes, I was going to say that BoundaryNormals always uses mesh order 1. Aug 27, 2017 at 1:13
• You could also use mesh = ToElementMesh[dom, "MeshOrder" -> 2]; coords = mesh["Coordinates"]; boundaryIDs = Flatten[ElementIncidents[mesh["PointElements"]]]; interiorIDs = Complement[Range[Length[coords]], boundaryIDs]; also note that "BoundaryNormals" is part of 11.1 and documented on the ElementMesh ref page. Aug 27, 2017 at 13:56

Look this post

region = DiscretizeRegion[dom];
Show[region, Graphics[{Red, PointSize[.02],
MeshPrimitives[DiscretizeRegion[dom], {0, "Interior"}]}]]


Using mesh region functions:

reg = DiscretizeRegion[dom];
boundary = RegionBoundary[reg];
coords = MeshCoordinates[reg];

pointindex =
Pick[Range[Length[coords]], RegionMember[boundary, coords], False];

MeshRegion[reg, MeshCellStyle -> {{0, pointindex} -> Red}]


If you want to get coordinates of interior points:

coords[[pointindex]]


or you can do from beginning:

Pick[coords, RegionMember[boundary, coords], False]