Impose electric current magnitude / Floating potential

variational formulation
5

I think the problem only really makes sense if \Gamma_1 is a strict subset of \partial\Omega. Otherwise, only I_0 = 0 would be compatible with a zero source term, and the solution would just be u \equiv \text{const} in the entire domain. I hacked together an ad hoc Nitsche-like formulation for the problem on a unit square, with a nonzero Dirichlet BC on the left, insulating (J=0) Neumann BCs on the top and bottom, and this special BC on the right, in terms of an average current density,

J_\text{avg} := -I_0 / \vert\Gamma_1\vert\text{ ,}

where \vert\Gamma_1\vert is the length of the right edge. The formulation is: Find (u,u_0)\in S\times\mathbb{R} such that, \forall(v,v_0)\in V\times\mathbb{R},

\int_\Omega\sigma\nabla u\cdot\nabla v\,d\Omega\\ - \int_{\Gamma_1}v\sigma\nabla u\cdot\mathbf{n}\,d\Gamma + \int_{\Gamma_1}(u-u_0)\sigma\nabla v\cdot\mathbf{n}\,d\Gamma \\ + \int_{\Gamma_1}(J_\text{avg} + \sigma\nabla u\cdot\mathbf{n})v_0\,d\Gamma = 0\text{ ,}

where u_0 is the free constant that u must equal on \Gamma_1, S is the trial space satisfying the Dirichlet BC on the left, and V is the corresponding test space. The solution in the following FEniCS implementation appears to converge toward satisfaction of both u = u_0 on \Gamma_1 and \int_{\Gamma_1}\mathbf{J}\cdot\mathbf{n}\,d\Gamma = -I_0, as the mesh is refined:

from dolfin import *
N = 64
mesh = UnitSquareMesh(N,N)
cell = mesh.ufl_cell()
uE = FiniteElement("CG",cell,1)
u0E = FiniteElement("Real",cell,0)
W = FunctionSpace(mesh,MixedElement([uE,u0E]))
x = SpatialCoordinate(mesh)

# u is the potential, while u0 is a global scalar unknown, which is the
# constant potential over the right face of the domain.
u,u0 = TrialFunctions(W)
v,v0 = TestFunctions(W)

# Apply some arbitrary Dirichlet BC to the potential on the left face:
bc = DirichletBC(W.sub(0),
                 Expression("sin(3*pi*x[1])",degree=2),
                 "x[0]<DOLFIN_EPS")

# Characteristic function for right face:
right_char = conditional(gt(x[0],Constant(1.0-DOLFIN_EPS)),1.0,Constant(0))

# Formulation, using a Nitsche-like enforcement of u=u0 on the right
# face of the domain, and weak enforcement of an average current.
sigma = Constant(1)
def J(u):
    return -sigma*grad(u)
n = FacetNormal(mesh)
I_0 = 1.0
meas_Gamma1 = assemble(right_char*ds)
J_avg = Constant(-I_0/meas_Gamma1)
res = inner(-J(u),grad(v))*dx \
    + right_char*(dot(J(u),n)*v - dot(J(v),n)*(u-u0)
                  + (J_avg - dot(J(u),n))*v0)*ds
w = Function(W)
solve(lhs(res)==rhs(res),w,[bc,])
u_out,u0_out = w.split()

# Check that u=u0 on the right boundary:
e = u_out-u0_out
print(sqrt(assemble(right_char*e*e*ds)))

# Check that the total current through the right boundary
# has average value J_avg:
print(abs(assemble(right_char*(dot(J(u_out),n) - J_avg)*ds)))

# Visualize output:
u_out.rename("u","u")
File("u.pvd") << u_out