Interior of discretized region

Question

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["NDSolve`FEM`"]
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.

C. E. · Accepted Answer · 2017-08-28T23:28:49.150

12

It can be done like this:

Needs["NDSolve`FEM`"]
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]
   }]
 ]

Mathematica graphics

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];

edited Aug 28 '17 at 23:28

answered Aug 26 '17 at 11:49

C. E.

70,533
6
140
264

Thank you, the second solution is actually very short and elegant! – user16320 Aug 26 '17 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? – user16320 Aug 26 '17 at 21:01
@user16320 Have you seen the BoundaryNormals property? It is documented in the ElementMesh documentation page. – C. E. Aug 26 '17 at 22:21
Thank you, however I get the error "There is no method BoundaryNormals for ElementMesh objects." even for the examples from the documentation page, so this option, I'm afraid, is of no use to me :( – user16320 Aug 26 '17 at 23:10
@user16320 The examples work for me in Mathematica 11.2. – C. E. Aug 26 '17 at 23:15
This is what it shows. – user16320 Aug 26 '17 at 23:25
@C.E. 11.2 is it a beta? You pretty lucky... ;) – PlatoManiac Aug 27 '17 at 00:24
@PlatoManiac Yes, I am lucky. Although it was unintentional to reveal this here, my point is simply that it works. – C. E. Aug 27 '17 at 00:39
Oh, I got it working. But I want to match the normal with the point on the boundary: by simply checking for size, First@bmesh is of a different size than the list of normals, thus it is not a one to one correspondence between points on the boundary and their corresponding normals (not to mention that the normals are oriented differently on the outer boundary than on the inner one - I expect all of them to be oriented either towards the domain, or the other way around, but not mixed). – user16320 Aug 27 '17 at 01:06
Never mind, I somehow managed it by setting MeshOrder back to 1, fiddling with MaxCellMeasure and adding SameTest to the Complement function you used. – user16320 Aug 27 '17 at 01:12
2

@user16320 Yes, I was going to say that BoundaryNormals always uses mesh order 1. – C. E. Aug 27 '17 at 01:13
3

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. – user21 Aug 27 '17 at 13:56

score 8 · Answer 2 · answered Aug 29 '17 at 01:22

8

Look this post

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

Mathematica graphics

answered Aug 29 '17 at 01:22

yode

26,686
4
62
167

score 6 · Answer 3 · answered Aug 28 '17 at 18:21

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]

Interior of discretized region

3 Answers3

Linked