Hello,
I am trying to solve for the current density vector of one-turn coil by specifying a source and sink, but the following Python code does not act as expected.
It is just a simple model but when the loop-coil shape is complicated, we would like to solve the laplace problem for the current densities.
Please advise us.
from ngsolve import *
from netgen.occ import *
p1 = Pnt(0.6, 0, -0.1)
p2 = Pnt(0.4, 0, -0.1)
p3 = Pnt(0.4, 0, 0.1)
p4 = Pnt(0.6, 0, 0.1)
rect1 = Face(Wire([Segment(p1, p2), Segment(p2, p3), Segment(p3, p4), Segment(p4, p1)]))
coil= Revolve(rect1, Axis(Pnt(0,0,0), Vec(0,0,1)), 360)
source1 = rect1.bc('plus')
rect2 = Face(Wire([Segment(p1, p2), Segment(p2, p3), Segment(p3, p4), Segment(p4, p1)]))
source2 = rect2.bc('minus')
coil = Glue([coil, source1, source2])
coil.maxh = 0.05
coil.mat('coil')
coil.faces[5].bc('minus')
coil.faces[0].bc('plus')
for n in range(len(coil.faces)):
print(coil.faces[n].name)
print(coil.faces[n].mass)
#for n in range(len(coil.solids)):
# print(coil.solids[n].name)
# print(coil.solids[n].mass)
#from netgen.webgui import Draw
#Draw(coil, clipping = { "pnt" : (0.0,-0.1,0), "vec" : (0,1,0) })
mesh = Mesh(OCCGeometry(coil).GenerateMesh(maxh=0.1)).Curve(2)
fes = H1(mesh, order=2, definedon="coil", dirichlet="minus")
phi,psi = fes.TnT()
a = BilinearForm(grad(phi)*grad(psi)*dx).Assemble()
inv = a.mat.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")
f = LinearForm(psi*ds("plus")).Assemble()
gfphi = GridFunction(fes)
gfphi.vec.data = inv * f.vec
J = -grad(gfphi)
#J = CoefficientFunction((1/(0.2*0.2)*y/sqrt(x**2+y**2), -1/(0.2*0.2)*x/sqrt(x**2+y**2), 0))
from ngsolve.webgui import Draw
Draw(J, mesh, "J", clipping = { "pnt" : (0.0,-0.1,0), "vec" : (0,1,0) })