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

Example 12/18 - Acoustic scattering by a sphere with single layer BEM

XLiFE++ is also able to deal with integral equation. This example illustrates the computation of the acoustic scattering by a sphere:

Using single layer potential leads to the integral equation, :

where $G$ is the Green function of the Helmhotz equation:

$G\left(x,y\right)=\frac{{e}^{ik|x-y|}}{4\pi |x-y|}.$

We deal with the variational formulation in $V={H}^{\frac{1}{2}}\left(S\right)$:

The solution $u$ is get from potential $p$ from the integral representation:

$u\left(x\right)={\int }_{S}G\left(x,y\right)p\left(y\right)dy.$

This example has been implemented in XLiFE++ using a ${P}^{0}$ Lagrange interpolation:

#include "xlife++.h"
using namespace xlifepp;

// incident plane wave
Complex uinc(const Point& p, Parameters& pa = defaultParameters)
{
Real kx=pa("kx"), ky=pa("ky"), kz=pa("kz");
Real kp=kx*p(1)+ky*p(2);
return exp(i_*kp);
}

int main(int argc, char** argv)
{
init(argc, argv, _lang=en); // mandatory initialization of xlifepp
// define parameters and functions
Parameters pars;
pars << Parameter(1., "k");  // wave number k
pars << Parameter(1., "kx") << Parameter(0., "ky") << Parameter(0., "kz"); // kx, ky, kz
pars << Parameter(1., "radius");  // disk radius
Kernel G = Helmholtz2dKernel(pars);  // load Helmholtz2D kernel
Function finc(uinc, pars);  // define right hand side function
Function scatSol(scatteredFieldDiskDirichlet, pars);  // exact solution

// meshing the unit disk
Number npa=16;  // nb of points by diameter of disk
Disk sp(_center=Point(0., 0.), _radius=1, _nnodes=npa, _domain_name="disk");
Mesh mS(sp, segment, 1 , gmsh);
Domain disk = mS.domain("disk");

// Lagrange P0 space and unknown
Space V1(_domain=disk, _interpolation=P1, _name="V1", _notOptimizeNumbering);
Unknown u1(V1, "u1");  TestFunction v1(u1, "v1");

// form definitions
IntegrationMethods ims(Duffy, 5, 0., Gauss_Legendre, 5, 1., Gauss_Legendre, 4, 2., Gauss_Legendre, 3 );
BilinearForm blf0=intg(disk, disk, u1*G*v1, ims);
LinearForm fv0 = -intg(disk, finc*v1);

// compute matrix and right hand side and solve system
TermMatrix A0(blf0, denseDualStorage, "A0");
TermVector B0(fv0, "B0");
TermVector U0 = directSolve(A0, B0);

// integral representation on x plane (far from disk), using P1 nodes
Number npp=20, npc=8*npp/10;
Real xm=4., eps=0.0001;
Point C1(0., -xm), C2(0., xm), C3(0., -xm);
SquareGeo sqx(_center=Point(0., 0.), _length=4., _nnodes=npp, _domain_name="Omega");
Disk dx(_center=Point(0., 0.), _radius=1.25, _nnodes=npc);
Mesh mx0(sqx-dx, triangle, 1 , gmsh);
Domain planx0 = mx0.domain("Omega");
Space Wx(_domain=planx0, _interpolation=P1, _name="Wx", _notOptimizeNumbering);
Unknown wx(Wx, "wx");
TermVector U0x0=integralRepresentation(wx, planx0, intg(disk, G*u1), U0);
TermMatrix Mx0(intg(planx0, wx*wx), "Mx0");

// compare to exact solution
TermVector solx0(wx, planx0, scatSol);
TermVector ec0x0=U0x0 - solx0;
theCout << "L2 error on x=0 plane: " << sqrt(abs((Mx0*ec0x0)|ec0x0)) << eol;

// export solution to file
saveToFile("U0", U0, vtu);
saveToFile("U0x0", U0x0, vtu);
return 0;
}