I have a 3-dimensional domain D, with 3 types of BC which I am trying to solve the Poisson equation on. $$-\nabla \cdot (\sigma(x,y,z)\nabla u)=0$$

Insulator (Neumann BC)

Electrode set at some potential (Dirichlet BC).

Additionally I have some electrodes $E_l$ which are set to have a fixed current $I_l$, so on these the boundary condition is $\int_{E_l} \sigma(x,y,z)(du/dn) dA=I_l$.

How would I go about incoprorating this third type of BC into FEM? Any description or references is what im looking for.

Additionaly or alternatively, I would like to know if there is some way to obtain accurate flux computations over surfaces in FEM. The normal way to solve Poisson equation using FEM does not give very accurate fluxes, however commercial software using FEM does, does anyone know how they do it? I am looking for references on how to compute such accurate fluxes.