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, GMRes
keys: PEC scattering, single and double layer potential, MoM, unknown Neumann trace (surface current density)
Maxwell PEC Direct method#
We consider a perfect conductor \(\Omega \subset \mathbb R^3\) and an electromagnetic signal that induces a scattered electric field which propagates into \(\Omega^c\). The scattered electric field \(\boldsymbol E\) solves the following Dirichlet boundary value problem:
A possible representation for \(\boldsymbol E\) reads
NGSolve Solution
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(3)
Next, we create finite element spaces for \(\boldsymbol H^{-\frac12}(\mathrm{curl}_\Gamma, \Gamma)\) (test functions) and \(\boldsymbol H^{-\frac12}(\mathrm{div}_\Gamma, \Gamma)\) (trial functions):
order = 3
fesHCurl = HCurl(mesh, order=order, complex=True)
fesHDiv = HDivSurface(mesh, order=order, 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
Boundary Integral Equation
We carefully apply the tangential trace \(\gamma_R\) on \(\boldsymbol E\) and obtain a boundary integral equation for the unknown trace \(\boldsymbol j\). The boundary integral equation is solved by the boundary element method, i.e. the numerical solution of the variational formulation
In the enineering community, the approximation scheme is also known as method of moments (MoM).
Define the incoming plane wave and compute the given boundary data \(\boldsymbol m = -\gamma_R \boldsymbol E_{\mathrm{inc}}\) :
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 = GridFunction(fesHCurl)
m.Set( -E_inc, definedon=mesh.Boundaries(".*"), dual=True) # Hcurl
kappa = 5.0034083602476045
The discretisation of the above variational formulation leads to a system of linear equations, ie $\( \mathrm{V} \, \mathrm j = \left( \frac12 \,\mathrm{M} - \mathrm{K} \right)\, \mathrm{m}\,,\)$ where
\(\mathrm V\) is the Maxwell single layer operator. \(V\) is a regular, symmetric matrix.
\(\mathrm{M}\) is the mass matrix.
\(\mathrm K\) is the Maxwell double layer operator. \(K\) is quadratic:
The assembly of the potential operators is declarative and close to their actual implementations. All contributions to the Maxwell layer potentials are derivates of the Helmholtz layer potentials and thus, we use the Helmholtz operators as building blocks.
Details
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
The two contributions are derivates of the Helmholtz single layer potential operator and may be implemented as such:
V = kappa * (HelmholtzSL(u*ds, kappa)*v*ds - 1/kappa * HelmholtzSL(div(u)*ds, kappa)*div(v)*ds
For a trial function \(\boldsymbol u_j \in H^{-\frac12}(\mathrm{curl}_\Gamma,\Gamma)\) and a test function \(\boldsymbol v_i \in H^{-\frac12}(\mathrm{div}_\Gamma,\Gamma)\) the Maxwell double layer potential operator entry \(K_{ij}\) can formally written as
Giving transparency to its relation to the Helmholtz double layer potential operator, you may assemble it as follows:
K = HelmholtzDL( u.Operator("rotated_trace") * ds, kappa) * v * ds
intorder = 2 * order + 6
with TaskManager():
M = BilinearForm( uHCurl.Trace()*vHDiv.Trace()*ds(bonus_intorder=3)).Assemble() # <Hcurl, Hdiv>
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
K = HelmholtzDL(uHCurl.Operator("rotated_trace") * ds(bonus_intorder=6), kappa) * vHDiv.Trace() * ds
# solve direct formulation for j, i.e., Neumann data of scattered electric field E
j = GridFunction(fesHDiv)
pre = BilinearForm( uHDiv.Trace() * vHDiv.Trace() *ds).Assemble().mat.Inverse(freedofs=fesHDiv.FreeDofs())
with TaskManager():
rhs = ( (0.5 * M.mat - K.mat ) * m.vec).Evaluate()
GMRes(A = V, pre=pre, b = rhs, x=j.vec, tol=1e-8, maxsteps=5000, printrates=True)
GMRes iteration 1, residual = 5.875912781730529
GMRes iteration 2, residual = 3.895835799831684
GMRes iteration 3, residual = 3.4282185918274277
GMRes iteration 4, residual = 2.929630898281488
GMRes iteration 5, residual = 2.5206022969868114
GMRes iteration 6, residual = 2.226877635311094
GMRes iteration 7, residual = 1.7019221835989653
GMRes iteration 8, residual = 1.5120679099879302
GMRes iteration 9, residual = 1.4524220108643635
GMRes iteration 10, residual = 1.2479901603252486
GMRes iteration 11, residual = 1.0930969688820054
GMRes iteration 12, residual = 0.9882770233858005
GMRes iteration 13, residual = 0.7870043002718301
GMRes iteration 14, residual = 0.6221401838875233
GMRes iteration 15, residual = 0.5914503367566595
GMRes iteration 16, residual = 0.5638398602421663
GMRes iteration 17, residual = 0.42960188925807286
GMRes iteration 18, residual = 0.3209168882275205
GMRes iteration 19, residual = 0.2891982644417119
GMRes iteration 20, residual = 0.2799262353738335
GMRes iteration 21, residual = 0.2212017993039132
GMRes iteration 22, residual = 0.19259890970425872
GMRes iteration 23, residual = 0.180538524145343
GMRes iteration 24, residual = 0.15799732740064976
GMRes iteration 25, residual = 0.15308371735086684
GMRes iteration 26, residual = 0.14580605602953603
GMRes iteration 27, residual = 0.14150894669009723
GMRes iteration 28, residual = 0.1134049285029418
GMRes iteration 29, residual = 0.09387117927904735
GMRes iteration 30, residual = 0.09334128307924996
GMRes iteration 31, residual = 0.08274496785641715
GMRes iteration 32, residual = 0.07578010881472153
GMRes iteration 33, residual = 0.0754967614745199
GMRes iteration 34, residual = 0.0646960305223182
GMRes iteration 35, residual = 0.05111578158201664
GMRes iteration 36, residual = 0.04995982556483856
GMRes iteration 37, residual = 0.04536842850060138
GMRes iteration 38, residual = 0.03548716200117304
GMRes iteration 39, residual = 0.035073423793167906
GMRes iteration 40, residual = 0.03047940878465369
GMRes iteration 41, residual = 0.0297873074317586
GMRes iteration 42, residual = 0.024783357428607978
GMRes iteration 43, residual = 0.020982709735705853
GMRes iteration 44, residual = 0.0194420994378084
GMRes iteration 45, residual = 0.018521842571608528
GMRes iteration 46, residual = 0.017833204362742196
GMRes iteration 47, residual = 0.01682181536354052
GMRes iteration 48, residual = 0.01632577408488297
GMRes iteration 49, residual = 0.014317035106190928
GMRes iteration 50, residual = 0.01267842371518941
GMRes iteration 51, residual = 0.011981029843644046
GMRes iteration 52, residual = 0.010967706618433241
GMRes iteration 53, residual = 0.009309156931245
GMRes iteration 54, residual = 0.009288987773763156
GMRes iteration 55, residual = 0.008149129583021247
GMRes iteration 56, residual = 0.007035833406100308
GMRes iteration 57, residual = 0.007017532164598678
GMRes iteration 58, residual = 0.0065337705392920055
GMRes iteration 59, residual = 0.006399350334820928
GMRes iteration 60, residual = 0.006244976226740058
GMRes iteration 61, residual = 0.006036746267875387
GMRes iteration 62, residual = 0.005968553466414637
GMRes iteration 63, residual = 0.005870004038640157
GMRes iteration 64, residual = 0.005396458991840618
GMRes iteration 65, residual = 0.005214591759801279
GMRes iteration 66, residual = 0.004971150900201493
GMRes iteration 67, residual = 0.00493420599160983
GMRes iteration 68, residual = 0.004680312998550088
GMRes iteration 69, residual = 0.0041821908903025415
GMRes iteration 70, residual = 0.00408309957021825
GMRes iteration 71, residual = 0.004015118800209142
GMRes iteration 72, residual = 0.0038440097600087315
GMRes iteration 73, residual = 0.003809953833034546
GMRes iteration 74, residual = 0.0036290847809485327
GMRes iteration 75, residual = 0.0029954558509393886
GMRes iteration 76, residual = 0.002798211535613743
GMRes iteration 77, residual = 0.0027193186349833056
GMRes iteration 78, residual = 0.0024585449090583802
GMRes iteration 79, residual = 0.002432679781440819
GMRes iteration 80, residual = 0.0021657149968882894
GMRes iteration 81, residual = 0.002069593795126662
GMRes iteration 82, residual = 0.0020205837270987808
GMRes iteration 83, residual = 0.0018028678628945098
GMRes iteration 84, residual = 0.0017599187845828278
GMRes iteration 85, residual = 0.0015945154456765227
GMRes iteration 86, residual = 0.0014591801843920726
GMRes iteration 87, residual = 0.001299579110272271
GMRes iteration 88, residual = 0.0011880010780607176
GMRes iteration 89, residual = 0.001182086051278911
GMRes iteration 90, residual = 0.0010697657201590351
GMRes iteration 91, residual = 0.0009793526024751487
GMRes iteration 92, residual = 0.0009520396697023597
GMRes iteration 93, residual = 0.0009436223371359537
GMRes iteration 94, residual = 0.0008999431845302399
GMRes iteration 95, residual = 0.0008785231164369504
GMRes iteration 96, residual = 0.0008703378728335613
GMRes iteration 97, residual = 0.0007991368103559334
GMRes iteration 98, residual = 0.0007700442633720919
GMRes iteration 99, residual = 0.000746469335141632
GMRes iteration 100, residual = 0.0007173064920481095
GMRes iteration 101, residual = 0.000688535291023884
GMRes iteration 102, residual = 0.0005802684023144156
GMRes iteration 103, residual = 0.0005629195440887707
GMRes iteration 104, residual = 0.0005524669143940528
GMRes iteration 105, residual = 0.0005029512017530035
GMRes iteration 106, residual = 0.0004967056146980039
GMRes iteration 107, residual = 0.0004859902458742612
GMRes iteration 108, residual = 0.00047891912364841455
GMRes iteration 109, residual = 0.0004175700741425474
GMRes iteration 110, residual = 0.00038661735540090143
GMRes iteration 111, residual = 0.0003866117157685571
GMRes iteration 112, residual = 0.00035858370573274296
GMRes iteration 113, residual = 0.0003530736489263808
GMRes iteration 114, residual = 0.0003276343941135416
GMRes iteration 115, residual = 0.0003176175284724877
GMRes iteration 116, residual = 0.00030288182503549393
GMRes iteration 117, residual = 0.00029133867128002113
GMRes iteration 118, residual = 0.00028826255373812405
GMRes iteration 119, residual = 0.00028288112787082956
GMRes iteration 120, residual = 0.00028032698450459895
GMRes iteration 121, residual = 0.0002677949170141268
GMRes iteration 122, residual = 0.00025414373980248205
GMRes iteration 123, residual = 0.00025157624858484655
GMRes iteration 124, residual = 0.0002459893082780905
GMRes iteration 125, residual = 0.00023639095364446307
GMRes iteration 126, residual = 0.00023155078026727025
GMRes iteration 127, residual = 0.0002283405614754116
GMRes iteration 128, residual = 0.00022704805907451192
GMRes iteration 129, residual = 0.000212388675129082
GMRes iteration 130, residual = 0.00019447197586166886
GMRes iteration 131, residual = 0.0001861725191606631
GMRes iteration 132, residual = 0.0001848331812886134
GMRes iteration 133, residual = 0.00017724305525666724
GMRes iteration 134, residual = 0.00016882838864563433
GMRes iteration 135, residual = 0.000167944802562362
GMRes iteration 136, residual = 0.00016090213712957813
GMRes iteration 137, residual = 0.00015953638564446642
GMRes iteration 138, residual = 0.0001499743731936391
GMRes iteration 139, residual = 0.00014638040715799375
GMRes iteration 140, residual = 0.0001361879162593781
GMRes iteration 141, residual = 0.00011252153374050801
GMRes iteration 142, residual = 0.00010617180270402915
GMRes iteration 143, residual = 0.00010315255458024677
GMRes iteration 144, residual = 8.531737471014759e-05
GMRes iteration 145, residual = 8.150429090115636e-05
GMRes iteration 146, residual = 8.016371916699445e-05
GMRes iteration 147, residual = 6.505360842527056e-05
GMRes iteration 148, residual = 5.721760292549816e-05
GMRes iteration 149, residual = 5.640397178624766e-05
GMRes iteration 150, residual = 5.2034430372560116e-05
GMRes iteration 151, residual = 5.133427981651785e-05
GMRes iteration 152, residual = 4.901457425944337e-05
GMRes iteration 153, residual = 4.494845008502004e-05
GMRes iteration 154, residual = 4.476958900203222e-05
GMRes iteration 155, residual = 4.237311111608042e-05
GMRes iteration 156, residual = 4.202699667437448e-05
GMRes iteration 157, residual = 3.855313388278495e-05
GMRes iteration 158, residual = 3.3871352177276296e-05
GMRes iteration 159, residual = 3.381671701441948e-05
GMRes iteration 160, residual = 3.0624329777488964e-05
GMRes iteration 161, residual = 2.8878360259979553e-05
GMRes iteration 162, residual = 2.6875461222728113e-05
GMRes iteration 163, residual = 1.7851820628116935e-05
GMRes iteration 164, residual = 1.5915446874566823e-05
GMRes iteration 165, residual = 1.5314034380429553e-05
GMRes iteration 166, residual = 1.0186814278645395e-05
GMRes iteration 167, residual = 8.88612403459561e-06
GMRes iteration 168, residual = 8.736347132262719e-06
GMRes iteration 169, residual = 6.781907378404527e-06
GMRes iteration 170, residual = 5.658910178871896e-06
GMRes iteration 171, residual = 5.395132995698412e-06
GMRes iteration 172, residual = 5.291943673007458e-06
GMRes iteration 173, residual = 4.6911461431183416e-06
GMRes iteration 174, residual = 4.478941808656404e-06
GMRes iteration 175, residual = 4.462801189496685e-06
GMRes iteration 176, residual = 4.194299342011604e-06
GMRes iteration 177, residual = 4.0021147428679e-06
GMRes iteration 178, residual = 3.9316960013104005e-06
GMRes iteration 179, residual = 3.6387533389867317e-06
GMRes iteration 180, residual = 3.559488626641445e-06
GMRes iteration 181, residual = 3.4588486409531523e-06
GMRes iteration 182, residual = 3.0208703220289396e-06
GMRes iteration 183, residual = 2.390532828058634e-06
GMRes iteration 184, residual = 2.258039617380226e-06
GMRes iteration 185, residual = 2.2380232793576227e-06
GMRes iteration 186, residual = 1.8791952005091292e-06
GMRes iteration 187, residual = 1.551440348254456e-06
GMRes iteration 188, residual = 1.516630503101824e-06
GMRes iteration 189, residual = 1.4410088767557652e-06
GMRes iteration 190, residual = 1.0790059298565328e-06
GMRes iteration 191, residual = 9.351840748502907e-07
GMRes iteration 192, residual = 9.334951644970278e-07
GMRes iteration 193, residual = 7.870640070854497e-07
GMRes iteration 194, residual = 7.5180550140001e-07
GMRes iteration 195, residual = 6.82969975861739e-07
GMRes iteration 196, residual = 5.829685079048732e-07
GMRes iteration 197, residual = 5.823338713350085e-07
GMRes iteration 198, residual = 5.122747796466715e-07
GMRes iteration 199, residual = 4.722577494054486e-07
GMRes iteration 200, residual = 4.7174777813855475e-07
GMRes iteration 201, residual = 4.416518254314164e-07
GMRes iteration 202, residual = 3.7925900332659293e-07
GMRes iteration 203, residual = 3.5693911637104844e-07
GMRes iteration 204, residual = 3.513716561170613e-07
GMRes iteration 205, residual = 2.7696830570621527e-07
GMRes iteration 206, residual = 2.077943138328978e-07
GMRes iteration 207, residual = 2.0601933391679668e-07
GMRes iteration 208, residual = 1.8031553009777958e-07
GMRes iteration 209, residual = 1.500425165710449e-07
GMRes iteration 210, residual = 1.500018485530224e-07
GMRes iteration 211, residual = 1.3248833932891248e-07
GMRes iteration 212, residual = 1.1633883690330599e-07
GMRes iteration 213, residual = 1.140760389941352e-07
GMRes iteration 214, residual = 1.0677926400222113e-07
GMRes iteration 215, residual = 9.225819903393542e-08
GMRes iteration 216, residual = 8.883795703579575e-08
GMRes iteration 217, residual = 8.722200697256501e-08
GMRes iteration 218, residual = 8.306351030809432e-08
GMRes iteration 219, residual = 7.191546658880918e-08
GMRes iteration 220, residual = 7.120073253203576e-08
GMRes iteration 221, residual = 6.27347867871844e-08
GMRes iteration 222, residual = 5.2025381940705835e-08
GMRes iteration 223, residual = 4.942284759901129e-08
GMRes iteration 224, residual = 4.939376246256086e-08
GMRes iteration 225, residual = 4.5349996267990806e-08
GMRes iteration 226, residual = 3.630784638955651e-08
GMRes iteration 227, residual = 3.4227854411951615e-08
GMRes iteration 228, residual = 3.3873686131320834e-08
GMRes iteration 229, residual = 2.9282237256151098e-08
GMRes iteration 230, residual = 2.6954668540759703e-08
GMRes iteration 231, residual = 2.694686544832879e-08
GMRes iteration 232, residual = 2.4604406700942968e-08
GMRes iteration 233, residual = 2.2479431207793948e-08
GMRes iteration 234, residual = 2.2461550924840198e-08
GMRes iteration 235, residual = 2.130685052028368e-08
GMRes iteration 236, residual = 1.9170173242881253e-08
GMRes iteration 237, residual = 1.909389900425931e-08
GMRes iteration 238, residual = 1.779998760587915e-08
GMRes iteration 239, residual = 1.7016602626756542e-08
GMRes iteration 240, residual = 1.6997793502085226e-08
GMRes iteration 241, residual = 1.593398634138124e-08
GMRes iteration 242, residual = 1.525404615457302e-08
GMRes iteration 243, residual = 1.523813877024295e-08
GMRes iteration 244, residual = 1.4872245733356945e-08
GMRes iteration 245, residual = 1.3699780985750833e-08
GMRes iteration 246, residual = 1.2666666396436208e-08
GMRes iteration 247, residual = 1.2616108263932282e-08
GMRes iteration 248, residual = 1.201160069110699e-08
GMRes iteration 249, residual = 1.0182887894344217e-08
GMRes iteration 250, residual = 9.859050478562424e-09
Draw(j[1].real, mesh, draw_vol=False, order=3);
TODO Evaluate the Solution
To evaluate the solution inside \(\Omega\) we have to evaluate the representation formula, i.e.,
The kernel function that appears in the Maxwell layer potentials is the same as for Helmholtz. Indead, it turns out that the Maxwell layer potentials are derivates of Helmholtz layer potentials. We make this relation clear by using the Helmholtz layer potentials and layer potential operators as building blocks for Maxwell operators. Note that there are formulations such as the low frequency stabilization that even require an explicit split of the summands of the Maxwell single layer potential (References below).
References (theoretical and numerical results):