Orientation of boundary domain normals

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:

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;

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)}]]]