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")