I think we are able to bring in the additional dof through the NumberSpaces, e.g.,

V = H1(mesh, order = 1, dirichlet = “Outer”)

N = NumberSpace(mesh)

X = FESpace([V,N])

I wonder however if there is a simpler way to setup the bilinear form for the circuit equation, which requires the division by the CoilArea:

a += tau/dt * A * I_ / CoilArea * dx(“Coil”)

Any suggestions ? I have attached the file in case anybody is interested (or finds any improvements / bugs).

-------- Adding file here as attaching is not working for my browser

from netgen.geom2d import CSG2d, Circle, Rectangle

CoilRadius = 0.01

circle = Circle( center=(0,0), radius=CoilRadius, mat=“Coil”, bc=“Boundary” )

rect = Rectangle(pmin=(-0.05,-0.05), pmax=(0.05,0.05), mat=“Air”, bc=“Outer” )

air = rect - circle

geo = CSG2d()

geo.Add(circle)

geo.Add(air)

from ngsolve import *

mesh = Mesh(geo.GenerateMesh(maxh=2.e-2, quad_dominated=True))

nu = 1./1.257e-6

U = 5.

R = 2.

Turns = 20

dt = 5.e-2

t = 0

Tend = 2.

##################

V = H1(mesh, order = 1, dirichlet = “Outer”)

N = NumberSpace(mesh)

X = FESpace([V,N])

(A,I), (A_,I_) = X.TnT()

CoilArea = Integrate(CoefficientFunction(1), mesh, definedon=mesh.Materials(“Coil”))

tau = 1. / CoilArea * Turns

a = BilinearForm(X, symmetric=False)

a += nu*grad(A)*grad(A_) * dx

a += -I * tau * A_ * dx(“Coil”)

a += tau/dt * A * I_ / CoilArea * dx(“Coil”)

a += R * I * I_ / CoilArea * dx(“Coil”)

c = Preconditioner(a, “direct”)

a.Assemble()

n = 0

Current = [ 0.]

time = [ 0.]

Aold = GridFunction(X)

while t < Tend:

```
f = LinearForm(X)
f += U * I_ / CoilArea * dx("Coil")
f += tau/dt * Aold.components[0] * I_ / CoilArea * dx("Coil")
S = GridFunction(X)
with TaskManager():
f.Assemble()
bvp = BVP(bf=a, lf=f, gf=S, pre=c)
bvp.Do()
Aold.vec.data = S.vec.data
t += dt
CoilCurrent = S.components[1].vec.data
n = n + 1
print(f"Time: {t} Current: {CoilCurrent}")
Current.extend(CoilCurrent)
time.extend([t])
Draw(S.components[0], mesh, "A")
```

import pylab

import numpy

pylab.plot(time, Current)

pylab.plot(time, (U/R) * numpy.ones_like(time))

pylab.xlabel(“Time”)

pylab.ylabel(“Current”)

pylab.ylim((0, 3))

pylab.savefig(“Current.png”)