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

Boundary Value Problem

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}\]

Representation Formula

A possible representation for \(\boldsymbol E\) reads

\[ \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)}} \,.\]

Solution with NGSolve

Define the geometry of the perfect conductor 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

Boundary Integral Equation

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

\[\begin{split} \begin{array}{l c c c c} \forall \, \boldsymbol v\in H^{-\frac12}(\mathrm{div}_\Gamma, \Gamma): & & \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 \\[2ex] & \Rightarrow & \mathrm{V}\,\mathrm j &=& \mathrm{M}\, \mathrm m \,. \end{array} \end{split}\]

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.

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

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

where the single layer operator \(V\) on the left is a regular and symmetric matrix. Its assembly 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

\[ V_{ij} = \kappa \, \int\limits_\Gamma \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 \mathrm{d} s_x - \frac{1}{\kappa} \int\limits_\Gamma \int\limits_\Gamma \displaystyle{ \frac{1}{4\,\pi}\, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \mathrm{div}_\Gamma \boldsymbol u_i(y)\, \mathrm{div}_\Gamma \boldsymbol v_i(x)\, \mathrm{d}\sigma_y\,\mathrm d s_x\,. \]

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

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

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

j = GridFunction(fesHDiv)
pre = BilinearForm( uHDiv.Trace() * vHDiv.Trace() *ds ).Assemble().mat.Inverse(freedofs=fesHDiv.FreeDofs()) 
with TaskManager():
    CG(mat = V, 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);

Evaluation of the Scattered Field

Define the geometry of the receiver, generate a mesh and evaluate the scattered field \(\boldsymbol E\).

• consider as receiver a sphere around the PEC body:

# sphere 
order = 3
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)

V = MaxwellSingleLayerPotentialOperator(fesHDiv, kappa, intorder=8 + 2 * order)
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);

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

# 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);

References

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