Maxwell solver for Pec using direct formulation

Maxwell solver for Pec using direct formulation#

keys: Maxwell double layer potential, PEC scattering, MoM, Neumann trace, surface current density

from netgen.occ import *
import netgen.meshing as meshing
from ngsolve import *
from ngsolve.webgui import Draw
from libbem import *
from ngsolve import Projector, Preconditioner
from ngsolve.krylovspace import CG, GMRes

Loading ngbem library

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:

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

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}\sigma_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}\sigma_y }_{\displaystyle{\mathrm{SL}(\boldsymbol j)} } + \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}\sigma_y }_{\displaystyle{ \mathrm{DL} (\boldsymbol n \times \boldsymbol{m}) } } \,.\]

We carefully apply the tangential trace $\gamma_R$ on $\boldsymbol E$ and obtain a boundary integral equation for the unknown trace $\boldsymbol j$. Note that $\boldsymbol m \in H^{-\frac12}(\mathrm{curl}_\Gamma, \Gamma)$ is the given Dirichlet data. The boundary integral equation is solved by the boundary element method, i.e. the numerical solution of the variational formulation

\[ \forall \, \boldsymbol v\in H^{-\frac12}(\mathrm{div}_\Gamma, \Gamma): \quad \left\langle \mathrm{SL} (\boldsymbol j), \boldsymbol v \right\rangle_{-\frac12} = \left\langle \boldsymbol m, \boldsymbol v\right\rangle_{-\frac12} - \left\langle \mathrm{DL}(\boldsymbol n \times \boldsymbol{m}), \boldsymbol v\right\rangle_{-\frac12} \,. \]

In the enineering community, the approximation scheme is also known as method of moments (MoM).

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

# V: Hdiv is TnT; K: Hdiv as test and HCurl as trial space with rotated evaluator   
with TaskManager(): 
    V = MaxwellSingleLayerPotentialOperator(fesHDiv, kappa, 
                                            intorder=10, leafsize=40, eta=3., eps=1e-4, method="aca", testhmatrix=False)
    K = MaxwellDoubleLayerPotentialOperator(fesHCurl, fesHDiv, kappa, 
                                            intorder=12, leafsize=40, eta=3., eps=1e-6, method="aca", testhmatrix=False)

# 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(): 
    M = BilinearForm( uHCurl.Trace() * vHDiv.Trace()* ds(bonus_intorder=3)).Assemble() # <Hcurl, Hdiv>  
    rhs = ( (0.5 * M.mat - K.mat ) * m.vec).Evaluate() 
    #GMRes(A = V.mat, pre=pre, b = rhs, x=j.vec, tol=1e-8, maxsteps=5000, printrates=True)
    CG(mat = V.mat, pre=pre, rhs = rhs, sol=j.vec, tol=1e-8, maxsteps=500, initialize=False, printrates=False)

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