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 , we want to solve the Maxwell problem on the exterior domain :
where 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 in the space
such that,
where 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 :
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 FE−RT1 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 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 () and a quadrature method of order 3 when the two triangles are far ().
Note that the unknowns in RT approximation are the normal fluxes on the edge of the triangulation. In order to plot the potential , we have to move to a P1 representation, say . This can be done using a L2 projection from to :
This is what is done by the XLiFE++ function projection.
We obtain the following potential:
On the following figures, we show the approximated electric field and the exact electric field. The component is not shown because it is zero.