3D Maxwell problem using EFIE#

Keywords: PDE:Maxwell Method:BEM FE:Lagrange FE:RaviartThomas Solver:Direct

Solving diffraction of an electromagnetic plane wave on a obstacle using BEM is more intricate. Indeed, it is a vector problem and it involves Raviart-Thomas elements. We show how XLiFE++ can deal easily with.

Let \(\Gamma\) be the boundary of a bounded domain \(\Omega\) of \(\mathbb{R}^3\), we want to solve the Maxwell problem on the exterior domain \(\Omega_e\):

\[\begin{split}\DeclareMathOperator{\curl}{curl} \begin{cases} \curl \mathbf{E} - i k \mathbf{H} = 0 & \mathrm{in\;} \Omega_e \\ \curl \mathbf{H} + i k \mathbf{E} = 0 & \mathrm{in\;} \Omega_e \\ \mathbf{E} \times \mathbf{n} = 0 & \mathrm{on\;}\Gamma\\ \displaystyle \lim_{|x| \to \infty} \left( \left(\mathbf{H} - \mathbf{H}_{inc} \right) \times \dfrac{x}{|x|} - \left( \mathbf{E}-\mathbf{E}_{inc}\right)\right) = 0 & \mathrm{(Silver-Muller\;condition)} \end{cases}\end{split}\]

where \((\mathbf{E}_{inc},\mathbf{H}_{inc})\) is an incident field (a solution of Maxwell equation in free space), for instance a plane wave.

The EFIE (Electric Field Integral Equation) consists in finding the potential \(\mathbf{J}\) in the space

\[\DeclareMathOperator{\divop}{div} H_{\divop}(\Gamma) = \left\{\mathbf{V} \in L^2(\Gamma)^3, \mathbf{V} \cdot \mathbf{n}=0, \divop_{\Gamma} \mathbf{V} \in L^2(\Gamma) \right\}\]

such that:

\[\DeclareMathOperator{\divop}{div} \forall \mathbf{V}\in H_{\divop}(\Gamma), k \int_{\Gamma} \int_{\Gamma} \mathbf{J}(y) * G(x,y) * \mathbf{V}(x) - \dfrac{1}{k}\int_{\Gamma} \int_{\Gamma} \div_{\Gamma} \mathbf{J}(y) * G(x,y) * \div_{\Gamma} \mathbf{V}(x) = -\int_{\Gamma} \mathbf{E}_{inc} \cdot \mathbf{V}\]

where \(G\) is the Green function related to the Helmholtz 3D problem in free space.

This equation has a unique solution, except for a discrete set of wavenumbers corresponding to the resonance frequencies of the cavity \(\Omega\).

Using the Stratton-Chu representation formula, the scattered electric field may be reconstructed in \(\Omega_e\):

\[\DeclareMathOperator{\divop}{div} \mathbf{E}(x) = \mathbf{E}_{inc}(x) + \dfrac{1}{k}\int_\Gamma \nabla_x G(x,y) \divop_{\Gamma} \mathbf{J}(y) + k \int_\Gamma G(x,y) * \mathbf{J}(y).\]

This problem is implemented in XLiFE++ as follows:

#include "xlife++.h"
using namespace xlifepp;
Vector<complex_t> data_incField(const Point& P, Parameters& pars)
{
  Vector<real_t> incPol(3, 0.); incPol(1)=1.; Point incDir(0., 0., 1.) ;
  Real k = pars("k");
  return  incPol * exp(i_*k * dot(P, incDir));
}

Vector<complex_t> uinc(const Point& P, Parameters& pars)
{
  Vector<real_t> incPol(3, 0.); incPol(1)=1.; Point incDir(0., 0., 1.) ;
  Real k = pars("k");
  return  incPol*exp(i_*k * dot(P, incDir));
}

int main(int argc, char** argv)
{
  init(argc, argv, _lang=en);
  // define parameters and functions
  Real k= 1, R=1.; Parameters pars;
  pars << Parameter(k, "k") << Parameter(R, "radius");
  Kernel H = Helmholtz3dKernel(pars);          // load Helmholtz3D kernel
  Function Einc(data_incField, pars);          // define right hand side 
  Function Uex(scatteredFieldMaxwellExn, pars);
  // meshing the unit sphere
  Number npa=15; Point O(0, 0, 0);                          
  Sphere sphere(_center=O, _radius=R, _nnodes=npa, _domain_name="Gamma");
  Mesh meshSh(sphere, _shape=triangle, _generator=gmsh);
  Domain Gamma = meshSh.domain("Gamma");
  // define FE-RT1 space and unknown
  Space V_h(_domain=Gamma, _interpolation=RT_1, _name="Vh");
  Unknown U(V_h, _name="U"); TestFunction V(U, _name="V");
  // compute BEM system and solve it
  IntegrationMethods ims(_method=SauterSchwabIM, _order=4, _bound=0.,
                         _quad=defaultQuadrature, _order=5, _bound=2.,
                         _quad=defaultQuadrature, _order=3);
  BilinearForm auv = k*intg(Gamma, Gamma, U*H|V, ims)
                   -(1./k)*intg(Gamma, Gamma, div(U)*H*div(V), ims);
  TermMatrix A(auv, _name="A");
  TermVector B(-intg(Gamma, Einc|V));
  TermVector J = directSolve(A, B);
  // get P1 representation of solution and export it to vtu file
  Space L_h(_domain=Gamma, _interpolation=P1, _name="Lh");
  Unknown U3(L_h, _name="U3", _dim=3); TestFunction V3(U3, _name="V3");
  TermVector JP1=projection(J, L_h, 3, _L2Projector);
  saveToFile("JP1", JP1(U3[1]), vtu);
  // integral representation on y=0 plane (excluding sphere), using P1 nodes
  Number npp=30, npc=5;
  Square sqx(_center=O, _length=20., _nnodes=npp, _domain_name="Omega");
  Disk dx(_center=O, _radius=1.2*R, _nnodes=npc);
  Mesh mx0(sqx-dx, _shape=triangle, _generator=gmsh);
  mx0.rotate3d(1., 0., 0., pi_/2);
  Domain py0 = mx0.domain("Omega");
  Space Vy0(_domain=py0, _interpolation=P1, _name="Vy0", _notOptimizeNumbering);
  Unknown W(Vy0, _name="W", _dim=3);
  IntegrationMethods im(_quad=defaultQuadrature, _order=10, _bound=1., _quad=defaultQuadrature, _order=5);
  TermVector Uext=(1./k)*integralRepresentation(W, py0, intg(Gamma, grad_x(H)*div(U), _method=im), J)
                  + k*integralRepresentation(W, py0, intg(Gamma, H*U, _method=im), J);
  saveToFile("Uext", Uext, _format=vtu);
  // build exact solution, export to vtu file and compute error
  TermVector Uexa(W, py0, Uex);
  saveToFile("Uexa", Uexa, _format=vtu);
  TermMatrix M(intg(py0, W|W));
  TermVector E=Uext-Uexa;
  theCout << "L2 error = " << sqrt(real(M*E|E)) << eol;
  return 0;
}

In order to build an approximated space of \(\DeclareMathOperator{\divop}{div} H_{\divop}(\Gamma)\) we use the Raviart-Thomas element of order 1.

As the integrals involved in bilinear form are singular, we use here the Sauter-Schwab method to compute them when two triangles are adjacent, a quadrature method of order 5 if the two triangles are close (\(0<d(T1,T2)<2h\)) and a quadrature method of order 3 when the two triangles are far (\(d(T1,T2)>=2h\)).

Note that the unknowns in RT approximation are the normal fluxes on the edge of the triangulation. In order to plot the potential \(\mathbf{J}\), we have to move to a P1 representation, say \(\widetilde{\mathbf{J}}\). This can be done using a L2 projection from \(\DeclareMathOperator{\divop}{div} H_{\divop}(\Gamma)\) to \(L^2(\Gamma)\):

\[\forall V\in L^2(\Gamma), \int_{\Gamma} \widetilde{\mathbf{J}} \cdot \mathbf{V} = \int_{\Gamma} \mathbf{J} \cdot \mathbf{V}\]

This is what is done by the XLiFE++ function projection().

We obtain the following potential:

../_images/J_real.png — Fig. 28 3D Maxwell problem on the unit sphere, using EFIE, potential#

On the following figures, we show the approximated electric field and the exact electric field. The component \(E_y\) is not shown because it is zero.

../_images/Ex_real.png — Fig. 29 3D Maxwell problem on the unit sphere, using EFIE, x component#

../_images/Ez_real.png — Fig. 30 3D Maxwell problem on the unit sphere, using EFIE, y component#