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:

{ Δu + k2u = 0 in Ω = ℝ3∕B(0,R) u = −uinc  on S

Using single layer potential leads to the integral equation, :

∫ SG(x,y)p(x)dx = −uinc on S

where G is the Green function of the Helmhotz equation:

G(x,y) = eik|x−y| 4π|x − y|.

We deal with the variational formulation in V = H1 2 (S):

| find p ∈ V  such that  ∫ S ∫ Sp(x)G(x,y)q¯(y)dxdy = −∫ Suincq¯∀ ⁡q ∈ V.

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

u(x) = ∫ SG(x,y)p(y)dy.

This example has been implemented in XLiFE++ using a P0 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 
  numberOfThreads(2); 
// 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; 
}

pict

Figure 8.1: Solution of the 3D Helmholtz problem using single layer BEM on the unit sphere