# Impose electric current magnitude / Floating potential

**URL:** <https://fenicsproject.discourse.group/t/impose-electric-current-magnitude-floating-potential/8022>\
**Category:** variational formulation\
**Created:** [April 11, 2022, 12:13pm UTC](https://fenicsproject.discourse.group/t/impose-electric-current-magnitude-floating-potential/8022 "2022-04-11T12:13:43Z")\
**Posts on this page:** 1\
**Showing post:** 5

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [April 21, 2022, 11:01pm UTC](https://fenicsproject.discourse.group/t/impose-electric-current-magnitude-floating-potential/8022/5 "2022-04-21T23:01:37Z")

</div>

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:

```python
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

```

---

_[View the full topic](https://fenicsproject.discourse.group/t/impose-electric-current-magnitude-floating-potential/8022)._