Maxwell solver for PEC bodies (EFIE)

keys: Maxwell single layer potential, EFIE, indirect ansatz, MoM, PEC scattering

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

Loading ngbem library

We consider a perfect conducter \(\Omega \subset \mathbb R^3\) and a plane wave \(\boldsymbol E_{\mathrm{inc}} \) with tangential trace \(\gamma_R \boldsymbol E_{\mathrm{inc}} = -\boldsymbol m \in \boldsymbol H^{-\frac12}\left( \mathrm{curl}_\Gamma, \Gamma\right)\,.\) The incoming wave thus induces a scattered electric field \(\boldsymbol E\) which propagates into \(\Omega^c\). The scattered electric field solves the following boundary value problem:

\[\begin{split} \left\{ \begin{array}{rcl l} \mathbf{curl} \, \mathbf{curl}\, \boldsymbol E - \kappa^2 \, \boldsymbol E &=& \boldsymbol 0, \quad &\textnormal{in } \Omega^c \subset \mathbb R^3\,,\\ \gamma_R \,\boldsymbol E &=& \boldsymbol m, \quad & \textnormal{on }\Gamma \\ \left\| \mathbf{curl} \, \boldsymbol E( x) - i\,\omega\,\epsilon \, \boldsymbol E( x)\right\| &=& \mathcal O\left( \displaystyle \frac{1}{\| x\|^2}\right), &\textnormal{for} \; \|x\| \to \infty\,.\end{array} \right. \end{split}\]

and a possible representation reads

\[ \boldsymbol E(x) = \underbrace{ \kappa \, \int_\Gamma \displaystyle{ \frac{1}{4\,\pi} \, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \boldsymbol j(y)\, \mathrm{d}\sigma_y + \frac{1}{\kappa} \nabla \int_\Gamma \displaystyle{ \frac{1}{4\,\pi}\, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \mathrm{div}_\Gamma \boldsymbol j(y)\, \mathrm{d}\sigma_y }_{\displaystyle{ \mathrm{SL} (\boldsymbol j) } } \,.\]

The density \(\boldsymbol j\) is determined by the boundary element method, i.e. the numerical solution of the variational formulation

\[ \begin{array}{c} \forall \, \boldsymbol v\in H^{-\frac12}(\mathrm{div}_\Gamma, \Gamma): \quad \left\langle \,\mathrm{SL} (\boldsymbol j),\, \boldsymbol v \right\rangle_{-\frac12} = \left\langle \boldsymbol m , \boldsymbol v\right\rangle_{-\frac12} \,. \end{array}\]

In the enineering community, the approximation scheme is also known as method of moments (MoM) and the resulting equation for \(\boldsymbol j\) is called the EFIE, the electric field integral equation.

Define the geometry of the perfect conductor \(\Omega\) and create a mesh:

sp = Sphere( (0,0,0), 1)
mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=1, perfstepsend=meshing.MeshingStep.MESHSURFACE)).Curve(4)

Next, we create finite element spaces \(\boldsymbol H^{-\frac12}(\mathrm{curl}_\Gamma, \Gamma)\) and \(\boldsymbol H^{-\frac12}(\mathrm{div}_\Gamma, \Gamma)\):

fesHCurl = HCurl(mesh, order=3, complex=True)
fesHDiv = HDivSurface(mesh, order=3, complex=True)
uHCurl,vHCurl = fesHCurl.TnT() # H(curl_Gamma) trial for Dirichlet data ( (nxE)xn )
uHDiv,vHDiv = fesHDiv.TnT() # H(div_Gamma) trial space for Neumann data ( nx curlE ) and test space for BIE

print ("ndof HCurl = ", fesHCurl.ndof)
print ("ndof HDiv = ", fesHDiv.ndof)

ndof HCurl =  1568
ndof HDiv =  1568

Define the incoming plane wave and compute the given Dirichlet data \(\boldsymbol m\):

eps0 = 8.854e-12 
mu0 = 4*pi*1e-7
omega = 1.5e9
kappa = omega*sqrt(eps0*mu0)
print("kappa = ", kappa)

E_inc = CF((1,0,0))*exp( -1j * kappa * z )
n = specialcf.normal(3)
m = - Cross( Cross(n, E_inc), n) # m = (nxE)xn in H(curl_Gamma)
Draw(Norm(m), mesh, draw_vol=False, order=2) ;

kappa =  5.0034083602476045

Compute the right hand side vector:

rhs = LinearForm( m * vHDiv.Trace() *ds(bonus_intorder=3) ).Assemble()  # <H(curl_Gamma), H(div_Gamma)> 

The discretisation of the above variational formulation leads to a system of linear equations, ie

\[ V\, \mathrm j = \mathrm{rhs}\,,\]

where the single layer operator \(V\) on the left is a regular and symmetric matrix.

j = GridFunction(fesHDiv)
pre = BilinearForm( uHDiv.Trace() * vHDiv.Trace() *ds ).Assemble().mat.Inverse(freedofs=fesHDiv.FreeDofs()) 
with TaskManager(): 
    V = MaxwellSingleLayerPotentialOperator(fesHDiv, kappa, 
                                            intorder=10, leafsize=40, eta=3., eps=1e-4, method="aca", testhmatrix=False)
    CG(mat = V.mat, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, maxsteps=500, initialize=False, printrates=False)

Draw (Norm(j), mesh, draw_vol=False, order=3);

Draw( j, mesh, draw_vol=False, order=3, min=-2, max=2, animate_complex=True);

postprocessing: evaluate the scattered field

Define the geometry of the receiver, generate a mesh and evaluate the scattered field \(\boldsymbol E\) in each node:

# sphere 
screen = Sphere( (0,0,0), 20)
mesh_screen = Mesh( OCCGeometry(screen).GenerateMesh(maxh=1, perfstepsend=meshing.MeshingStep.MESHSURFACE)).Curve(4)
fes_screen = HCurl(mesh_screen, order=3, complex=True)
E_screen = GridFunction(fes_screen)
E_screen.Set (V.GetPotential(j)[0:3], definedon=mesh_screen.Boundaries(".*"))

We obtain the typical far field pattern of the scattered field:

Draw(kappa*300*Norm(E_screen)*n, mesh_screen, deformation=True);

# plane screen 
screen = WorkPlane(Axes( (0,0,-3.5), Z, X)).RectangleC(20,20).Face()
mesh_screen = Mesh( OCCGeometry(screen).GenerateMesh(maxh=1)).Curve(1)
fes_screen = HCurl(mesh_screen, order=3, complex=True)
E_screen = GridFunction(fes_screen)
E_screen.Set (V.GetPotential(j)[0:3], definedon=mesh_screen.Boundaries(".*"))

Draw( E_screen, mesh_screen, draw_vol=False, order=3, min=0, max=0.1, animate_complex=True);

For details on convergence rates and low frequency stabilisation of the EFIE look here.