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

keys: PEC scattering, single layer potential, EFIE, MoM

Maxwell PEC Indirect Method

We consider a perfect conductor \(\Omega \subset \mathbb R^3\) emitting an electric field into \(\Omega^c\). The scattered electric field \(\boldsymbol E\) solves the following Dirichlet 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}\]

The electric field \(\boldsymbol E\) is given by

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

Define the geometry of the perfect conductor \(\Omega\), create a mesh and test and trial functions according to the mesh:

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

fesHCurl = HCurl(mesh, order=3, complex=True)
uHCurl,vHCurl = fesHCurl.TnT() # H(curl_Gamma) trial for Dirichlet data ( (nxE)xn )

fesHDiv = HDivSurface(mesh, order=3, complex=True)
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

Boundary Integral Equation

We carefully apply the tangential trace to (1) and obtain the electric field integral equation (EFIE). The EFIE is a boundary integral equation for \(\boldsymbol j\):

\[ \forall \, \boldsymbol v\in H^{-\frac12}(\mathrm{div}_\Gamma, \Gamma): \quad \displaystyle \int\limits_\Gamma \mathrm{SL}(\boldsymbol j) \cdot \boldsymbol v \, \mathrm d s = \displaystyle \int\limits_\Gamma \boldsymbol m \cdot \boldsymbol v \, \mathrm d s\,. \]

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

\[ \mathrm V\, \mathrm j = \mathrm{M} \, \mathrm m\,, \]

where

• \(\mathrm V\) is the Maxwell single layer potential operator.

• \(\mathrm M\) is a mass matrix.

The Maxwell single layer operator \(\mathrm V\) is a derivate of the Helmholtz single layer operator (details given below) and thus, in NGSBEM we use the Helmholtz operators as building blocks for Maxwell operator:

V = kappa * HelmholtzSL(u*ds, kappa)*v*ds
    - 1/kappa * HelmholtzSL(div(u)*ds, kappa)*div(v)*ds

rhs = LinearForm( m * vHDiv.Trace()
                 *ds(bonus_intorder=3) ).Assemble()

with TaskManager(): 
    V1 = HelmholtzSL( 
            uHDiv.Trace()*ds(bonus_intorder=6), 
            kappa
        )*vHDiv.Trace()*ds(bonus_intorder=6)
    V2 = HelmholtzSL( 
            div(uHDiv.Trace())*ds(bonus_intorder=6), 
            kappa
        )*div(vHDiv.Trace())*ds(bonus_intorder=6)
    V = kappa * V1.mat - 1/kappa * V2.mat

pre = BilinearForm( 
        uHDiv.Trace() * vHDiv.Trace()*ds 
                  ).Assemble().mat.Inverse(freedofs=fesHDiv.FreeDofs()) 

j = GridFunction(fesHDiv)
with TaskManager():
    CG(mat = V, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, 
       maxsteps=500, initialize=False, printrates=True)

CG iteration 1, residual = 0.8878706955909835     

CG iteration 2, residual = 0.7675484948921711     

CG iteration 3, residual = 0.6522908847399078     

CG iteration 4, residual = 0.5438772154771219     

CG iteration 5, residual = 0.32577579742123597     

CG iteration 6, residual = 0.38285182057836126     

CG iteration 7, residual = 0.46451571786011586     

CG iteration 8, residual = 0.5016221047913395     

CG iteration 9, residual = 0.4829068615796905     

CG iteration 10, residual = 0.5087289590428536     

CG iteration 11, residual = 0.3124684395022449     

CG iteration 12, residual = 0.29831571689106484     

CG iteration 13, residual = 1.1125732551773957     

CG iteration 14, residual = 0.3039874045721664     

CG iteration 15, residual = 0.09063205902539472     

CG iteration 16, residual = 0.07227404776232635     

CG iteration 17, residual = 0.37870365032344594     

CG iteration 18, residual = 0.061296805106396915     

CG iteration 19, residual = 0.12918056268698355     

CG iteration 20, residual = 0.2184497915360951     

CG iteration 21, residual = 0.08357025497105079     

CG iteration 22, residual = 0.09975574237240499     

CG iteration 23, residual = 0.0899104312491359     

CG iteration 24, residual = 0.2478338879601313     

CG iteration 25, residual = 0.07719944563982481     

CG iteration 26, residual = 0.0746870863731881     

CG iteration 27, residual = 0.04013173979562398     

CG iteration 28, residual = 0.012430995371250532     

CG iteration 29, residual = 0.011834162914665804     

CG iteration 30, residual = 0.015615968411721918     

CG iteration 31, residual = 0.08045099702406802     

CG iteration 32, residual = 0.1670677184298732     

CG iteration 33, residual = 0.27422610913927864     

CG iteration 34, residual = 0.1748296060265163     

CG iteration 35, residual = 0.16658275177320758     

CG iteration 36, residual = 0.1158642875612813     

CG iteration 37, residual = 0.14497092246790735     

CG iteration 38, residual = 0.3733772697315533     

CG iteration 39, residual = 0.014658401141416375     

CG iteration 40, residual = 0.009788582578760608     

CG iteration 41, residual = 0.02069206181196535     

CG iteration 42, residual = 0.012942468692378327     

CG iteration 43, residual = 0.014620752826192726     

CG iteration 44, residual = 0.01897407098137275     

CG iteration 45, residual = 0.03433029321966666     

CG iteration 46, residual = 0.127581842921377     

CG iteration 47, residual = 0.008621139345777542     

CG iteration 48, residual = 0.009198636036351798     

CG iteration 49, residual = 0.01049940016712042     

CG iteration 50, residual = 0.08355288826292584     

CG iteration 51, residual = 0.007881785961896983     

CG iteration 52, residual = 0.02253646991625219     

CG iteration 53, residual = 0.0018846345931802636     

CG iteration 54, residual = 0.0017049091642830662     

CG iteration 55, residual = 0.0029496748695508365     

CG iteration 56, residual = 0.0020120874123977857     

CG iteration 57, residual = 0.003952284126905005     

CG iteration 58, residual = 0.005271139161398289     

CG iteration 59, residual = 0.0012182973376681822     

CG iteration 60, residual = 0.0011498854352651208     

CG iteration 61, residual = 0.001384358655438875     

CG iteration 62, residual = 0.001614220555336213     

CG iteration 63, residual = 0.0054668972721325735     

CG iteration 64, residual = 0.004179666615966757     

CG iteration 65, residual = 0.002414549270890931     

CG iteration 66, residual = 0.0015495867311260958     

CG iteration 67, residual = 0.0005506876201262699     

CG iteration 68, residual = 0.00037660293082354175     

CG iteration 69, residual = 0.0034338782327507133     

CG iteration 70, residual = 0.0009471436537970873     

CG iteration 71, residual = 0.0014335818876241675     

CG iteration 72, residual = 0.0007283371372882074     

CG iteration 73, residual = 0.0005571482212848026     

CG iteration 74, residual = 0.0012157298897071244     

CG iteration 75, residual = 0.0007906937337094413     

CG iteration 76, residual = 0.0008794257320741882     

CG iteration 77, residual = 0.0003573931097134489     

CG iteration 78, residual = 0.00036376154502061993     

CG iteration 79, residual = 0.00011300801443305539     

CG iteration 80, residual = 0.00011047679215769697     

CG iteration 81, residual = 0.0002621509941934966     

CG iteration 82, residual = 0.0005983862184124033     

CG iteration 83, residual = 0.017819006324953864     

CG iteration 84, residual = 0.00028640420629043035     

CG iteration 85, residual = 0.00013746058209406902     

CG iteration 86, residual = 0.00011563976332848164     

CG iteration 87, residual = 0.00017042544893911964     

CG iteration 88, residual = 0.00025077530903512134     

CG iteration 89, residual = 0.00022189693818589547     

CG iteration 90, residual = 0.00020923611729538357     

CG iteration 91, residual = 0.0001235104815938112     

CG iteration 92, residual = 0.0001566595966006561     

CG iteration 93, residual = 0.00011844280859338196     

CG iteration 94, residual = 0.00017389182990587709     

CG iteration 95, residual = 0.00016887106547974323     

CG iteration 96, residual = 0.00018952996702171237     

CG iteration 97, residual = 0.00015974400684772805     

CG iteration 98, residual = 0.00014035005288671497     

CG iteration 99, residual = 0.0004398378078936469     

CG iteration 100, residual = 0.0001143854414982216     

CG iteration 101, residual = 0.0001256947253932601     

CG iteration 102, residual = 0.00014133330428478638     

CG iteration 103, residual = 4.237199529532093e-05     

CG iteration 104, residual = 4.1467757407425625e-05     

CG iteration 105, residual = 4.051746901657601e-05     

CG iteration 106, residual = 7.819385067735611e-05     

CG iteration 107, residual = 3.741501755319228e-05     

CG iteration 108, residual = 2.1940649927933642e-05     

CG iteration 109, residual = 3.140249468387033e-05     

CG iteration 110, residual = 4.368064954972929e-05     

CG iteration 111, residual = 0.00014788397151472202     

CG iteration 112, residual = 3.2932701038636595e-05     

CG iteration 113, residual = 0.0001617408482121544     

CG iteration 114, residual = 1.5008397391182214e-05     

CG iteration 115, residual = 1.1179809062833575e-05     

CG iteration 116, residual = 1.4795379855386945e-05     

CG iteration 117, residual = 1.7240022921987243e-05     

CG iteration 118, residual = 2.0853703674160856e-05     

CG iteration 119, residual = 4.2984169972410184e-05     

CG iteration 120, residual = 2.228812927053112e-05     

CG iteration 121, residual = 2.7306021354902898e-05     

CG iteration 122, residual = 3.877343126777903e-05     

CG iteration 123, residual = 4.8128536569088115e-05     

CG iteration 124, residual = 1.4811088400367004e-05     

CG iteration 125, residual = 6.913400109768796e-06     

CG iteration 126, residual = 4.705607554562901e-06     

CG iteration 127, residual = 4.640261600686317e-06     

CG iteration 128, residual = 9.111615957124439e-06     

CG iteration 129, residual = 1.1812665177556363e-05     

CG iteration 130, residual = 1.0094893849176135e-05     

CG iteration 131, residual = 6.44467411675827e-06     

CG iteration 132, residual = 2.746949334142758e-06     

CG iteration 133, residual = 3.6279682777650116e-06     

CG iteration 134, residual = 6.349778574774744e-06     

CG iteration 135, residual = 9.717864575247483e-06     

CG iteration 136, residual = 2.070139390402094e-06     

CG iteration 137, residual = 3.034341892590199e-06     

CG iteration 138, residual = 1.8487287102176112e-06     

CG iteration 139, residual = 3.4580320505720104e-06     

CG iteration 140, residual = 1.4177322797807481e-05     

CG iteration 141, residual = 3.8779036823823085e-06     

CG iteration 142, residual = 7.394657235204936e-06     

CG iteration 143, residual = 3.893264467080583e-06     

CG iteration 144, residual = 3.746213378117184e-06     

CG iteration 145, residual = 2.799397181502366e-06     

CG iteration 146, residual = 4.513990755235639e-06     

CG iteration 147, residual = 6.0916878470123024e-06     

CG iteration 148, residual = 4.931980695957252e-06     

CG iteration 149, residual = 1.0007061921678747e-05     

CG iteration 150, residual = 9.251734061600442e-06     

CG iteration 151, residual = 1.7367445299434495e-06     

CG iteration 152, residual = 1.601620001946759e-05     

CG iteration 153, residual = 1.3430003313269798e-06     

CG iteration 154, residual = 6.09375017244625e-06     

CG iteration 155, residual = 1.3334705911915656e-06     

CG iteration 156, residual = 3.5330039506940004e-06     

CG iteration 157, residual = 2.3016504146103222e-06     

CG iteration 158, residual = 8.88835495793434e-07     

CG iteration 159, residual = 1.3229761000661084e-06     

CG iteration 160, residual = 3.5956745810159612e-06     

CG iteration 161, residual = 8.057589106262411e-06     

CG iteration 162, residual = 7.749772168122223e-06     

CG iteration 163, residual = 5.512733990976164e-06     

CG iteration 164, residual = 2.017822549460402e-06     

CG iteration 165, residual = 1.8155623877012802e-06     

CG iteration 166, residual = 1.2547874209443942e-06     

CG iteration 167, residual = 3.7291165371356567e-06     

CG iteration 168, residual = 6.575521966783579e-06     

CG iteration 169, residual = 5.296062209493131e-06     

CG iteration 170, residual = 3.094999485786576e-06     

CG iteration 171, residual = 4.177698832552698e-06     

CG iteration 172, residual = 1.8578858663060847e-06     

CG iteration 173, residual = 1.5710926923333698e-06     

CG iteration 174, residual = 1.6309346495504602e-06     

CG iteration 175, residual = 2.599638475279503e-06     

CG iteration 176, residual = 6.3179774822161685e-06     

CG iteration 177, residual = 3.0607465135132463e-06     

CG iteration 178, residual = 3.996107793901329e-06     

CG iteration 179, residual = 1.3416502435095672e-06     

CG iteration 180, residual = 3.5717036058137423e-06     

CG iteration 181, residual = 8.949629394973365e-07     

CG iteration 182, residual = 2.871642190441616e-06     

CG iteration 183, residual = 4.822820456504448e-07     

CG iteration 184, residual = 3.799975020065045e-07     

CG iteration 185, residual = 4.984054587733439e-07     

CG iteration 186, residual = 4.599843172716155e-07     

CG iteration 187, residual = 6.584580173181138e-07     

CG iteration 188, residual = 1.5286457061810682e-06     

CG iteration 189, residual = 8.88828355569564e-07     

CG iteration 190, residual = 8.355506838006993e-07     

CG iteration 191, residual = 6.988165084487581e-07     

CG iteration 192, residual = 1.2835197230970112e-07     

CG iteration 193, residual = 1.2144887752276157e-07     

CG iteration 194, residual = 2.7192428405666655e-07     

CG iteration 195, residual = 1.1701564248129498e-06     

CG iteration 196, residual = 1.2943262031274369e-06     

CG iteration 197, residual = 1.093194601093173e-06     

CG iteration 198, residual = 1.2702860869467277e-06     

CG iteration 199, residual = 3.429579847737158e-07     

CG iteration 200, residual = 6.696254894359548e-07     

CG iteration 201, residual = 2.9360063385869755e-07     

CG iteration 202, residual = 7.067743438723994e-07     

CG iteration 203, residual = 1.1533909082965132e-07     

CG iteration 204, residual = 3.051098203810484e-07     

CG iteration 205, residual = 9.627569857119783e-08     

CG iteration 206, residual = 9.372965487799278e-08     

CG iteration 207, residual = 1.2497383959743207e-07     

CG iteration 208, residual = 1.7686286756547529e-07     

CG iteration 209, residual = 6.094742988840986e-08     

CG iteration 210, residual = 9.925209950113864e-08     

CG iteration 211, residual = 1.7642786334224647e-07     

CG iteration 212, residual = 8.051446194056417e-08     

CG iteration 213, residual = 5.0157366786485836e-08     

CG iteration 214, residual = 1.1877201728167736e-07     

CG iteration 215, residual = 1.0353285487827848e-07     

CG iteration 216, residual = 4.4415896103953655e-07     

CG iteration 217, residual = 7.507194201808819e-08     

CG iteration 218, residual = 4.470647733978831e-08     

CG iteration 219, residual = 4.465788871105328e-08     

CG iteration 220, residual = 5.6756713576630624e-08     

CG iteration 221, residual = 8.945833638541662e-08     

CG iteration 222, residual = 7.23101464067061e-08     

CG iteration 223, residual = 9.512742781510455e-08     

CG iteration 224, residual = 2.3822268397277802e-07     

CG iteration 225, residual = 2.7459436244400853e-08     

CG iteration 226, residual = 5.7005074368670836e-08     

CG iteration 227, residual = 3.519943737143964e-08     

CG iteration 228, residual = 2.1699135422678475e-08     

CG iteration 229, residual = 3.36538373805336e-08     

CG iteration 230, residual = 2.1215503435753443e-08     

CG iteration 231, residual = 4.898983092748078e-08     

CG iteration 232, residual = 3.584118835146107e-08     

CG iteration 233, residual = 2.642485858362718e-08     

CG iteration 234, residual = 4.83277771535435e-08     

CG iteration 235, residual = 2.916678719407789e-08     

CG iteration 236, residual = 2.9076237839396328e-08     

CG iteration 237, residual = 2.5201955348632453e-08     

CG iteration 238, residual = 2.7365862116911207e-08     

CG iteration 239, residual = 1.4708896913680368e-08     

CG iteration 240, residual = 2.10930551266155e-08     

CG iteration 241, residual = 1.5941959844731054e-08     

CG iteration 242, residual = 2.4482146590114438e-08     

CG iteration 243, residual = 1.7445183462121708e-08     

CG iteration 244, residual = 1.5984276632655868e-08     

CG iteration 245, residual = 3.744927928830659e-09     

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

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

Warning in Webgui.Draw: function has more than 3 components, drawing Norm(cf) instead

Evaluation of the Solution

Define the geometry of the receiver, generate a mesh and evaluate the scattered field \(\boldsymbol E\) as given by the representation formula, i.e.,

\[ \textnormal{ for } \quad x \quad \textnormal{ on receiver, it holds } \quad \boldsymbol E(x) = \underbrace{ \kappa \, \int\limits_\Gamma \displaystyle{ \frac{1}{4\,\pi} \, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \boldsymbol j(y)\, \mathrm{d}s_y + \frac{1}{\kappa} \, \nabla \int\limits_\Gamma \displaystyle{ \frac{1}{4\,\pi}\, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \mathrm{div}_\Gamma \boldsymbol j(y)\, \mathrm{d}s_y}_{\displaystyle{ \mathrm{SL}(\boldsymbol j)}} \,.\]

• consider as receiver a sphere around the PEC body:

# define receiver, i.e., a sphere around pec
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)

# evaluate represenation formula
repformula = kappa * HelmholtzSL(uHDiv.Trace()*ds, kappa)(j) + 1/ kappa * grad(HelmholtzSL(div(uHDiv.Trace())*ds, kappa)(j))
E_screen.Set (repformula, definedon=mesh_screen.Boundaries(".*"))
Draw( E_screen, mesh_screen, draw_vol=False, order=3, min=0, max=0.1, animate_complex=True);

Warning in Webgui.Draw: function has more than 3 components, drawing Norm(cf) instead

# far field pattern
Draw(kappa*300*Norm(E_screen)*n, mesh_screen, deformation=True);

• consider as receiver a plane screen in the back of the PEC body:

# define the receiver, i.e., plane screen (note that the potential evaluation does not yet evaluate the full SL)
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)

# evaluate represenation formula
repformula = kappa * HelmholtzSL(uHDiv.Trace()*ds, kappa)(j) + 1/ kappa * grad(HelmholtzSL(div(uHDiv.Trace())*ds, kappa)(j))
E_screen.Set (repformula, definedon=mesh_screen.Boundaries(".*"))

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

Warning in Webgui.Draw: function has more than 3 components, drawing Norm(cf) instead

Details: Explicit Representation of the Maxwell Single Layer Potential

For a trial function \(\boldsymbol u_j \in H^{-\frac12}(\mathrm{div}_\Gamma,\Gamma)\) and a test function \(\boldsymbol v_i \in H^{-\frac12}(\mathrm{div}_\Gamma,\Gamma)\) the Maxwell single layer potential operator entry \(V_{ij}\) can formally written as

\[\begin{split} \begin{array}{rcl} V_{ij} &=& \kappa \, \displaystyle \int\limits_\Gamma \underbrace{ \displaystyle \int\limits_\Gamma \displaystyle{ \frac{1}{4\,\pi} \, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \boldsymbol u_j(y) \cdot \boldsymbol v_i(x) \, \mathrm{d}s_x}_{\displaystyle{\mathrm{HelmholtzSL}(\boldsymbol v_i) } } \mathrm{d} s_x \\ && \quad - \; \frac{1}{\kappa} \displaystyle \int\limits_\Gamma \underbrace{ \displaystyle \int\limits_\Gamma \displaystyle{ \frac{1}{4\,\pi}\, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \mathrm{div}_\Gamma \boldsymbol u_j(y) }_{\displaystyle{ \mathrm{HelmholtzSL}(\mathrm{div}_\Gamma(\boldsymbol u_j) } } \, \mathrm{div}_\Gamma \boldsymbol v_i(x)\, \mathrm{d}\sigma_y\,\mathrm d s_x\,. \end{array} \end{split}\]

References

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