# Orientation of boundary domain normals

Bug introduced in 11.2.0 or earlier, and fixed in 11.3.0

I would like to use BoundaryNormals provided by the ElementMesh object. To demonstrate my problem I will take my favorite domain:

dom = ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= (1/4)^2, {{x, 0, 1}, {y, 0, 1}}];
bmesh = ToBoundaryMesh[dom];
mesh = ToElementMesh[dom];
normals = bmesh["BoundaryNormals"];


I copied the example to show the orientation of normals:

mean = Mean /@ GetElementCoordinates[bmesh["Coordinates"], #] & /@
ElementIncidents[bmesh["BoundaryElements"]];
Show[ListPlot[First@mesh, AspectRatio -> 1, PlotRange -> {{-0.1, 1.1}, {-0.1, 1.1}}],
Arrow[{#1, #2}] &, {Join @@ mean, Join @@ (normals/15 + mean)}]]]


Now it's obvious, that something's off: outer boundary normals are oriented away from the domain, whereas inner (circle) normals are oriented toward the domain: I would like to have all normals oriented either towards the domain or away from domain, but not like this, inconsistently. This could be done manually by checking the dot product with the closest non-boundary point, however, I suspect it has something to do how Mathematica circles around the domain: if it's clockwise, then the said normal is oriented away, if counter-clockwise, it's oriented the other way around.

• Congrats on finding a solution to your problem! For the benefit of this site, definitely encourage posting your solution as an answer to your question. This way when people search they'll see this question as having an answer (accepted or not). Right now it looks like XX votes up, 0 answers in search results. – John Joseph M. Carrasco Aug 29 '17 at 6:43
• Okay, thanks, will do. – user16320 Aug 29 '17 at 10:48
• I filed this as a bug - all normals should point outward (opposite to what you unified them to in yout answer) – user21 Oct 17 '17 at 12:38
• This will be fixed in the next release (post V11.2) and the normal will point outward. – user21 Nov 12 '17 at 23:21
• That's really weird considering there are also Neumann BCs (besides Dirichlet which really doesn't need normal vectors). How it's the solution to the Neumann problem implemented internally if it doesn't utilize the normal vector(s)? :O – user16320 Nov 13 '17 at 7:56

I somehow found a solution. It's very straightforward and not a bit too elegant.

First we define several things to get started with (look here and here to understand details, in short: the interior and boundary coordinates need to be separate and normals must correspond to exact points on the boundary to proceed with the flipping):

Needs["NDSolveFEM"]
d = 0.0005;
dom = ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= (1/4)^2, {{x, 0,
1}, {y, 0, 1}}];
grid = ToElementMesh[dom, "MeshOrder" -> 1, MaxCellMeasure -> d][
"Coordinates"];
bmesh = ToBoundaryMesh[dom, "MeshOrder" -> 1, MaxCellMeasure -> d];
boundary = Partition[Flatten@bmesh["Coordinates"], 2];
normals = Partition[Flatten@bmesh["BoundaryNormals"], 2];
interior = Complement[grid, boundary, SameTest -> (Norm[#1 - #2] < d &)];
e = bmesh["BoundaryElements"];
elements =
Join @@ (GetElementCoordinates[boundary, #] & /@
ElementIncidents[e]);
n = Table[
Normalize[
Mean /@ Transpose@
normals[[First@
Transpose@Position[elements, boundary[[i]]]]]], {i, 1,
Length@boundary}];


Now the normals are correctly paired with the boundary coordinates, but they can point either toward or away from the domain. To fix that, just check the dot product with the closest interior point and if it's negative, flip the normal:

For[i = 1, i <= Length@n, i++,
v = First@Drop[Nearest[interior, boundary[[i]], 2], 1] -
boundary[[i]];
If[v.n[[i]] < 0, n[[i]] = -n[[i]]];
]


This produces the result I've been looking for: Feel free to clean up the code.

EDIT: user John Joseph M. Carrasco provided me with those so much needed improvements that are explained in this short article. I am very thankful to him for explaining how the syntax works. In the end, this:

n = Table[
Normalize[
Mean /@ Transpose@
normals[[First@
Transpose@Position[elements, boundary[[i]]]]]], {i, 1,
Length@boundary}];


can be replaced with:

Normalize@Mean@normals[[First /@ Position[elements, #]]] & /@ boundary;


and the flipping:

For[i = 1, i <= Length@n, i++,
v = First@Drop[Nearest[interior, boundary[[i]], 2], 1] -
boundary[[i]];
If[v.n[[i]] < 0, n[[i]] = -n[[i]]];
]


with a much shorter version:

n[[Select[Range@Length@n, n[[#]].((Nearest[interior, #, 2][] - #) &@boundary[[#]]) < 0 &]]] *= -1;

• just for fun: Instead of your table def'n of $n$, can use: n = Normalize@Mean@normals[[First /@ Position[elements, #]]] & /@ boundary; also, instead of using For loop to flip signs: n[[Select[Range@Length@n, n[[#]].((Nearest[interior, #, 2][] - #) &@boundary[[#]]) < 0 &]]] *= -1; – John Joseph M. Carrasco Aug 29 '17 at 11:52
• No, no it's not just for fun! This is EXACTLY the kind of corrections I need! What I wrote is not a proper Mathematica code, it's clumsy and c++-like. I will adjust the code so it's more elegant due to your suggestion. Thank you very much ^_^ Although...as I look at your code, as much as it works, I do not understand it :( I can understand what @, @@ and @@@ is, what /@ does, but there are too many strange characters... – user16320 Aug 29 '17 at 11:54
• I'll make some notes to walk from your expressions to the ones above and place a link here (relatively soon). – John Joseph M. Carrasco Aug 29 '17 at 13:42
• "At first you make it "c++ way" and then gradually shape it into Mathematica compressed form with @#&/, or you directly code in that way?" - it takes (lots of) practice (and unlearning habits obtained from using other languages), but you do eventually get used to programming functionally instead of programming procedurally. – J. M.'s discontentment Aug 29 '17 at 23:34
• 1) Re: "c++ vs mma" after refactoring code for efficiency you just start to internalize heuristic aesthetics around what mma likes and what grinds its gears. It's ok to sketch in your comfort language -- it will evolve. 2) I think refactoring questions and why-didn't-this-work (especially when you start with something that does work like your For loop) are great actual questions (not just comment Q's) for this site + you'll get a wide spectrum of answers. In this case you want something like: {grid, interior} = Join[#, boundary + offset n] & /@ {grid, interior}; – John Joseph M. Carrasco Aug 30 '17 at 5:18

This is fixed in Version 11.3:

Needs["NDSolveFEM"]
dom = ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= (1/4)^2, {{x, 0,
1}, {y, 0, 1}}];
bmesh = ToBoundaryMesh[dom];
mesh = ToElementMesh[dom];
normals = bmesh["BoundaryNormals"];
mean = Mean /@ GetElementCoordinates[bmesh["Coordinates"], #] & /@
ElementIncidents[bmesh["BoundaryElements"]];
Show[ListPlot[First@mesh, AspectRatio -> 1,
PlotRange -> {{-0.1, 1.1}, {-0.1, 1.1}}], 