I think you are right : it is very probable that there is a mistake somewhere.

At the moment I don't know where.

It is likely that my manipulations with Grad and coordinates transformations are not valid at the mathematical level (nothing to do with mathematica)

I did not find a mistake either.

Show[
 ContourPlot[solCyl[r, p], {r, p} \[Element] solCyl["ElementMesh"], ColorFunction -> "Temperature", Contours -> 20, PlotLegends -> Automatic],
 StreamPlot[functionGivingAPair[r, p], {r, 0, 1}, {p, 0, 2 Pi}, StreamStyle -> Black]]

If one has a look at the electric field lines before transformation one can see that they enter (and leave) the left boundary where r=0. After transformation this will be the center.

This is strange because a Neumann value of 0 means that the scalar product of the boundary normals and the gradient has to be zero. Thus, the electric field lines are only allowed to be parallel to the boundary.

I guess I found the problem. The coordinate transformation of a vector field is more involved than I anticipated.

This works:

electricField[r_, p_] = -Grad[solCyl[r, p], {r, p}, "Polar"];
electricField2[x_, y_] = TransformedField["Polar" -> "Cartesian", electricField[r, p + \[Pi]], {r, p} -> {x, y}];
Show[
 ContourPlot[solCyl[Sqrt[x^2 + y^2], ArcTan[x, y] + \[Pi]], {x, -1, 1}, {y, -1, 1}, ColorFunction -> "Temperature", Contours -> 20, PlotLegends -> Automatic],
 StreamPlot[electricField2[x, y], {x, -1, 1}, {y, -1, 1}, StreamStyle -> Black]]