9
$\begingroup$

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}}], 
Graphics[MapThread[
   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:

enter image description here

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.

Improve this question
$\endgroup$
10
  • 4
    $\begingroup$ 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. $\endgroup$
    – John Joseph M. Carrasco
    Aug 29, 2017 at 6:43
  • $\begingroup$ Okay, thanks, will do. $\endgroup$
    – user16320
    Aug 29, 2017 at 10:48
  • 2
    $\begingroup$ I filed this as a bug - all normals should point outward (opposite to what you unified them to in yout answer) $\endgroup$
    – user21
    Oct 17, 2017 at 12:38
  • $\begingroup$ This will be fixed in the next release (post V11.2) and the normal will point outward. $\endgroup$
    – user21
    Nov 12, 2017 at 23:21
  • 1
    $\begingroup$ 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 $\endgroup$
    – user16320
    Nov 13, 2017 at 7:56

2 Answers 2

Reset to default
6
$\begingroup$

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

enter image description here

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][[2]] - #) &@boundary[[#]]) < 0 &]]] *= -1;
Improve this answer
$\endgroup$
13
  • 2
    $\begingroup$ 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][[2]] - #) &@boundary[[#]]) < 0 &]]] *= -1; $\endgroup$
    – John Joseph M. Carrasco
    Aug 29, 2017 at 11:52
  • 1
    $\begingroup$ 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... $\endgroup$
    – user16320
    Aug 29, 2017 at 11:54
  • 2
    $\begingroup$ I'll make some notes to walk from your expressions to the ones above and place a link here (relatively soon). $\endgroup$
    – John Joseph M. Carrasco
    Aug 29, 2017 at 13:42
  • 2
    $\begingroup$ "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. $\endgroup$
    – J. M.'s missing motivation ♦
    Aug 29, 2017 at 23:34
  • 2
    $\begingroup$ 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}; $\endgroup$
    – John Joseph M. Carrasco
    Aug 30, 2017 at 5:18
4
$\begingroup$

This is fixed in Version 11.3:

Needs["NDSolve`FEM`"]
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}}], 
 Graphics[MapThread[
   Arrow[{#1, #2}] &, {Join @@ mean, Join @@ (normals/15 + mean)}]]]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Keep up the good job! $\endgroup$
    – user16320
    Mar 9, 2018 at 21:04
  • $\begingroup$ @user16320, thanks for the encouragement! $\endgroup$
    – user21
    Mar 12, 2018 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.