This is a bug that only comes into play with the refinement function. Let me try to explain what happens and then offer workarounds.

Let's look at a particular example of a query point (qp) and a (blue) boundary mesh element and a neighboring interior element.

qp = {0.03821081952950365`, 0.25970403013985693`};
m = sol["ElementMesh"];
c = m["Coordinates"];
eleNr = 57;
(Join @@ m["ElementConnectivity"])[[eleNr]]
ele1 = (Join @@ ElementIncidents[m["MeshElements"]])[[eleNr]];
ele2 = (Join @@ ElementIncidents[m["MeshElements"]])[[497]];
uele = Union[ele1, ele2];
Graphics[{
  {Red, Point[qp]},
  {Blue, Line[c[[ele1]][[{1, 4, 2, 5, 3, 6, 1}]]]}, 
  Line[c[[ele2]][[{1, 4, 2, 5, 3, 6, 1}]]],
  MapThread[Text, {uele, c[[uele]]}]
  }]

[![enter image description here][1]][1]

What we see here is that the blue element (which is at the boundary) is (legitimately) curved (at node 8209) but also the black interior element is curved (at node 8210) which is not allowed. The interpolation algorithm tests if qp is in the blue element and concludes that it's just not inside. The picture here is not ideal as the black line with a kink is actually a smooth curve. In any case the interpolation code then tries the neighboring element. For that element, however, there is an optimization that the element can not be curved (as it's not an element on the region boundary) and uses a much cheaper linear inclusion test, which also fails. So the interpolation concludes that this point is outside region which is not correct.

One workaround you could experiment with is to set:

SetSystemOptions[
 "FiniteElementOptions" -> 
  "DefaultExtrapolationHandler" -> {Automatic, 
    "WarningMessage" -> False}]

How did this interior curved element get into the mesh in the first place? This is probably best explained by another workaround:

nr = ToNumericalRegion[rd];
mesh = ToElementMesh[rd, MeshRefinementFunction -> cf, 
   "MeshOrder" -> 1];
m2 = ImproveBoundaryPosition[mesh, nr];
m3 = MeshOrderAlteration[m2, 2];
m4 = ImproveBoundaryPosition[m3, nr];
(*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} \[Element] m4];*)

The point to realize in the above picture is that it's not the mid-side node 8210 on the black line that causes problems but the movement of the blue element corner node 296! The idea for a workaround is to generate a first order mesh and improve the boundary points. Next, we make a second order mesh. This means that for interior elements the introduced the mid side notes will be at the right position. After the second order mesh generation we improve again. In the second improvement case the element corner nodes do not move (since they have been fixed in the first round) and only the mid side nodes of the boundary elements are moved.

Other things you can try:

mesh = ToElementMesh[rd, MeshRefinementFunction -> cf, 
  "MaxBoundaryCellMeasure" -> 0.05]

This works because the boundary mesh is fine enough that the problem does not come up in the first place. Giving the refinement function alone does not pre-partition the boundary mesh (which maybe it should).

This also works (at the cost of may elements):

mesh = ToElementMesh[rd, MaxCellMeasure -> 0.001]

For the same reason as the previous example does.

If a first order mesh is fine for you than that's also an option. [1]: https://i.sstatic.net/VwQUf.png