Volume integration

Main Forum Old-Forum - do not post here
1

I’ve written a script (based on the examples) to get the electric field between conductors: a cylindrical rod and a surrounding larger cylindrical shield. I would like to integrate the square of the field over the volume mesh to give the stored energy, from which I can calculate capacitance. I’m interested to see how accurate the FEM value will be.
I integrate with:
func= E*E
integ = Integrate(func, mesh, VOL)
print (integ)
but this gives NaN.
I’m probably making an elementary mistake and would appreciate any help.
Andrew.

https://ngsolve.org/media/kunena/attachments/1012/proc_cyl.py
https://ngsolve.org/media/kunena/attachments/1012/coax.vol

Attachment: proc_cyl.py

Attachment: coax.vol

2

Hi Andrew,

your definition of epsl

epsl = CoefficientFunction( [1.0006] )

yields a singular matrix a as epsl is only defined on the first material domain. The sparsecholesky solver does not say if the system is singular, but integrating over the resulting u gives already NaN ( the UMFPACK solver complains that it is singular).

Defining

epsl = 1.0006

solves the problem. If you want to have different values of epsl on the domains you have to write

epsl = CoefficientFunction( [value first domain, value second domain, ...] )

Best
Michael

3

Hi Michael,

Thank you for the help. I’m very impressed by the speed and accuracy (the calculated capacitance compares well to results I’ve calculated with three other programs).

If I can trouble you further, I would like to make the adaptive meshing work. It is not actually needed for this problem, but it would help me to gain understanding. I’ve followed the recipe set out in the tutorial - but it does not quite work. See attachment (same .vol file needed).

Regards,

Andrew.

https://ngsolve.org/media/kunena/attachments/1012/proc_cyl3.py

Attachment: proc_cyl3.py

4

Hi Andrew,

you need the autoupdate=True flag also for die GridFunction u1

u1 = GridFunction(fes, name="extension", autoupdate=True)

Then, you also need to set the boundary data in every solve stage

def SolveBVP(): a.Assemble() f = u.vec.CreateVector() u1.Set(extendone, BND) f.data = a.mat * u1.vec u.vec.data -= a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * f Redraw (blocking=True)

Best
Michael

https://ngsolve.org/media/kunena/attachments/889/proc_cyl3.py

Attachment: proc_cyl3.py

5

Hi Michael,

I’ve edited the script as you suggested, but I’m still having problems with the refinement. I would appreciate your help again. The expected capacitance is about 2.326 pF for the dimensions given.

I get:
maxerr = 0.026982948684526564
maxerr = 0.9045386846912055
maxerr = 0.7627369747838466
maxerr = 0.9376416794932843
Capacitance = 24.875 pF

The revised script imports a module (meshcoaxcap.py) which generates the .vol file (line 14). If I comment out the line (#MakeCoaxMesh(…) ) and use last .vol file I get different output.
maxerr = 0.0269829486845266
maxerr = nan
maxerr = nan
etc.

That’s rather strange, because line 14 is at the top of the script and is called only once.

I note that the outer boundary shown in Netgen changes from red (1V) to blue (0V) after the first refinement. This suggests that there might be a problem with refinement on the conductor surfaces. Or perhaps some of the metadata associated with the mesh (surfnr etc.) is not correct?

I have 6.2.2001-0~ubuntu18.04.1 installed.

Regards,
Andrew.

https://ngsolve.org/media/kunena/attachments/1012/proc_cyl3.py
https://ngsolve.org/media/kunena/attachments/1012/meshcoaxcap.py

Attachment: proc_cyl3_2020-02-15.py

Attachment: meshcoaxcap.py

6

Hi andrew,

the problem with the different boundary condition after the first refinement was that you first set u to u1 and then set the right boundary condition for u1. You have to switch these lines to

def SolveBVP():
    a.Assemble()
    #first set u1
    u1.Set(voltage, definedon=mesh.Boundaries("outer"))# Boundary sources
    u.vec.data = u1.vec

I also observed different behavior when commented out the mesh generation function. I am not an expert of merging meshing, so I don’t really have a clue where this behavior comes from.

However, for this geometry, you can simplifying the mesh generation by only using CSG geometries

cube = OrthoBrick( Pnt(-outer_rad,-outer_rad,-leng/2), Pnt(outer_rad,outer_rad,2.0*quarterleng) ) cylb = Cylinder ( Pnt(0.0, 0.0, -leng), Pnt(0.0, 0.0, 0.0), outer_rad) cube2 = OrthoBrick( Pnt(-inner_rad,-inner_rad,-leng/4), Pnt(inner_rad,inner_rad,quarterleng) ) cyli = Cylinder ( Pnt(0.0, 0.0, -leng), Pnt(0.0, 0.0, 0.0), inner_rad) inner = (cyli*cube2).bc('inner') outer = (cylb*cube-inner).bc('outer') geo = CSGeometry() geo.Add(outer) #geo.Add(inner) mesh = Mesh(geo.GenerateMesh (maxh=h))

Is there a reason why you mesh also the inner part? The solution seems to be always zero there.
With this, you can now use strings to identify your boundaries and you can also Curve your mesh.

With the attached code I get the following output:
ndof = 10423
maxerr = 0.13384372761117735 , toterr = 2.267302450023096
ndof = 15180
maxerr = 0.027621924443630222 , toterr = 1.809503014134304
ndof = 26375
maxerr = 0.007188620406700055 , toterr = 1.3764282981704608
ndof = 52185
maxerr = 0.007488492893794685 , toterr = 0.8460445676985052
ndof = 53249
maxerr = 0.001808892090986848 , toterr = 0.8097602634799458
ndof = 83743
maxerr = 0.0005241508456320649 , toterr = 0.64365894648496
ndof = 141903
maxerr = 0.00018470303727374993 , toterr = 0.5030248548701931
ndof = 245948
maxerr = 9.058668270776927e-05 , toterr = 0.37027924480707936
ndof = 326343
maxerr = 4.618238944138729e-05 , toterr = 0.30930491094782586
ndof = 418213
maxerr = 2.0765179180659515e-05 , toterr = 0.2665229929125358
ndof = 594250
maxerr = 1.1146830833658861e-05 , toterr = 0.22296841251806954
Capacitance = 2.3274 pF

which is quite near to your expected capacitance of 2.326 pF

Best,
Michael

https://ngsolve.org/media/kunena/attachments/889/proc_cyl3_2020-02-17.py

Attachment: proc_cyl3_2020-02-17.py