Maxwell Demo 2

Maxwell Demo 2#

keys: exterior Dirichlet bvp, Maxwell single layer potential, EFIE, PEC scattering

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

Loading ngbem library

Dirichlet Boundary Value Problem		Maxwell Single Layer Potential		Variational Formulation
\( \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. \)	\(\quad \Rightarrow \quad\)	\( \boldsymbol E(x) = \mathrm{SL}(j) \)	\(\quad \Rightarrow \quad\)	\(\left\langle \,\mathrm{SL} (\boldsymbol j),\, \boldsymbol v \right\rangle_{-\frac12} = \left\langle \boldsymbol m , \boldsymbol v\right\rangle_{-\frac12} \)

				\(\mathrm{V} \, \mathrm{j} = \mathrm{M} \, \mathrm{m} \)

NG-BEM Python interface	symbol	FE trial space	FE test space
`MaxwellSingleLayerPotentialOperator`	\(\mathrm V\)	`HDivSurface`	`HDivSurface`
`MaxwellDoubleLayerPotentialOperator`	\(\mathrm K\)	\(\gamma_R\) `HCurl`	`HDivSurface`
`MaxwellSingleLayerPotentialOperatorCurl`	\(\mathrm D\)	\(\gamma_R\) `HCurl`	\(\gamma_R\) `HCurl`
`MaxwellDoublelayerPotentialOperator`	\(\mathrm K'\)	`HDivSurface`	\(\gamma_R\) `HCurl`

Mesh

mesh = Mesh("resources/ship_1.stl")
mesh.Curve(1);

Trail and Test Functions

order = 3
fesHDiv = HDivSurface(mesh, order=order, complex=True)
uHDiv,vHDiv = fesHDiv.TnT() 
fesHCurl = HCurl(mesh, order=order, complex=True)
uHCurl, vHCurl = fesHCurl.TnT() 

Dirichlet Data

a = -0.758
b = -0.453
c = -0.524
E_inc = CF( (-c, 0., a) ) * exp( -1j * (a * x + b * y + c * z))
kappa = sqrt( a*a + b*b + c*c)
n = specialcf.normal(3)
m_exa = -Cross(Cross(n, E_inc), n)
m = GridFunction(fesHCurl)
m.Set(-E_inc, definedon=mesh.Boundaries(".*"), dual=True)
Draw(m[0].real, mesh, draw_vol=False, order=3, min=-0.5, max=0.5);

Right Hand Side

rhs = LinearForm( - E_inc* vHDiv.Trace() * ds(bonus_intorder=20)).Assemble()

System Matrix \(\,\mathrm{V}\)

V = MaxwellSingleLayerPotentialOperator(fesHDiv, kappa, intorder=16, leafsize=100, eta=3, eps=1e-8)

Solve \(\, \mathrm{V} \mathrm j = \mathrm{rhs}\)

j = GridFunction(fesHDiv)
with TaskManager():
    pre = BilinearForm(uHDiv.Trace() * vHDiv.Trace() * ds).Assemble().mat.Inverse(freedofs=fesHDiv.FreeDofs()) 
    GMRes(A=V.mat, pre=pre, b=rhs.vec, x=j.vec, tol=1e-8, maxsteps=5000, printrates=False)
Draw (j[0].real, mesh, draw_vol=False, order=3, min=-0.15, max=0.15);

Note: more demos are available here.