Maxwell solver with mixed BC

Maxwell solver with mixed BC#

keys: mixed bvp, MoM, Maxwell equations, Dirichlet data, Neumann data, manufactured solution

from netgen.occ import *
import netgen.meshing as meshing
from ngsolve import *
from ngsolve.webgui import Draw
from libbem import *
from ngsolve import Projector, Preconditioner
from ngsolve.krylovspace import CG, GMRes

Loading ngbem library

Consider the mixed problem:


\( \begin{array}{l rcl r} & \mathrm{curl}\,\mathrm{curl} \boldsymbol E - \kappa^2 \boldsymbol E &=& \boldsymbol 0 &\mathrm{in}\; \Omega^c \,,\\ \textnormal{Dirichlet condition} & \gamma_R \boldsymbol E &=& \boldsymbol m & \mathrm{on}\; \Gamma_0\,,\\ \textnormal{Neumann condition} & \gamma_N \boldsymbol E &=& \boldsymbol j & \mathrm{on}\; \Gamma_1\,, \\ \textnormal{radiation condition} & \| \mathrm{curl}\boldsymbol E \times \dfrac{x}{\|x\|} - i \omega\epsilon \boldsymbol E(x)\| &=& \mathcal O(1/\|x\|^2), & \|x\|\to \infty\,. \end{array} \)	\(\quad\quad\quad\)

lower half sphere: Dirichlet boundary with given \(\gamma_D \, E = \boldsymbol m\)
upper half sphere: Neumann boundary with given \(\gamma_N \, E = \boldsymbol j\)

Define geometry and mesh:

topsphere = Sphere((0,0,0), 1) * Box((-1,-1,0),(1,1,1))
botsphere = Sphere((0,0,0), 1) - Box((-1,-1,0),(1,1,1))
topsphere.faces.name = "neumann" # 1/kappa n x curl E
botsphere.faces.name = "dirichlet" # nxE
shape = Fuse( [topsphere,botsphere] )

h = 0.5
order = 3
mesh = Mesh( OCCGeometry(shape).GenerateMesh(maxh=h, perfstepsend=meshing.MeshingStep.MESHSURFACE)).Curve(order)
#Draw (mesh);

Define a manufactured solution \(\boldsymbol E\):

# manufactured solution: kernel with source in point s and direction e
kappa = 1. / h
s = CF((0., 0., 0.))
e = (1.192, -0.189, 2.745)
e = CF(e) / sqrt(e[0] ** 2 + e[1] ** 2 + e[2] ** 2)
xms = CF((x, y , z)) - s
r = Norm(xms)
f = (-kappa ** 2 * r ** 2 - 1j * kappa * r + 1.) / r ** 3
g = (e * xms) * (-kappa ** 2 * r ** 2 - 3. * 1j * kappa * r + 3.) / r ** 5
E = exp(1j * kappa * r) * (f * e - g * xms)
curlE = CF((E[2].Diff(y) - E[1].Diff(z), E[0].Diff(z) - E[2].Diff(x), E[1].Diff(x) - E[0].Diff(y)))

print("kappa = ", kappa)

kappa =  2.0

Define the finite element spaces, i.e., \(H^{-\frac12}(\mathrm{curl}_\Gamma, \Gamma)\) and \(H^{-\frac12}(\mathrm{div}_\Gamma, \Gamma)\) conforming spaces:

order = 3
fesHDiv = HDivSurface(mesh, order=order, complex=True)
uHDiv,vHDiv = fesHDiv.TnT() 
fesHCurl = HCurl(mesh, order=order, complex=True)
uHCurl,vHCurl = fesHCurl.TnT() 

Get the Dirichlet trace \(\boldsymbol m\) and exact Neumann data \(\boldsymbol j\) as grid functions defined everywhere on \(\Gamma\):

n = specialcf.normal(3)
m = GridFunction(fesHCurl) 
m.Set( Cross( Cross(n, E), n) , definedon=mesh.Boundaries(".*"), dual=True) # Hcurl, gamma_R = (nxE)xn
j = GridFunction(fesHDiv) 
j.Set( 1/kappa*Cross(n,curlE) , definedon=mesh.Boundaries(".*"), dual=True) # Hdiv, gamma_N = 1/k nx curl E

Define Dirichlet boundary condition on the lower half sphere \(\Gamma_0\) and have a look at it:

mDir = GridFunction(fesHCurl)
mDir.Set( m, definedon=mesh.Boundaries("dirichlet"), dual=True) # given Dirichlet data
Draw (Norm(mDir), mesh, draw_vol=False, order=3);

Define Neumann boundary condition on the upper half sphere \(\Gamma_1\) and have a look at it:

jNeu = GridFunction(fesHDiv)
jNeu.Set( j, definedon=mesh.Boundaries("neumann"), dual=True) # given Neumann data
Draw (Norm(jNeu), mesh, draw_vol=False, order=3);

Compute boundary integral operators \(\mathrm{V}\), \(\mathrm{K}\), \(\mathrm{D}\) and the mass matrix \(\mathrm M\)

# M, V, K and D 
intorder = 2 * order + 6
with TaskManager(): 
    M = BilinearForm( uHCurl.Trace() * vHDiv.Trace()* ds(bonus_intorder=3)).Assemble() # <Hcurl, Hdiv>  
    V = MaxwellSingleLayerPotentialOperator(fesHDiv, kappa, intorder=intorder, eps=1e-4)
    K = MaxwellDoubleLayerPotentialOperator(fesHCurl, fesHDiv, kappa, intorder=intorder, eps=1e-4)
    D = MaxwellSingleLayerPotentialOperatorCurl(fesHCurl, kappa, intorder=intorder, eps=1e-4)

Generate right hand side, block matrix and solve the linear system of equations:

eps = 1e-6
intorder = 2 * order + 6
with TaskManager():
    fd = ((0.5 * M.mat - K.mat) * mDir.vec - V.mat * jNeu.vec).Evaluate()
    fn = ((0.5 * M.mat.T + K.mat.T) * jNeu.vec + D.mat * mDir.vec).Evaluate()

fesHDivN = HDivSurface(mesh, order=order, dirichlet="neumann") # Dirichlet nodes free dofs
fesHCurlD = HCurl(mesh, order=order, dirichlet="dirichlet") # Neumann nodes free dofs

lhs = BlockMatrix([[V.mat, K.mat], [-K.mat.T, -D.mat]]) 
rhs = BlockVector([fd, fn])
preHDiv = BilinearForm( uHDiv.Trace() * vHDiv.Trace() *ds).Assemble()
preHCurl = BilinearForm( uHCurl.Trace() * vHCurl.Trace() *ds).Assemble()
preBlock = BlockMatrix ([[preHDiv.mat.Inverse(freedofs=fesHDivN.FreeDofs()), None], [None, preHCurl.mat.Inverse(freedofs=fesHCurlD.FreeDofs())]])

sol = GMRes(A=lhs, b=rhs, pre=preBlock, maxsteps=3000, tol=1e-8, printrates=False)

Check the numerical solution

Have a look at Neumann data on \(\Gamma_0\) and compute the \(L_2\)-error:

# compare with exact Dirichlet trace
gfdiv = GridFunction(fesHDiv)
gfdiv.vec[:] = sol[0]
print ("Neumann error =", sqrt(Integrate(Norm(gfdiv + jNeu - j)**2, mesh.Boundaries(".*"), BND)))
Draw (Norm(gfdiv), mesh.Boundaries("dirichlet"), order=3);

Neumann error = 0.021621506900978364

Have a look at Dirichlet data on \(\Gamma_1\) and compute the \(L_2\)-error:

# compare with the exact Neuman trace  
gfcurl = GridFunction(fesHCurl)
gfcurl.vec[:] = sol[1]
print ("Dirichlet error =", sqrt(Integrate(Norm(gfcurl +mDir - m)**2, mesh.Boundaries(".*"), BND)))
Draw (Norm(gfcurl), mesh.Boundaries("neumann"), draw_vol=False, order=3);

Dirichlet error = 0.029447011579737066

References (theory results):

Boundary Element Methods for Maxwell Transmission Problems in Lipschitz Domains

Note: the Calderon-projector is in the referenced paper is slightly different from ours as only trace space \(H^{-\frac12}(\mathrm{div}_\Gamma,\Gamma)\) is used). Theoretical results stay the same.