Maxwell PEC Direct method

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

\[\begin{split} \begin{array}{rcl}(1) \quad \quad \quad \boldsymbol E(x) &=& \underbrace{ \kappa \, \displaystyle \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)} \\[2ex] && \quad + \;\underbrace{ \nabla \times \int\limits_\Gamma \displaystyle{ \frac{1}{4\,\pi} \, \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \, \boldsymbol n \times \boldsymbol{m}(y)\, \mathrm{d}s_y }_{\displaystyle \mathrm{DL}(\boldsymbol n \times\boldsymbol m)}\,. \end{array}\end{split}\]

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

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

Consider as incidence an incoming plane wave and compute the given 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 = GridFunction(fesHCurl) 
m.Set( -E_inc, definedon=mesh.Boundaries(".*"), dual=True) # Hcurl

kappa =  5.0034083602476045

Boundary Integral Equation

We carefully apply the tangential trace to (1) and obtain 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 + \displaystyle \int\limits_\Gamma \mathrm{DL}(\boldsymbol n \times \boldsymbol{m}) \cdot \boldsymbol v \, \mathrm d s \]

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

\[ \mathrm{V} \, \mathrm j = \left( \frac12 \,\mathrm{M} - \mathrm{K} \right)\, \mathrm{m}\,,\]

where

  • \(\mathrm V\) is the Maxwell single layer operator.

  • \(\mathrm K\) is the Maxwell double layer operator.

  • \(\mathrm{M}\) is the mass matrix.

The Maxwell single and double layer potential operators are derivates of Helmholtz operators (details are given below) and thus, in NGSBEM we use the Helmholtz operators as building blocks for Maxwell operators:

V = kappa * HelmholtzSL(u*ds, kappa)*v*ds
    - 1/kappa * HelmholtzSL(div(u)*ds, kappa)*div(v)*ds
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=False)

Draw(j[1].real, mesh, draw_vol=False, order=3);

Details: Explicit Representation of the Maxwell Layer Potentials

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

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

\[ K_{ij} = \int\limits_\Gamma \underbrace{\displaystyle\int\limits_\Gamma \displaystyle{ \frac{1}{4\,\pi} \, \nabla_x \frac{e^{i\,\kappa\,\|x-y\|}}{\| x-y\|} } \times ( \boldsymbol n \times \boldsymbol{u}_j(y) ) \, \mathrm{d} s_y }_{ \displaystyle{\mathrm{HelmholtzDL}( \boldsymbol n \times \boldsymbol u_j) }}\cdot \boldsymbol v_i(x) \mathrm d s_x\,.\]

References (theoretical and numerical results):