//ff-mpirun -np 4 Helmholtz-2d-interfaceBCs-comp.edp -ns -wg

// If you have openmpi you may need to add the option --oversubscribe to allow more processes than the number of cores available on your computer

// for the make check:
// NBPROC 4
// PARAM -frequency 2.0

macro dimension 2// EOM            // 2D or 3D

load "bem"
include "ffddm.idp"

macro def(i)i// EOM                         // scalar field definition
macro init(i)i// EOM                        // scalar field initialization
func Pk = P2;

int Dirichlet = 1;
int Robin = 2;

int[int] labs=[Robin, Robin, Robin, Robin];
//int[int] labs=[Dirichlet, Dirichlet, Dirichlet, Dirichlet];

real freq = getARGV("-frequency",4.0);
func k = 2.*pi*freq;
real lambda = (1./freq);
if (mpirank == 0) cout << "lambda = " << lambda << endl;

real nppwl = getARGV("-nppwl",10.0);
real h = lambda/nppwl;

real length = 1.;
real lengthpml = getARGV("-lengthpml",1.0)*lambda;

real nloc = nppwl*length/lambda;
int npl = nloc*(length+2*lengthpml)/length;
mesh Th = square(npl,npl,[x*(length+2*lengthpml)-lengthpml, y*(length+2*lengthpml)-lengthpml], label = labs);

func f = 100*exp(-((x-length/2)^2+(y-length/2)^2)*k*10);

macro Grad(u) [dx(u),dy(u)] // EOM

macro Varf(varfName, meshName, PhName)
    varf varfName(u,v) =
  int2d(meshName)(-k^2*u*v+Grad(u)'*Grad(v))
				- int1d(meshName,Robin)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

// for the preconditioner
macro VarfOptRobin(varfName, meshName, PhName)
    varf varfName(u,v) =
  int2d(meshName)(-k^2*u*v+Grad(u)'*Grad(v))
				- int1d(meshName,Robin,10)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

macro Varfrhs(varfName, meshName, PhName)
    varf varfName(u,v) =
                -int2d(meshName)(f*v)
                + on(Dirichlet,u=0); // EOM

ffddmbuildDmesh(Hmesh,Th,mpiCommWorld)

macro HFEdefmplot(u)real(u)//

ffddmbuildDfespace(HFE,Hmesh,complex,def,init,Pk)

vsym = 2; // symmetric, not hermitian
vtgv = 1.e+30;
vtgvelim = 1.e+30;

ffddmsetupOperator(H,HFE,Varf)

ffddmprecond = "oras";

complex[int] rhs(1);
ffddmbuildrhs(H,Varfrhs,rhs)

complex[int] x0(rhs.n);
x0 = 0;
HFEVhi<complex> def(u);

{
NewMacro HRobinmyA()
func complex[int] HRobinA(complex[int] &l) {return HA(l);}
EndMacro

ffddmsetupinit(HRobin,HFE)

ffddmsetupPrecond(HRobin,VarfOptRobin)

if (mpirank == 0) cout << "Robin:" << endl;

u[] = HRobinfGMRES(x0, rhs, 1.e-6, 200, "right");

HRobinwritesummary
}

{
NewMacro HBEMmyA()
func complex[int] HBEMA(complex[int] &l) {return HA(l);}
EndMacro

meshL ThiL = extract(HmeshThi);
ThiL = OrientNormal(ThiL,unbounded=1);
ThiL = trunc(ThiL,region==10);

fespace VhiL(ThiL,Pk);
fespace Chi = HFEVhi*VhiL;

varf vc(<[ufem],[ubem]>,<[v1],[v2]>) = int2d(HmeshThi)(-k^2*ufem*v1+Grad(ufem)'*Grad(v1)) // F
                              - int1d(HmeshThi,Robin)(1i*k*ufem*v1)
                              + on(Dirichlet,ufem=0)
                              + int1dx1d(ThiL)(ThiL)(BEM(BemKernel("TDL",k=k),ubem,v1)) + int1d(ThiL)(0.5*ubem*v1) // TDL
                              + int1d(ThiL)(ufem*v2)                        // mass
                              + int1dx1d(ThiL)(ThiL)(BEM(-1*BemKernel("SL",k=k),ubem,v2));  // -SL

matrix<complex> HC = vc(Chi,Chi,commworld=mpiCommSelf,factorize=3);

complex[int] rhsi(HFEVhi.ndof+VhiL.ndof);
rhsi = 0;

HFEVhi<complex> ufem;
VhiL<complex> ubem;

NewMacro HBEMmyPREC()
func complex[int] HBEMPREC(complex[int] &l) {
  rhsi(0:HFEVhi.ndof-1) = l;
  complex[int] res = HC^-1*rhsi;
  [ufem[],ubem[]] = res;

  HFEupdate(ufem[],true);

  return ufem[];
}
EndMacro

ffddmsetupinit(HBEM,HFE)

if (mpirank == 0) cout << "BEM:" << endl;

u[] = HBEMfGMRES(x0, rhs, 1.e-6, 200, "right");

HBEMwritesummary
}

ffddmplot(HFE,real(u), "Global solution");


