I am solving the following finite element problem. When I try to take values from the solution I get Interpolation failures. 

    $Version
    "11.0.0 for Microsoft Windows (64-bit) (July 28, 2016)"

 Here is the finite element problem.

    Needs["NDSolve`FEM`"];
    Lp = 1; (* length of half plate *)
    tk = 0.5;(* plate thickness *)
    Ls = 9; (* from plate to boundary *)
    Lu = 5; (* upstream distance *)
    Ld = 5; (* downstream distance *)
    r1 = Rectangle[{-Lp, -Ld}, {Ls, Lu}];
    r2 = ImplicitRegion[x^2 + y^2 <= (tk/2)^2, {x, y}];
    r3 = Rectangle[{-Lp, -tk/2}, {0, tk/2}];
    rd = RegionDifference[r1, RegionUnion[r2, r3]];
    area1 = (tk/10)^2 Sqrt[3]/4;
    area2 = (Ls/10.)^2 Sqrt[3]/4;
    cf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
       Block[{r, r0 = tk/2},
        r = Norm@(Total[c]/3);
        If[a > area1 ( area2/area1)^((r - tk/2)/Ls), True, False]
        ]
       ];
    
    mesh = ToElementMesh[rd, MeshRefinementFunction -> cf];
    sol = NDSolveValue[{
        D[u[x, y], x, x] + D[u[x, y], y, y] ==
         NeumannValue[-1, -Lp <= x <= Ls && y == -Ld] + 
          NeumannValue[1, -Lp <= x <= Ls && y == Lu],
        DirichletCondition[u[x, y] == 0, x == -Lp && y == Lu]
        },
       u, {x, y} ∈ mesh];

I now generate some points where I want to calculate values. When I substitute the points into the solution I get, for example, "Input value {0.0369866,0.259881} lies outside the range of data in the \
interpolating function." With Indeterminate returned. 

    eps = 0.05;
    pts = Table[{(1 + eps) tk/2 Cos[(-s \[Pi])/2], (1 + eps) tk/
         2 Sin[(-s \[Pi])/2]}, {s, -1, 1, 0.001}];
    pp = sol[#[[1]], #[[2]]] & /@ pts;

I extract the Indeterminate points and check that they are within the mesh:

    ind = Flatten@Position[pp, Indeterminate];
    pp1 = pts[[ind]];
    Show[
     mesh["Wireframe"],
     mesh["Wireframe"["MeshElement" -> "PointElements", 
       "MeshElementStyle" -> Red]], PlotRange -> {{-tk, tk}, {-tk, tk}}, 
     Frame -> True, Epilog -> {Blue, Point[pp1]}]

![Mathematica graphics](https://i.sstatic.net/eaLqM.png)

The points where there are failurs are shown in blue. I have deliberately put in a value to take them off the boundary and into the mesh. My goal is the take the gradient in this region but it is too full of interpolation failures to work. I have tried different mesh sizes but that did not help. 

Is there a workaround? Thanks