Mixed Neumann-Dirichlet boundary conditions in 3d

Natural boundary conditions vanish in your variational formulation. The following works as expected.

problem Laplace(u, v, solver=CG)
= int3d(Th)(Grad(u)' * Grad(v))
+ on(1, u=20.0)
+ on(5, u=-13.1);