# How to find the direction of the outward normal in a point of the boundary of a 2D region?

Suppose we have a domain $\Omega$ such that RegionQ[\[Omega]] is True and RegionDimension[[\Omega]] == RegionEmbeddingDimension[\[Omega]] == 2 is True.

We also have a point $P$ such that RegionMember[RegionBoundary[\[Omega]], P] is True.

There is a way to find the direction of the outward normal to $\Omega$ in $P$?

I know how to mathematically obtain the outward normal if the region is described by a parametric boundary for example, but I don't know how to compute the outward normal for a generic Region.

The method should be reasonably efficient because I need to apply to many points lying on the boundary.

A different approach: Recall that the gradient of a potential is normal to equipotential surfaces. So if we solve a heat/voltage equation with internal sink/charge (e.g. $\nabla^2 u = 1$) and a fixed-temperature/conducting boundary (e.g. $u = 0$), the gradient of the solution will be an outward pointing normal.

(* Given *)
reg = RegionUnion[Disk[{0, 0}, 5/4], Rectangle[]];

(* Solve heat equation and take gradient *)
<< NDSolveFEM
mesh = ToElementMesh[reg, "MaxBoundaryCellMeasure" -> 0.01];
pot = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 1,
DirichletCondition[u[x, y] == 0, True]}, u, Element[{x, y}, mesh]];
norm = Grad[pot[x, y], {x, y}] // Normalize;

(* Some boundary points *)
sols = Solve[RegionMember[RegionBoundary@reg, {x, y}] && (x == 1/4 || y == 7/8), {x, y}];

(* Points and normals *)
bpts = {x, y} /. sols;
bnorms = norm /. sols;

(* Hooray *)
Show[
RegionPlot[reg],
Graphics@{MapThread[Arrow[{##}] &, {bpts, bpts + bnorms/2}]},
PlotRange -> All, AspectRatio -> Automatic
]


This works for any region, including ones with holes. Second example region:

reg = ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= (1/3)^2, {{x, 0, 1}, {y, 0, 1}}]; • +1, this is a neat solution. There is also the possibility to use the BoundaryNormal property of the mesh as is done here. This should now also work for holes, but only for Mathematica 11.2 and higher. – C. E. Mar 6 '18 at 8:37
• Thanks @C.E. I came across all of the related boundary normal questions earlier this week, but wasn't satisfied with BoundaryNormal because I have to wait for v11.2 and also it would mean writing a Nearest or similar thing to map points to actual boundary mesh midpoints, which isn't as clean and satisfying as a "continuum-style" solution. – user40265 Mar 6 '18 at 14:35
• Very nice! ${}$ – J. M.'s discontentment Mar 8 '18 at 0:05

Maybe this can get you started. It assumes the point {x, y} is on the boudary and that the region is described by inequalities involving either <= or >= that are returned by RegionMember.

normal[reg_] := Piecewise /@ Transpose[
Thread[{Normalize@D[First@#, {{x, y}}], # /. LessEqual -> Equal}] & /@
Cases[
Simplify[
RegionMember[reg, {x, y}], (x | y) ∈ Reals] /.
{a_ <= v_ <= b_ :> a - v <= 0 && v - b <= 0,
a_ <= b_ :> a - b <= 0, a_ >= b_ :> b - a <= 0},
_LessEqual,
Infinity]]

reg = RegionUnion[Disk[{0, 0}, 5/4], Rectangle[]];
sols1 = Solve[RegionMember[RegionBoundary@reg, {x, y}] && x == -1/2, {x, y}];
sols2 = Solve[RegionMember[RegionBoundary@reg, {x, y}] && y == 7/8, {x, y}];

Show[
RegionPlot[reg],
Graphics[{
Arrow /@ ({{x, y}, {x, y} + normal[reg]} /. Join[sols1, sols2])
}],
PlotRange -> All, AspectRatio -> Automatic
] • This is a good starting point, thanks, but I have some problem applying the method 1) to "numerical" regions with points "approximately" on the boundary of that region and 2) to region with holes. I'll come back on this problem within few days... – unlikely Dec 16 '14 at 21:06
• @unlikely Anticipating the indefinite variety of possible region-member conditions and handling them all seems like a big project. E.g. I didn't include all the possible inequalities. Discretized regions can probably be handled by using Cross on the boundary elements. – Michael E2 Dec 16 '14 at 23:34

I'm not sure it will be useful, but one approach is to flip the problem on its head by imagining a ball about the region r, and finding the nearest points from the ball to the region... these are perpendiculars, though you might need to search them to find a specific one:

r = Disk[{0, 0}, {2, 1}];
pts = Table[3 {Cos[k 2 \[Pi]/16], Sin[k 2 \[Pi]/16]}, {k, 0., 15}];
nst = RegionNearest[r, #] & /@ pts;
Graphics[{{Gray, r}, {Thin, Line[Thread[{pts, nst}]]}, Point[pts], Point[nst]}] • Should work if the region is convex. Might not be fast. Probably not effective for nonconvex regions though. – Michael E2 Dec 16 '14 at 1:18
• @Michael E2 -- Agreed. It would be nice to have a RegionNormal function that gives the perpendicular to any region at a point! – bill s Dec 16 '14 at 1:19