XLiFE++ website is now hosted here with up-to-date content. The current website, with potentially deprecated content, will be removed soon.

Example 15/18 - 3D Maxwell problem using EFIE

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 Γ be the boundary of a bounded domain Ω of 3, we want to solve the Maxwell problem on the exterior domain Ωe:

{ curlE ikH = 0 in Ωe curlH + ikE = 0 in Ωe E ×n = 0 on Γ lim |x| ((H Hinc) × x |x| (E Einc)) = 0(Silver-Muller condition)

where (Einc,Hinc) 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 J in the space

Hdiv(Γ) = {V L2(Γ)3,V.n = 0,div ΓV L2(Γ)}

such that, V Hdiv(Γ)

kΓΓJ(y)G(x,y).V(x) 1 kΓΓdivΓJ(y)G(x,y)divΓV(x) = ΓEinc.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 Ω.

Using the Stratton-Chu representation formula, the scattered electric field may be reconstructed in Ωe:

E(x) = Einc(x) + 1 kΓxG(x,y)divΓJ(y) + kΓG(x,y)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, triangle, 1, gmsh); 
  Domain Gamma = meshSh.domain("Gamma"); 
// define FERT1 space and unknown 
  Space V_h(_domain=Gamma, _interpolation=RT_1, _name="Vh"); 
  Unknown U(V_h, "U"); TestFunction V(U, "V"); 
// compute BEM system and solve it 
  IntegrationMethods ims(_SauterSchwabIM, 4, 0., _defaultRule, 5, 2., _defaultRule, 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, "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, "U3", 3); TestFunction V3(U3, "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, triangle, 1, 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, "W", 3); 
  IntegrationMethods im(_defaultRule, 10, 1., _defaultRule, 5); 
  TermVector Uext= 
   (1./k)*integralRepresentation(W, py0, intg(Gamma, grad_x(H)*div(U), im), J) 
      + k*integralRepresentation(W, py0, intg(Gamma, H*U, im), J); 
  saveToFile("Uext", Uext, vtu); 
// build exact solution, export to vtu file and compute error 
  TermVector Uexa(W, py0, Uex); 
  saveToFile("Uexa", Uexa, 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 Hdiv(Γ) 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 J, we have to move to a P1 representation, say J~. This can be done using a L2 projection from Hdiv(Γ) to L2(Γ):

ΓJ~|V =ΓJ|VV L2(Γ)

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

We obtain the following potential:


pict

Figure 11.1: 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 Ey is not shown because it is zero.


pict

Figure 11.2: 3D Maxwell problem on the unit sphere, using EFIE, x component


pict

Figure 11.3: 3D Maxwell problem on the unit sphere, using EFIE, y component