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.7675484949610056     

CG iteration 3, residual = 0.6522908847428317     

CG iteration 4, residual = 0.5438772154441428     

CG iteration 5, residual = 0.3257757977899446     

CG iteration 6, residual = 0.38285181787495864     

CG iteration 7, residual = 0.46451579607613397     

CG iteration 8, residual = 0.5016220896455411     

CG iteration 9, residual = 0.4829068660391469     

CG iteration 10, residual = 0.5087290260771855     

CG iteration 11, residual = 0.3124684494631325     

CG iteration 12, residual = 0.29831572677186863     

CG iteration 13, residual = 1.112573077357725     

CG iteration 14, residual = 0.3039873741663487     

CG iteration 15, residual = 0.09063208237073635     

CG iteration 16, residual = 0.0722740585129648     

CG iteration 17, residual = 0.37870343757146147     

CG iteration 18, residual = 0.06129680638073469     

CG iteration 19, residual = 0.12918060497854805     

CG iteration 20, residual = 0.2184497431617648     

CG iteration 21, residual = 0.08357019799517788     

CG iteration 22, residual = 0.09975556095539428     

CG iteration 23, residual = 0.08991033414317959     

CG iteration 24, residual = 0.24783452857692412     

CG iteration 25, residual = 0.07719947908764722     

CG iteration 26, residual = 0.07468648266863184     

CG iteration 27, residual = 0.040131409835212184     

CG iteration 28, residual = 0.012430972073777946     

CG iteration 29, residual = 0.011834129178266832     

CG iteration 30, residual = 0.015615949007801971     

CG iteration 31, residual = 0.08045098578908988     

CG iteration 32, residual = 0.1670677335189299     

CG iteration 33, residual = 0.27422615540708684     

CG iteration 34, residual = 0.1748296218171257     

CG iteration 35, residual = 0.16658272795787657     

CG iteration 36, residual = 0.11586427704119147     

CG iteration 37, residual = 0.14497102402910936     

CG iteration 38, residual = 0.3733764949751639     

CG iteration 39, residual = 0.014658403759615038     

CG iteration 40, residual = 0.009788555551363057     

CG iteration 41, residual = 0.020691944466949573     

CG iteration 42, residual = 0.01294247249019566     

CG iteration 43, residual = 0.014620745965870445     

CG iteration 44, residual = 0.018973872840476847     

CG iteration 45, residual = 0.03432920684552786     

CG iteration 46, residual = 0.12757098647571896     

CG iteration 47, residual = 0.008621347990286668     

CG iteration 48, residual = 0.009199023114811517     

CG iteration 49, residual = 0.010499993674001172     

CG iteration 50, residual = 0.0835379400763569     

CG iteration 51, residual = 0.007881965737981354     

CG iteration 52, residual = 0.02253659791288671     

CG iteration 53, residual = 0.001884664198260332     

CG iteration 54, residual = 0.001704931670081129     

CG iteration 55, residual = 0.0029496869953803215     

CG iteration 56, residual = 0.00201209385283071     

CG iteration 57, residual = 0.003952297534871801     

CG iteration 58, residual = 0.005271099941983281     

CG iteration 59, residual = 0.0012179945346108944     

CG iteration 60, residual = 0.0011483451162070573     

CG iteration 61, residual = 0.0013840449848485734     

CG iteration 62, residual = 0.0016143663506699528     

CG iteration 63, residual = 0.00546683466673916     

CG iteration 64, residual = 0.004181168647779056     

CG iteration 65, residual = 0.0024136557639763867     

CG iteration 66, residual = 0.0015276847100682894     

CG iteration 67, residual = 0.00048792401888450137     

CG iteration 68, residual = 0.0002957345787988464     

CG iteration 69, residual = 0.00045486447806173134     

CG iteration 70, residual = 0.0008464103180028784     

CG iteration 71, residual = 0.0027330142353865212     

CG iteration 72, residual = 0.0007392763604560441     

CG iteration 73, residual = 0.0005545083077893293     

CG iteration 74, residual = 0.0012032577594203472     

CG iteration 75, residual = 0.0008032621041076767     

CG iteration 76, residual = 0.0008842075482257568     

CG iteration 77, residual = 0.000377653140396701     

CG iteration 78, residual = 0.00018913550238964212     

CG iteration 79, residual = 0.00019389924305359595     

CG iteration 80, residual = 0.00013835353050295178     

CG iteration 81, residual = 0.0002983331823334408     

CG iteration 82, residual = 0.00121703579594617     

CG iteration 83, residual = 0.0026218977064665135     

CG iteration 84, residual = 0.00015631118782924813     

CG iteration 85, residual = 0.00029969910216110006     

CG iteration 86, residual = 9.844045781061686e-05     

CG iteration 87, residual = 0.0001708251924757036     

CG iteration 88, residual = 0.00020620114770152466     

CG iteration 89, residual = 0.00024048110638942862     

CG iteration 90, residual = 0.00014721122020132295     

CG iteration 91, residual = 0.0006563940696334723     

CG iteration 92, residual = 0.00048447221501472196     

CG iteration 93, residual = 0.0001498666834451406     

CG iteration 94, residual = 0.00029550729438994187     

CG iteration 95, residual = 0.00017815042157461042     

CG iteration 96, residual = 9.753640028906295e-05     

CG iteration 97, residual = 9.28823592708103e-05     

CG iteration 98, residual = 0.0001593829510918941     

CG iteration 99, residual = 0.00014869735609933195     

CG iteration 100, residual = 8.467715654545165e-05     

CG iteration 101, residual = 0.00012054193135396488     

CG iteration 102, residual = 0.00013416530014904933     

CG iteration 103, residual = 4.259099925681084e-05     

CG iteration 104, residual = 4.1408086205220534e-05     

CG iteration 105, residual = 4.0705042270415105e-05     

CG iteration 106, residual = 7.833777846955292e-05     

CG iteration 107, residual = 3.737570581488143e-05     

CG iteration 108, residual = 2.2105495163994154e-05     

CG iteration 109, residual = 3.1149679304083046e-05     

CG iteration 110, residual = 0.0001074762579451603     

CG iteration 111, residual = 7.165184529414234e-05     

CG iteration 112, residual = 4.154545027833018e-05     

CG iteration 113, residual = 9.810466774541883e-06     

CG iteration 114, residual = 1.9397857831425178e-05     

CG iteration 115, residual = 1.4440734947229271e-05     

CG iteration 116, residual = 2.1110299143862167e-05     

CG iteration 117, residual = 2.6973955236166043e-05     

CG iteration 118, residual = 2.0983918627436867e-05     

CG iteration 119, residual = 4.875277500178075e-05     

CG iteration 120, residual = 2.8723868382070037e-05     

CG iteration 121, residual = 3.179887920081959e-05     

CG iteration 122, residual = 3.9599678762537026e-05     

CG iteration 123, residual = 3.267718409775318e-05     

CG iteration 124, residual = 1.3211756357484176e-05     

CG iteration 125, residual = 6.030038733369286e-06     

CG iteration 126, residual = 1.1863749638935548e-05     

CG iteration 127, residual = 8.381252840968402e-06     

CG iteration 128, residual = 1.1425318785013613e-05     

CG iteration 129, residual = 1.819324463722697e-05     

CG iteration 130, residual = 6.04350642466939e-06     

CG iteration 131, residual = 2.613783843324706e-06     

CG iteration 132, residual = 4.699015899489147e-06     

CG iteration 133, residual = 2.709415478967247e-06     

CG iteration 134, residual = 3.141143603589289e-06     

CG iteration 135, residual = 2.0715276156063946e-06     

CG iteration 136, residual = 4.472502493259666e-06     

CG iteration 137, residual = 1.3181885621269223e-05     

CG iteration 138, residual = 1.8026672932039892e-06     

CG iteration 139, residual = 2.089580771588785e-06     

CG iteration 140, residual = 1.7690844794489872e-06     

CG iteration 141, residual = 1.9785689650333993e-06     

CG iteration 142, residual = 3.0773866555749024e-06     

CG iteration 143, residual = 2.190863818758248e-06     

CG iteration 144, residual = 4.048018916742281e-06     

CG iteration 145, residual = 2.913061139745785e-06     

CG iteration 146, residual = 1.3801895825077007e-05     

CG iteration 147, residual = 2.857064310445932e-06     

CG iteration 148, residual = 5.15047915941731e-06     

CG iteration 149, residual = 1.8708166515395321e-06     

CG iteration 150, residual = 1.460958825912483e-06     

CG iteration 151, residual = 2.3643934237110708e-06     

CG iteration 152, residual = 1.7284871989718608e-06     

CG iteration 153, residual = 7.120790061866152e-07     

CG iteration 154, residual = 8.367591486193599e-07     

CG iteration 155, residual = 7.968494270454812e-07     

CG iteration 156, residual = 1.5791473529732477e-06     

CG iteration 157, residual = 2.4630869994091505e-06     

CG iteration 158, residual = 2.886901215967435e-06     

CG iteration 159, residual = 3.889873190732558e-06     

CG iteration 160, residual = 2.938314986679838e-06     

CG iteration 161, residual = 1.9059196596153892e-06     

CG iteration 162, residual = 1.8560158676523714e-06     

CG iteration 163, residual = 1.5437328442567165e-06     

CG iteration 164, residual = 1.6199472335415415e-06     

CG iteration 165, residual = 2.070088900050721e-06     

CG iteration 166, residual = 2.3560988807223344e-06     

CG iteration 167, residual = 2.49404830927336e-06     

CG iteration 168, residual = 3.120285515377673e-05     

CG iteration 169, residual = 3.7915337001113784e-06     

CG iteration 170, residual = 4.351682164492531e-06     

CG iteration 171, residual = 2.8412590177363748e-05     

CG iteration 172, residual = 3.6530475010630547e-06     

CG iteration 173, residual = 3.011156360798155e-06     

CG iteration 174, residual = 2.5666330915351435e-06     

CG iteration 175, residual = 3.0798695892552687e-06     

CG iteration 176, residual = 2.039468653013079e-06     

CG iteration 177, residual = 1.7241735962693453e-06     

CG iteration 178, residual = 3.6593111967828437e-06     

CG iteration 179, residual = 7.697206769663619e-06     

CG iteration 180, residual = 1.9330236949495707e-06     

CG iteration 181, residual = 6.2852042798801e-07     

CG iteration 182, residual = 5.612494046032096e-07     

CG iteration 183, residual = 1.1075374528180186e-06     

CG iteration 184, residual = 1.2049147514362113e-06     

CG iteration 185, residual = 9.398440593779961e-07     

CG iteration 186, residual = 3.7980458358672633e-07     

CG iteration 187, residual = 2.5670523296503786e-07     

CG iteration 188, residual = 6.084088449139637e-07     

CG iteration 189, residual = 3.983694106613975e-07     

CG iteration 190, residual = 8.363100169294131e-07     

CG iteration 191, residual = 3.310925435285469e-07     

CG iteration 192, residual = 4.731034701404834e-07     

CG iteration 193, residual = 1.806126618630165e-06     

CG iteration 194, residual = 8.721849041932642e-07     

CG iteration 195, residual = 1.027003004516983e-06     

CG iteration 196, residual = 4.3657078689953596e-07     

CG iteration 197, residual = 1.5700410821813838e-06     

CG iteration 198, residual = 1.3644325724038912e-06     

CG iteration 199, residual = 8.841636419589532e-07     

CG iteration 200, residual = 5.618200349470977e-07     

CG iteration 201, residual = 5.6118875825585905e-08     

CG iteration 202, residual = 1.075763119871179e-07     

CG iteration 203, residual = 3.7867124772891555e-08     

CG iteration 204, residual = 3.755452404676314e-08     

CG iteration 205, residual = 8.890158386447584e-08     

CG iteration 206, residual = 4.4655064046432975e-08     

CG iteration 207, residual = 2.680927532310563e-07     

CG iteration 208, residual = 6.949609081783606e-08     

CG iteration 209, residual = 1.5829720724994856e-07     

CG iteration 210, residual = 1.331901320327998e-07     

CG iteration 211, residual = 1.6594524703698698e-07     

CG iteration 212, residual = 3.8369362040272903e-07     

CG iteration 213, residual = 7.146615462210559e-08     

CG iteration 214, residual = 6.507041434222231e-08     

CG iteration 215, residual = 6.17405563985919e-08     

CG iteration 216, residual = 5.6322806663179336e-08     

CG iteration 217, residual = 3.4729362616092196e-08     

CG iteration 218, residual = 6.018053077986717e-08     

CG iteration 219, residual = 1.641805705872491e-07     

CG iteration 220, residual = 3.454875596629985e-08     

CG iteration 221, residual = 5.413604536880722e-08     

CG iteration 222, residual = 1.4811447600133004e-08     

CG iteration 223, residual = 1.6554447376295076e-08     

CG iteration 224, residual = 2.7010189541250416e-08     

CG iteration 225, residual = 1.0446490347240108e-08     

CG iteration 226, residual = 4.308989495579487e-08     

CG iteration 227, residual = 6.837619794838597e-08     

CG iteration 228, residual = 7.866579935652998e-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.