Lost with boundary conditions in ME spaces

I’ve given more thoughts about this observation. My current conclusion is that the claim that vanishing Neumann boundary conditions are applied by default if no boundary conditions are explicitly specified, for a given finite element function, is not correct. And the above MWE is indeed a case where this happens.

If you go through the hand derivation of the weak form of the PDE \nabla \cdot \vec J (equivalently, if you integrate this equation over the whole volume and then apply Gauss theorem), you’ll find a different expression than the one I posted in my MWE. In fact, in my MWE, it is exactly the same problem except that I enforce a given electrical current. By doing so, \nabla V has to adjust such that the input and output currents match what I prescribe on the boundaries. And \nabla V does not vanish on all boundaries, in particular those where Dirichlet boundary conditions are specified for the other unknown (temperature in this case).

This also responds to Implementing a potential problem with a current condition - #21 by Daniel.Baker. I do not believe prescribing the current on the 2 boundaries (in and out current) is insufficient to fully specify the voltage (scalar field). I believe it is perfectly fine, keeping in mind that voltage (or electrostatics potential) is defined up to a constant.

I would appreciate if @kamensky could comment.