I’m attempting to implement a stray field solver within ngsolve, using ngbem to compute the double layer potential operator. The algorithm is from Redirecting, specifically Algorithm 12 =

I have some code that works for step (i) and (iii) using the standard methods of NGSolve. I am not convinced that step (ii) is functioning properly, but I am unsure of what the error is, or if my discrete formulation is correct.

Attached are four pieces of code:

- The
`strayfield.py`

file contains all of the actual code for the computations. - A test on a unit sphere with m=(0,0,1), for which the
`stray_field`

function should return (0,0,-1/3), and an energy of 2pi/9. - A test on a unit cube with m=(0,0,1). The stray field has no analytical form, but the energy should be exactly 1/6.
- An auxiliary file.

General_Functions.py (2.2 KB)

strayfield.py (8.4 KB)

strayfield_cube_test.ipynb (2.1 KB)

strayfield_sphere_test.ipynb (2.3 KB)

When I run these tests, I don’t see the convergence I would expect. It is “close”, so it cannot be too wrong. In the sphere test, we generally see convergence as h->0 (generating meshes with `maxh=1/n`

, for `n=1,2,4,...,32`

)

but it doesn’t seem to be fast enough.

|hmax | energy | L^2 error|

|1.986897527 | 0.6244154164 | 0.0006382110957|

|1.993270202 | 0.6343580189 | 0.0002835307397|

|0.6891751754 | 0.6780127155 | 2.15E-05|

|0.3123171932 | 0.6934280217 | 1.67E-06|

|0.1694168532 | 0.6969098422 | 7.87E-07|

|0.09040233825 | 0.6978044914 | 4.29E-07|

As for the cube test, it is even worse. Again, it gets closer to the expected value of 1/6, but eventually misses and goes past the target.

|hmax | energy|

|1.732050808 | 0.1350023107|

|1 | 0.1465266853|

|0.5097136556 | 0.1642313063|

|0.2726903446 | 0.1664041535|

|0.0688505752 | 0.1669602132|

|0.0372839163 | 0.1675898603|

Within the strayfield.py file, step (i) is dealt with in `compute_internal`

. We feed the output of this into `DL_boundary_projection`

to achieve step (ii). Then for step (iii), we use the gfu.Set method to put this onto the boundary of an `H1`

function and then apply the `BVP`

function. This is all combined in `stray_field`

.