//ff-mpirun -np 4 Helmholtz_circle_Dirichlet_all_indirect.edp -wg

// for the make check:
// NBPROC 4
// PARAM -k 2

load "bem"

include "getARGV.idp"

complex k = getARGV("-k",10.); // wavenumber

// incident wave
real[int] dir=[1,0];
func finc = exp(1i*k*(x*dir[0]+y*dir[1]));

//  Mesh
int n = 100;
border circle(t=0, 2*pi){x=cos(t); y=sin(t); z=0;}
meshL Th = buildmeshL(circle(n));

Th = OrientNormal(Th,unbounded=1);

// Dof
fespace Uh(Th,P1);

// Hmatrix for formulation of first kind
BemKernel ker1("SL",k=k);
varf vk1(u,v)=int1dx1d(Th)(Th)(BEM(ker1,u,v)) ;  
HMatrix<complex> HFirstKind = vk1(Uh,Uh,eta=10,eps=1e-3,maxleafsize=10);
if (mpirank == 0) cout << HFirstKind.infos << endl;
display(HFirstKind);

// Hmatrix for formulation of second kind
BemKernel ker2("DL",k=k);
varf vk2(u,v)=int1dx1d(Th)(Th)(BEM(ker2,u,v)) + int1d(Th)(0.5*u*v);  

HMatrix<complex> HSecondKind = vk2(Uh,Uh,eta=10,eps=1e-3,maxleafsize=10);
if (mpirank == 0) cout << HSecondKind.infos << endl;
display(HSecondKind);

// Hmatrix for combined formulation
BemKernel ker5=-1i*k*ker1+ker2;
varf vk3(u,v)=int1dx1d(Th)(Th)(BEM(ker5,u,v)) + int1d(Th)(0.5*u*v); 

HMatrix<complex> HCombined = vk3(Uh,Uh,eta=10,eps=1e-3,maxleafsize=10);

if (mpirank == 0) cout << HCombined.infos << endl;
display(HCombined);


// Solve formulation of first kind
Uh<complex> uFirstKind, bFirstKind;
varf vmassFirstKind(u,v) = int1d(Th)(-(finc)*v);
bFirstKind[] = vmassFirstKind(0,Uh);
uFirstKind[] = HFirstKind^-1*bFirstKind[];
plot(uFirstKind,dim=2,nbiso=20,value=1,cmm="First kind");

// Solve formulation of second kind 
Uh<complex> uSecondKind, bSecondKind;
bSecondKind[] = bFirstKind[];
uSecondKind[] = HSecondKind^-1*bSecondKind[];
plot(uSecondKind,dim=2,nbiso=20,value=1,cmm="Second kind");

// Solve combined formulation
Uh<complex> uCombined, bCombined;
bCombined[] = bFirstKind[];
uCombined[] = HCombined^-1*bCombined[];
plot(uCombined,dim=2,nbiso=20,value=1,cmm="Combined");


// Mesh output
int R = 4;
int np = 5*R*real(k);

border b1(t=-R, R){x=t; y=-R;}
border b2(t=-R, R){x=R; y=t;}
border b3(t=-R, R){x=-t; y=R;}
border b4(t=-R, R){x=-R; y=-t;}
border circleout(t=0, 2*pi){x=1.1*cos(t); y=1.1*sin(t); z=0;}

mesh ThOut = buildmesh(b1(np)+b2(np)+b3(np)+b4(np)+circleout(-np*pi/R));

// Incident wave in volume
fespace UhOut(ThOut,P1);
UhOut<complex> vinc;
vinc = finc;

// HMatrix for single layer potential
BemPotential SLPot("SL",k=k);
varf vpSL(u,v)=int1d(Th)(POT(SLPot,u,v)) ;  
HMatrix<complex> HSLPot = vpSL(Uh,UhOut,eta=10,eps=1e-3,maxleafsize=10);

if (mpirank == 0) cout << HSLPot.infos << endl;
display(HSLPot);

// HMatrix for double layer potential
BemPotential DLPot("DL",k=k);
varf vpDL(u,v)=int1d(Th)(POT(DLPot,u,v)) ;  
HMatrix<complex> HDLPot = vpDL(Uh,UhOut,eta=10,eps=1e-3,maxleafsize=10);

if (mpirank == 0) cout << HDLPot.infos << endl;
display(HDLPot);

// Solution for formulation of first kind
UhOut<complex> vFirstKind;
vFirstKind[] = HSLPot*uFirstKind[];
vFirstKind[] += vinc[];
UhOut vFirstKindabs = abs(vFirstKind);

plot(vFirstKindabs, dim=1, fill=1, value=1, nbiso=20,cmm="First kind");

// Solution for formulation of second kind
UhOut<complex> vSecondKind;
vSecondKind[] = HDLPot*uSecondKind[];
vSecondKind[] += vinc[];
UhOut vSecondKindabs = abs(vSecondKind);

plot(vSecondKindabs, dim=1, fill=1, value=1, nbiso=20,cmm="Second kind");

// Solution for combined formulation
UhOut<complex> vCombined1,vCombined2;
vCombined1[] = HSLPot*uCombined[];
vCombined2[] = HDLPot*uCombined[];
UhOut vCombinedabs = abs(-1i*k*vCombined1+vCombined2+vinc);

plot(vCombinedabs, dim=1, fill=1, value=1, nbiso=20,cmm="Combined");
