A sphere in a homogenous field

For a 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 permittivity $\epsilon=15$. The mesh for the sphere we generate 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.525731, 0.850651, 0.0]      
 [0.525731, 0.850651, 0.0]       
 [-0.525731, -0.850651, 0.0]     
 [0.525731, -0.850651, 0.0]      
 [0.0, -0.525731, 0.850651]      
 [0.0, 0.525731, 0.850651]       
 [0.0, -0.525731, -0.850651]     
 [0.0, 0.525731, -0.850651]      
 [0.850651, 0.0, -0.525731]      
 [0.850651, 0.0, 0.525731]       
 ⋮                               
 [0.780204, 0.62024, -0.0811419] 
 [0.702907, 0.711282, 0.0]       
 [0.780204, 0.62024, 0.0811419]  
 [0.840178, 0.519258, -0.156434] 
 [0.891007, 0.386187, 0.238677]  
 [0.912982, 0.399607, 0.0823236] 
 [0.987688, 0.133071, -0.0822425]
 [0.963861, 0.266405, 0.0]       
 [0.912982, 0.399607, -0.0823236]

At the moment users need to take the derivative of the potential themselves which can be done in multiple ways. It is also possible that a user might wish to calculate the field due to the surface current for which 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}:
  5.7758016744828055e-18
  4.910019280439825e-18 
  1.664433787138662e-18 
  2.1086403756558402e-18
  0.15070798014297607   
  0.15070798014297587   
 -0.15070798014297612   
 -0.15070798014297615   
 -0.09314265410420339   
  0.09314265410420361   
  ⋮                     
 -0.014383516223370783  
  9.687533947885387e-18 
  0.01438351622337079   
 -0.027719944834048084  
  0.04231338661486813   
  0.014597289655776991  
 -0.01459344178082      
  1.399702275284999e-17 
 -0.01459728965577698

which solves a regularized boundary integral equation with BLAS. That then, for example, can be used to calculate the 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 do a finite differentiation on the calculated potential for calculating tangential field components. That can be easily achieved with method tangentderivatives:

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

642-element Array{GeometryTypes.Point{3,Float64},1}:
 [-1.25456e-16, -7.75363e-17, 0.177355]
 [-1.17003e-16, 7.23119e-17, 0.177355] 
 [-6.72952e-17, 4.15907e-17, 0.177355] 
 [3.56034e-17, 2.20041e-17, 0.177355]  
 [8.94127e-16, 0.0793155, 0.0490197]   
 [1.66959e-17, -0.0793155, 0.0490197]  
 [2.00769e-16, -0.0793155, 0.0490197]  
 [1.29338e-16, 0.0793155, 0.0490197]   
 [0.0793155, 9.59317e-18, 0.128335]    
 [-0.0793155, 4.04283e-16, 0.128335]   
 ⋮                                     
 [0.0112242, 0.00891038, 0.176034]     
 [1.38199e-17, -1.36571e-17, 0.177254] 
 [-0.0112242, -0.00891038, 0.176034]   
 [0.0232566, 0.0144422, 0.172845]      
 [-0.0376494, -0.0163521, 0.167007]    
 [-0.0133089, -0.00586792, 0.176082]   
 [0.0144029, 0.0018833, 0.176018]      
 [4.13775e-18, -1.49705e-17, 0.177351] 
 [0.0133089, 0.00586792, 0.176082]

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

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

642-element Array{Float64,1}:
  4.3478609925839566e-17
 -5.661778784752401e-18 
  2.1903125410541584e-17
  7.549456373568602e-18 
  0.15020381397983737   
  0.15020381397983681   
 -0.1502038139798374    
 -0.1502038139798376    
 -0.09283106227940627   
  0.09283106227940639   
  ⋮                     
 -0.014356536563115673  
 -3.663712594986743e-18 
  0.014356536563115663  
 -0.027715463803127716  
  0.04222118203520266   
  0.014574843900638273  
 -0.014505718232074822  
 -5.152964225180785e-18 
 -0.014574843900638228

We can compare numerics with analytics easily by the use of azimuthal symmetry in the z-direction. The normal field is equal to potential since vertex positions coincide with normals which is thus not shown twice. And so the comparison:

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.17547812860197165 	 -0.17527596928348
-0.8626684804161862 	 -0.15223957291810866 	 -0.1528390473513706 	 -0.15259803176313602
-0.7071067811865475 	 -0.12478679448610892 	 -0.12528225109358154 	 -0.12509420513088565
-0.6015009550075456 	 -0.10614998760126844 	 -0.10658154929292295 	 -0.1063437682264094
-0.4539904997395468 	 -0.0801180538738164 	 -0.08043043818581963 	 -0.08024964268345759
-0.30901699437494745 	 -0.0545338288300314 	 -0.05474912580827211 	 -0.05466764176232228
-0.16245984811645317 	 -0.028670130478926876 	 -0.0288099470257459 	 -0.028737128592528013
-0.08108629344330352 	 -0.01430971800124824 	 -0.014381050442757769 	 -0.014330146478530944
0.0822424652793623 	 0.014513753630852578 	 0.01459344178082001 	 0.014505718232074843
0.2201170274732924 	 0.038845191420892995 	 0.03899063423723775 	 0.038958549715552626
0.35822879348657893 	 0.0632184897969791 	 0.06345968851389137 	 0.06327240001131593
0.48444164206066787 	 0.0854919246098835 	 0.08581995775709064 	 0.08562690297002988
0.6156420208736806 	 0.10864553470532934 	 0.10908898616739926 	 0.10881356703605717
0.7579354200477765 	 0.13375678753431802 	 0.13430149412697107 	 0.13405923202953704
0.8649293358632748 	 0.15263855778368335 	 0.15324065900681055 	 0.15287738583542854