Circulant integral on cut surface

General
1

I would like to compute
I = \oint {{\bf H} \cdot d{\bf l}}
on HCurl elements.
Our current implementation is shown below. It may not work for some shapes, is there a better way to implement it?

def integrate(mesh, gf):
    sum=0.
    it=mesh.BBoundaries("inedge").Elements()
    for e in it:
#        print(e.vertices)
        start=e.vertices[0].nr
        end=e.vertices[1].nr
        dir=1
        if start> end: dir=-1
        for i in range(len(e.edges)): 
            k=gf.space.GetDofNrs(e.edges[i])
            sum = sum+gf.vec[k[0]]*dir
#            print(gf.vec[k[0]])

    return sum;
2 
Integrate(gf*dir, mesh.BBoundaries("inedge"))

should work.

Best

1 Like
3

Thank you. The following script seems to be working except for ths signs.
How do I control the signs of normals and tangentials?

from ngsolve import *
from netgen.occ import *
from netgen.webgui import Draw

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)

rect = Wire([Segment(p1, p2), Segment(p2, p3), Segment(p3, p4), Segment(p4, p1)])
coil= Revolve(rect, Axis(Pnt(0,0,0), Vec(0,0,1)), 360)
source = Face(rect).bc('source')
source.edges[0].name = 'edges'
source.edges[1].name = 'edges'
source.edges[2].name = 'edges'
source.edges[3].name = 'edges'
source.maxh = 0.01
coil = Glue([coil, source])
coil.maxh = 0.05
sphere = Sphere(Pnt(0, 0, 0), 1.0).bc('outer')
shape = Glue([sphere,coil])
shape.solids[0].mat('air')
shape.solids[0].faces.col = (1,0,0)
shape.solids[1].mat('coil')
shape.solids[1].faces.col = (0,0,1)
#Draw(shape, clipping = { "pnt" : (0.0,-0.1,0), "vec" : (0,1,0) })
mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.1))
mesh.Curve(3)

print(mesh.GetMaterials())
print(mesh.GetBoundaries())
print(mesh.GetBBoundaries())

fes = VectorH1(mesh, order=2)
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))

normal = specialcf.normal(3)
print(f'Surface Integral = {Integrate(InnerProduct(normal,J), mesh, definedon=mesh.Boundaries("source")):.16f}A')

fes = HCurl(mesh, order=3, dirichlet="outer", nograds = True)
A,w = fes.TnT()

a = BilinearForm(fes, symmetric=True)
a += curl(A)*curl(w)*dx + 1e-6*A*w*dx

c = Preconditioner(a, type="bddc")
f = LinearForm(fes)
f += J * w * dx("coil")

gfA= GridFunction(fes)

with TaskManager():
    a.Assemble()
    f.Assemble()
    solver = CGSolver(mat=a.mat, pre=c.mat)
    gfA.vec.data = solver * f.vec

H= GridFunction(fes)
H.Set(curl(gfA))
tangent = specialcf.tangential(3)
print(f'Line integral = {Integrate(InnerProduct(tangent, H), mesh.BBoundaries("edges")):.16f}A')

from ngsolve.webgui import Draw
Draw (H, mesh, clipping = { "pnt" : (0.0,-0.1,0), "vec" : (0,1,0) })