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

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 potential 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

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 Solution

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 (note that the potential evaluation does not yet evaluate the full SL)
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 (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)
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);

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.