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

CG iteration 2, residual = 0.7675484948921725     

CG iteration 3, residual = 0.6522908847399106     

CG iteration 4, residual = 0.5438772154771251     

CG iteration 5, residual = 0.3257757974212374     

CG iteration 6, residual = 0.3828518205783512     

CG iteration 7, residual = 0.46451571786004864     

CG iteration 8, residual = 0.5016221047913626     

CG iteration 9, residual = 0.48290686157949725     

CG iteration 10, residual = 0.508728959042658     

CG iteration 11, residual = 0.31246843950220715     

CG iteration 12, residual = 0.29831571689100556     

CG iteration 13, residual = 1.1125732551757428     

CG iteration 14, residual = 0.3039874045704071     

CG iteration 15, residual = 0.09063205902538926     

CG iteration 16, residual = 0.07227404776234742     

CG iteration 17, residual = 0.37870365032132386     

CG iteration 18, residual = 0.06129680510639545     

CG iteration 19, residual = 0.12918056268691255     

CG iteration 20, residual = 0.21844979153596294     

CG iteration 21, residual = 0.08357025497080074     

CG iteration 22, residual = 0.09975574237281183     

CG iteration 23, residual = 0.08991043125014153     

CG iteration 24, residual = 0.2478338879683376     

CG iteration 25, residual = 0.07719944563793023     

CG iteration 26, residual = 0.07468708638611218     

CG iteration 27, residual = 0.04013173980426076     

CG iteration 28, residual = 0.012430995371741926     

CG iteration 29, residual = 0.011834162915436595     

CG iteration 30, residual = 0.015615968412647762     

CG iteration 31, residual = 0.08045099702397933     

CG iteration 32, residual = 0.16706771842989868     

CG iteration 33, residual = 0.2742261091389619     

CG iteration 34, residual = 0.17482960602654074     

CG iteration 35, residual = 0.16658275177362245     

CG iteration 36, residual = 0.11586428756141803     

CG iteration 37, residual = 0.14497092246849158     

CG iteration 38, residual = 0.3733772697290971     

CG iteration 39, residual = 0.014658401141594139     

CG iteration 40, residual = 0.009788582577771803     

CG iteration 41, residual = 0.020692061805050965     

CG iteration 42, residual = 0.012942468692453099     

CG iteration 43, residual = 0.014620752826879508     

CG iteration 44, residual = 0.018974070980016574     

CG iteration 45, residual = 0.03433029321570887     

CG iteration 46, residual = 0.12758184280329118     

CG iteration 47, residual = 0.008621139343384451     

CG iteration 48, residual = 0.00919863693928698     

CG iteration 49, residual = 0.010499384520971785     

CG iteration 50, residual = 0.0835531456032219     

CG iteration 51, residual = 0.00788180039388755     

CG iteration 52, residual = 0.022536712039589516     

CG iteration 53, residual = 0.0018846561801928064     

CG iteration 54, residual = 0.0017049216885456967     

CG iteration 55, residual = 0.002949689330026594     

CG iteration 56, residual = 0.0020120928596884077     

CG iteration 57, residual = 0.003952298039004946     

CG iteration 58, residual = 0.005271102095278247     

CG iteration 59, residual = 0.0012179963909683877     

CG iteration 60, residual = 0.0011483472238320532     

CG iteration 61, residual = 0.0013840424595872008     

CG iteration 62, residual = 0.0016143550394456829     

CG iteration 63, residual = 0.005466806835382151     

CG iteration 64, residual = 0.004181173129576876     

CG iteration 65, residual = 0.0024137312707454136     

CG iteration 66, residual = 0.0015277024552144275     

CG iteration 67, residual = 0.0004880507254271947     

CG iteration 68, residual = 0.00029577159853513396     

CG iteration 69, residual = 0.00045498397582712896     

CG iteration 70, residual = 0.0008468523798601628     

CG iteration 71, residual = 0.0027340718179700366     

CG iteration 72, residual = 0.0007395697225260665     

CG iteration 73, residual = 0.0005542792391691084     

CG iteration 74, residual = 0.0012025103228978803     

CG iteration 75, residual = 0.0007935879593495375     

CG iteration 76, residual = 0.0008796495835036756     

CG iteration 77, residual = 0.0003582060732353074     

CG iteration 78, residual = 0.00035474341965780035     

CG iteration 79, residual = 0.00011369013940790966     

CG iteration 80, residual = 0.0001108242835509134     

CG iteration 81, residual = 0.00026273364388183616     

CG iteration 82, residual = 0.000602104846986161     

CG iteration 83, residual = 0.018886684691044528     

CG iteration 84, residual = 0.00028849708417135546     

CG iteration 85, residual = 0.0001379319659754317     

CG iteration 86, residual = 0.00011373992415663969     

CG iteration 87, residual = 0.00016861540650487777     

CG iteration 88, residual = 0.00021525045996288508     

CG iteration 89, residual = 0.0002514389866132597     

CG iteration 90, residual = 0.00015269796019450884     

CG iteration 91, residual = 0.0008578520446322545     

CG iteration 92, residual = 0.0003115578411544488     

CG iteration 93, residual = 0.0001741316553337495     

CG iteration 94, residual = 0.00019951835883993103     

CG iteration 95, residual = 0.00020214589455647323     

CG iteration 96, residual = 0.0001063144416258087     

CG iteration 97, residual = 9.29874908261698e-05     

CG iteration 98, residual = 0.00016231989836625274     

CG iteration 99, residual = 0.00014746638098150756     

CG iteration 100, residual = 8.43993460451938e-05     

CG iteration 101, residual = 0.00011974224747985178     

CG iteration 102, residual = 0.00013414002967049876     

CG iteration 103, residual = 4.2585165176743254e-05     

CG iteration 104, residual = 4.134792664171689e-05     

CG iteration 105, residual = 4.070038520291852e-05     

CG iteration 106, residual = 7.826568368018298e-05     

CG iteration 107, residual = 3.737525155413043e-05     

CG iteration 108, residual = 2.210500976527113e-05     

CG iteration 109, residual = 3.114926757589349e-05     

CG iteration 110, residual = 0.00010747456277509236     

CG iteration 111, residual = 7.165188342281352e-05     

CG iteration 112, residual = 4.154560570411353e-05     

CG iteration 113, residual = 9.810597140909305e-06     

CG iteration 114, residual = 1.9397399947415383e-05     

CG iteration 115, residual = 1.4438666785627173e-05     

CG iteration 116, residual = 2.1095984936138085e-05     

CG iteration 117, residual = 2.694609856444525e-05     

CG iteration 118, residual = 2.0979915894302752e-05     

CG iteration 119, residual = 4.856117485661302e-05     

CG iteration 120, residual = 2.8953145319404095e-05     

CG iteration 121, residual = 3.162908906620334e-05     

CG iteration 122, residual = 3.936152621918116e-05     

CG iteration 123, residual = 3.163505004765991e-05     

CG iteration 124, residual = 1.3925859948707997e-05     

CG iteration 125, residual = 6.128690155594763e-06     

CG iteration 126, residual = 1.1474840837005365e-05     

CG iteration 127, residual = 7.270720753240743e-06     

CG iteration 128, residual = 1.1530735870345444e-05     

CG iteration 129, residual = 7.229900874544562e-05     

CG iteration 130, residual = 4.945952064974844e-06     

CG iteration 131, residual = 2.757417235800297e-06     

CG iteration 132, residual = 4.9577542138908915e-06     

CG iteration 133, residual = 1.6956809362238109e-06     

CG iteration 134, residual = 1.3632654709851112e-06     

CG iteration 135, residual = 2.144743673898861e-06     

CG iteration 136, residual = 4.212037860707773e-06     

CG iteration 137, residual = 7.891421042252774e-06     

CG iteration 138, residual = 1.9170343696706787e-06     

CG iteration 139, residual = 6.8185106339528925e-06     

CG iteration 140, residual = 1.6934600399349452e-06     

CG iteration 141, residual = 1.8041487389312564e-06     

CG iteration 142, residual = 2.6494524218985833e-06     

CG iteration 143, residual = 2.137154168576765e-06     

CG iteration 144, residual = 3.3344822035607256e-06     

CG iteration 145, residual = 3.153477479113954e-06     

CG iteration 146, residual = 7.261638693621376e-06     

CG iteration 147, residual = 3.240959083666438e-06     

CG iteration 148, residual = 5.734503547594529e-06     

CG iteration 149, residual = 3.1017054759510646e-06     

CG iteration 150, residual = 3.502731842010674e-06     

CG iteration 151, residual = 1.5608914466437834e-06     

CG iteration 152, residual = 4.334589926873434e-06     

CG iteration 153, residual = 1.38068739674066e-06     

CG iteration 154, residual = 4.794618771082254e-06     

CG iteration 155, residual = 1.3528313752681319e-06     

CG iteration 156, residual = 1.6195081805467902e-06     

CG iteration 157, residual = 8.119964751901109e-07     

CG iteration 158, residual = 2.901165960425894e-06     

CG iteration 159, residual = 4.499810708714811e-06     

CG iteration 160, residual = 3.4299594611611387e-06     

CG iteration 161, residual = 4.010135903416973e-06     

CG iteration 162, residual = 1.2021425587674848e-06     

CG iteration 163, residual = 1.8225350938966835e-06     

CG iteration 164, residual = 1.4285985855667896e-06     

CG iteration 165, residual = 1.687379850791429e-06     

CG iteration 166, residual = 4.285955134391863e-06     

CG iteration 167, residual = 3.1143717768366783e-06     

CG iteration 168, residual = 7.158013149403422e-06     

CG iteration 169, residual = 2.671197890194019e-06     

CG iteration 170, residual = 2.930624311130495e-06     

CG iteration 171, residual = 3.7827531388511606e-06     

CG iteration 172, residual = 3.4785939971851688e-06     

CG iteration 173, residual = 4.647737968491314e-06     

CG iteration 174, residual = 3.173825883662576e-06     

CG iteration 175, residual = 4.269020716564197e-06     

CG iteration 176, residual = 1.965011400883924e-06     

CG iteration 177, residual = 1.7424944249969577e-06     

CG iteration 178, residual = 3.943902199317897e-06     

CG iteration 179, residual = 4.45798775592135e-06     

CG iteration 180, residual = 1.9394716804116428e-06     

CG iteration 181, residual = 8.25952487013108e-07     

CG iteration 182, residual = 7.407196257709326e-07     

CG iteration 183, residual = 9.119352301559292e-07     

CG iteration 184, residual = 1.032757770036132e-06     

CG iteration 185, residual = 9.344962122439767e-07     

CG iteration 186, residual = 4.810692702151693e-07     

CG iteration 187, residual = 2.6984308369092647e-07     

CG iteration 188, residual = 1.212556637874085e-06     

CG iteration 189, residual = 5.202993469019153e-07     

CG iteration 190, residual = 4.226653829467434e-07     

CG iteration 191, residual = 4.926266247721214e-07     

CG iteration 192, residual = 2.965772631592379e-07     

CG iteration 193, residual = 1.9549332552075186e-07     

CG iteration 194, residual = 2.8581346956606045e-07     

CG iteration 195, residual = 1.2036314114635542e-06     

CG iteration 196, residual = 1.473748851869797e-06     

CG iteration 197, residual = 4.5594174815750263e-07     

CG iteration 198, residual = 3.330905539995466e-07     

CG iteration 199, residual = 6.087064757537714e-07     

CG iteration 200, residual = 4.085566532141036e-07     

CG iteration 201, residual = 2.5882134118511506e-07     

CG iteration 202, residual = 6.496047347408687e-08     

CG iteration 203, residual = 1.6156293934341624e-07     

CG iteration 204, residual = 9.253814317575873e-08     

CG iteration 205, residual = 1.1893572849631256e-07     

CG iteration 206, residual = 8.479410056902864e-08     

CG iteration 207, residual = 1.3201022247356923e-07     

CG iteration 208, residual = 5.3453487001086745e-08     

CG iteration 209, residual = 7.799493426449284e-08     

CG iteration 210, residual = 7.515853736056125e-08     

CG iteration 211, residual = 2.351123841779232e-07     

CG iteration 212, residual = 2.3545739627857754e-07     

CG iteration 213, residual = 1.1494145213529988e-07     

CG iteration 214, residual = 1.1776544513152334e-07     

CG iteration 215, residual = 4.147183895835142e-08     

CG iteration 216, residual = 4.240922889251228e-08     

CG iteration 217, residual = 4.380045850858469e-08     

CG iteration 218, residual = 8.183761786643345e-08     

CG iteration 219, residual = 4.0254760958949176e-08     

CG iteration 220, residual = 4.4540218685145354e-08     

CG iteration 221, residual = 1.107013587687246e-07     

CG iteration 222, residual = 2.5387251935130698e-08     

CG iteration 223, residual = 2.3987700056118546e-08     

CG iteration 224, residual = 1.5105761904675412e-08     

CG iteration 225, residual = 1.0225551762495869e-08     

CG iteration 226, residual = 2.6908927823366423e-08     

CG iteration 227, residual = 1.318413784963539e-08     

CG iteration 228, residual = 1.1422500606963993e-08     

CG iteration 229, residual = 2.1748918299869495e-08     

CG iteration 230, residual = 2.382423933184447e-08     

CG iteration 231, residual = 1.9085039648145576e-08     

CG iteration 232, residual = 2.7130744126992957e-08     

CG iteration 233, residual = 1.503762046432712e-08     

CG iteration 234, residual = 5.134807688298857e-08     

CG iteration 235, residual = 7.739079440682652e-08     

CG iteration 236, residual = 2.1460932186319042e-08     

CG iteration 237, residual = 1.5611845095367398e-08     

CG iteration 238, residual = 1.4125255348047326e-08     

CG iteration 239, residual = 1.0015222506633392e-08     

CG iteration 240, residual = 3.770767551731907e-08     

CG iteration 241, residual = 8.145855052349764e-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.