Sphere in a homogenous field

For an three axial ellipsoid in a constant field we have an analytic solution expressed in terms of demagnetization coefficients $n_i$:

\[H_i = \frac{H_{0i}}{(1 + (\mu-1)n_i)}\]

whre $\vec H$ is the interior field and $\vec H_{0}$ is the constant field. That makes ellipsoid a perfect subject to test the numerics.

using QuadGK

function ellipsoid_demagnetization_coefficients(a,b,c)

    UP_BOUND = 1000

    Ru2(u) = (u+a^2)*(u+b^2)*(u+c^2)

    nx = 1/2 * a*b*c * quadgk(s -> 1/(s+a^2)/sqrt(Ru2(s)), 0, UP_BOUND)[1]
    ny = 1/2 * a*b*c * quadgk(s -> 1/(s+b^2)/sqrt(Ru2(s)), 0, UP_BOUND)[1]
    nz = 1/2 * a*b*c * quadgk(s -> 1/(s+c^2)/sqrt(Ru2(s)), 0, UP_BOUND)[1]

    return [nx, ny, nz]
end

function EllipsoidField(a,b,c,mu,H0)

    H0x, H0y, H0z = H0
    nx, ny, nz = ellipsoid_demagnetization_coefficients(a,b,c)

    Hix = H0x/(1 + (mu-1)*nx)
    Hiy = H0y/(1 + (mu-1)*ny)
    Hiz = H0z/(1 + (mu-1)*nz)

    return [Hix,Hiy,Hiz]
end;

EllipsoidField (generic function with 1 method)

We define that we have a field of strength $H_0 = e_z$. A corresponding potential for that is $\psi_0 = z$, which we need to pass to the solver.

using LaplaceBIE
using LinearAlgebra

H0 = [0,0,1]

ψ0(x) = dot(x,H0)
∇ψ0(x) = H0;

∇ψ0 (generic function with 1 method)

For simplicity we consider a sphere with relative permitivity $\epsilon=15$. The sphere we can generate mesh from subdivisions of icosahedron given in sphere.jl. We also need normal vectors which we get from the vertex positions with high accuracy.

ϵ = 15

include("sphere.jl")
msh = unitsphere(3)
vertices, faces = msh.vertices, msh.faces
n = vertices;

642-element Array{GeometryTypes.Point{3,Float64},1}:
 [-0.5257311121191336, 0.85065080835204, 0.0]
 [0.5257311121191336, 0.85065080835204, 0.0]
 [-0.5257311121191336, -0.85065080835204, 0.0]
 [0.5257311121191336, -0.85065080835204, 0.0]
 [0.0, -0.5257311121191336, 0.85065080835204]
 [0.0, 0.5257311121191336, 0.85065080835204]
 [0.0, -0.5257311121191336, -0.85065080835204]
 [0.0, 0.5257311121191336, -0.85065080835204]
 [0.85065080835204, 0.0, -0.5257311121191336]
 [0.85065080835204, 0.0, 0.5257311121191336]
 ⋮
 [0.7802043707101266, 0.620239582613447, -0.08114185161993964]
 [0.7029070304877733, 0.7112817349622164, 0.0]
 [0.7802043707101266, 0.620239582613447, 0.08114185161993964]
 [0.8401778853271388, 0.5192584897281836, -0.15643446504023087]
 [0.8910065241883679, 0.38618738558759214, 0.2386769303195932]
 [0.912982492932299, 0.3996070517018512, 0.08232358003196016]
 [0.9876883405951378, 0.13307110414059134, -0.0822424652793623]
 [0.9638612634676226, 0.26640470113456743, 0.0]
 [0.912982492932299, 0.3996070517018512, -0.08232358003196016]

At the moment user needs to take the derivtive of the potential himself as there are so many ways to do that. It is also possible that user might wish to calculate the field due to the surface current for whic normalderivatives method could be useful.

Now we can proceed with calculation. To calculate surface potential everywhere on the surface one executes a method

ψ = surfacepotential(vertices,n,faces,ϵ,ψ0);

642-element Array{Float64,1}:
  1.8173340690005643e-17
  1.4137469307473297e-17
  1.5617380028648863e-17
  1.2021090333602364e-17
  0.1507079801429761
  0.15070798014297598
 -0.15070798014297596
 -0.1507079801429762
 -0.09314265410420328
  0.09314265410420361
  ⋮
 -0.01438351622337075
  1.4313488932057568e-17
  0.014383516223370809
 -0.027719944834048073
  0.042313386614868165
  0.014597289655777019
 -0.014593441780819992
  1.0398631304400936e-17
 -0.014597289655776974

which solves a regularized boundary integral equation with BLAS. That can then for example then can be used to calculate energy of the field.

Usually, however, one wants to know the field on the surface for force calculations due to $(M \cdot \nabla)H$. The library offers to use the calculated pontential for calculating tangential field components by a finite differences. That cna be easally achieved with method tangentderivatives:

P∇ψ = tangentderivatives(vertices,n,faces,ψ);

642-element Array{GeometryTypes.Point{3,Float64},1}:
 [-4.490946472340537e-17, -2.775557561562891e-17, 0.1773548873930608]
 [-1.263652073119082e-16, 7.809799311418602e-17, 0.17735488739306082]
 [-3.5427028040325176e-17, 2.1895107449316557e-17, 0.17735488739306066]
 [8.507771110527064e-17, 5.2580917148101646e-17, 0.17735488739306057]
 [6.132247487578013e-16, 0.07931551687054156, 0.04901968526126036]
 [-5.945316313203069e-16, -0.07931551687053988, 0.04901968526125933]
 [5.997554053931126e-16, -0.07931551687054056, 0.049019685261259736]
 [4.0592529337857286e-16, 0.0793155168705407, 0.04901968526125983]
 [0.07931551687054057, -2.1510571102112408e-16, 0.12833520213180033]
 [-0.07931551687054081, 7.720443752279785e-16, 0.1283352021318007]
 ⋮
 [0.011224197497333072, 0.00891038280492936, 0.17603418916161745]
 [4.9643519731222405e-17, -4.905901181208928e-17, 0.17725358302099356]
 [-0.011224197497332986, -0.008910382804929453, 0.17603418916161734]
 [0.02325660282856346, 0.01444222439731941, 0.17284510166164144]
 [-0.03764935424192934, -0.01635206746694934, 0.16700726958268122]
 [-0.013308919162044732, -0.005867920376971907, 0.1760816591143559]
 [0.014402891346747864, 0.0018832963793923255, 0.17601831527664324]
 [-2.090112744140055e-18, 7.562098948617905e-18, 0.17735078933190623]
 [0.013308919162044608, 0.005867920376972166, 0.17608165911435578]

And lastly to calculate the normal derivatives we can use Biot-Savarat law treating tangential field as raising from a surface current:

n∇ψ = normalderivatives(vertices,n,faces,P∇ψ,ϵ,∇ψ0);

642-element Array{Float64,1}:
  6.943952626461679e-17
 -4.3007029068577255e-17
 -1.803786798515189e-18
  2.452446043550493e-17
  0.1502038139798372
  0.1502038139798374
 -0.15020381397983693
 -0.15020381397983792
 -0.09283106227940602
  0.09283106227940634
  ⋮
 -0.014356536563115582
 -1.427595420399888e-17
  0.014356536563115672
 -0.02771546380312765
  0.04222118203520274
  0.01457484390063835
 -0.014505718232074822
 -6.506659544310088e-17
 -0.01457484390063827

To compare it is usefull to use the azimuthal symetry we have in the $z$ direction. Also the normal field is equal to potential since vertex positions coincide with normals. And so the comparisson:

Hin = EllipsoidField(1,1,1,ϵ,H0)

x = [v[3] for v in vertices]
sp = sortperm(x)

for xkey in sp[2:43:end]
    psit = dot(Hin,vertices[xkey])
    println("$(vertices[xkey][3]) \t $psit \t $(ψ[xkey]) \t $(n∇ψ[xkey])")
end

-0.9904388819568619 	 -0.17478787716673896 	 -0.1754781286019715 	 -0.17527596928348
-0.8626684804161862 	 -0.15223957291810866 	 -0.15283904735137047 	 -0.1525980317631358
-0.7071067811865475 	 -0.12478679448610892 	 -0.1252822510935816 	 -0.12509420513088595
-0.6015009550075456 	 -0.10614998760126844 	 -0.10658154929292299 	 -0.10634376822640942
-0.4539904997395468 	 -0.0801180538738164 	 -0.08043043818581969 	 -0.08024964268345774
-0.30901699437494745 	 -0.0545338288300314 	 -0.054749125808272105 	 -0.054667641762322305
-0.16245984811645317 	 -0.028670130478926876 	 -0.028809947025745908 	 -0.02873712859252811
-0.08108629344330352 	 -0.01430971800124824 	 -0.014381050442757757 	 -0.014330146478530993
0.0822424652793623 	 0.014513753630852578 	 0.014593441780820027 	 0.014505718232074869
0.2201170274732924 	 0.038845191420892995 	 0.038990634237237726 	 0.03895854971555256
0.35822879348657893 	 0.0632184897969791 	 0.06345968851389136 	 0.06327240001131587
0.48444164206066787 	 0.0854919246098835 	 0.08581995775709064 	 0.08562690297002987
0.6156420208736806 	 0.10864553470532934 	 0.10908898616739937 	 0.10881356703605721
0.7579354200477765 	 0.13375678753431802 	 0.1343014941269712 	 0.1340592320295373
0.8649293358632748 	 0.15263855778368335 	 0.15324065900681047 	 0.1528773858354282