Neumann boundary conditions in NDSolve over nontrivial region

Here is what I think the issue is: Let's look at what NDSolve parses.

Needs["NDSolve`FEM`"]
{state} = 
  NDSolve`ProcessEquations[{op == 0, bc}, u, {x, y} ∈ reg, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {MaxCellMeasure -> 0.005}}];
femData = state["FiniteElementData"];
femData["PDECoefficientData"]["All"]

{{{{0}}, {{{{0}, {0}}}}}, {{{{{-1, 
      0}, {0, -1}}}}, {{{{0}, {0}}}}, {{{{2 E^(-x^2 - y^2) x, 
      2 E^(-x^2 - y^2) y}}}}, {{4 E^(-x^2 - y^2) - 
     4 E^(-x^2 - y^2) x^2 - 4 E^(-x^2 - y^2) y^2}}}, {{{0}}}, {{{0}}}}

So, there is the diffusion term ({{-1,0},{0,-1}}) the convection term ({2 E^(-x^2 - y^2) x, 2 E^(-x^2 - y^2) y}) and a reaction term. This means NDSolve is modeling something like:

∇⃗ ⋅[∇⃗ u(x,y)] + v(x,y)∇⃗ u(x,y) + u(x,y)∇⃗ v(x,y) = 0

If we manually set the PDE coefficients like this, to be a diffusion coefficient and a conservative convection coefficient (with an added minus sign, as noted in the comments)

methodData = femData["FEMMethodData"];
bcData = femData["BoundaryConditionData"];
sd = state["SolutionData"][[1]];
vd = methodData["VariableData"];

pdeData = 
 InitializePDECoefficients[vd, sd, 
  "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, 
  "ConservativeConvectionCoefficients" -> {{-{2 E^(-x^2 - y^2) x, 
      2 E^(-x^2 - y^2) y}}}];

This should be more like

∇⃗ ⋅[∇⃗ u(x,y) + u(x,y)∇⃗ v(x,y)]= 0

And discretize and solve the PDE with:

dpde = DiscretizePDE[pdeData, methodData, sd];
dbcs = DiscretizeBoundaryConditions[bcData, methodData, sd];
{l, s, d, m} = dpde["SystemMatrices"];
DeployBoundaryConditions[{l, s, d, m}, dbcs]
lsol = LinearSolve[s, l];
mesh = NDSolve`SolutionDataComponent[sd, "Space"]["ElementMesh"];
sol = {u -> ElementMeshInterpolation[{mesh}, lsol]};

We see that the difference between the in and out flux is relatively small.

Plot[{Evaluate[((Derivative[0, 1][u][x, 5] + 
        u[x, 5] Derivative[0, 1][v][x, 5]) - (Derivative[0, 1][u][
         x, -5] + u[x, -5] Derivative[0, 1][v][x, -5])) /. 
    sol]}, {x, -5, 5}]

Update

In 10.0.2 a better way to input the PDE and get the expected result is made available:

reg = 
  ImplicitRegion[-5 <= x <= 5 && -5 <= y <= 5 && x^2 + y^2 >= 1^2, {x,
     y}];
v = Function[{x, y}, -E^(-x^2 - y^2)];
bc = {DirichletCondition[u[x, y] == 1, y == -5], 
   DirichletCondition[u[x, y] == 0, y == 5]};
α = Grad[v[x, y], {x, y}];
op = Inactive[Div][
   Inactive[Plus][Inactive[Grad][u[x, y], {x, y}], 
    Inactive[Times][α, u[x, y]]], {x, y}];
solI = NDSolve[{op == 0, bc}, u, {x, y} ∈ reg];
Plot[Evaluate[Function[y,
    (Derivative[0, 1][u][x, y] + 
       u[x, y] Derivative[0, 1][v][x, y]) /. 
     solI] /@ {5, -5}], {x, -5, 5}, PlotRange -> {-0.2, 0}]

I would say that you misread the documentation for NeumannValue. I found the documentation fairly difficult to fathom. This was mainly my fault because what they do is substantially different than what I expected. I was trying to cram the round peg that is NeumannValue into the square hole that is my brain.

If you can write the PDE as

Div J = 0

(the form your PDE is in) and you use NeumannValue with its first argument set to zero you are then getting n.J=0. That is if you write your pde as

 Div J = NeumannValue[q - r u,pred]

where u is the dependent variable you are solving for and where pred is true for coordinates that are in the boundaries and false for all other coordinate values then NeumannValue is imposing n.J = q-ru on the boundaries for which pred is true. Thus if you set q =r=0 so that the first argument to NeumannValue is zero then you would in fact be getting the boundary conditions you wanted. If your PDE can't be written as the divergence of a flux (plus time derivatives possibly) or if you want some other boundary conditions then you can't use NeumannValue as far as I know.