Hey
I want to solve the Laplace equation with Neumann boundary conditions on all boundaries, The solution seems random when I do not include any DirichletConditions though.
Here is an example of the Laplace in cylindrical coordinates (with cylindrical symmetry). "mesh" and "region" are define on the region 0 <= r <= 20 and 0 <= z <= 30.
T = NDSolveValue[
op == NeumannValue[r, z == 0] + NeumannValue[r, z == 30],
u, {r, z} \[Element] mesh,
Method -> {"FiniteElement",
"MeshOptions" -> {"BoundaryMeshGenerator" -> "Continuation"}}];
Show[
ContourPlot[T[r, z], {r, z} \[Element] region , Mesh -> None,
ColorFunction -> "TemperatureMap", Contours -> 40,
AspectRatio -> Automatic, PlotRange -> All,
Method -> {"GridLinesInFront" -> True}]
]
"op" is: 
And I get the result: 
What do I do wrong?