//ff-mpirun -np 4 Helmholtz-2d-FEM-BEM-coupling.edp -ffddm_partitioner 0 -ffddm_gmres_restart 1000 -wg

// for the make check:
// NBPROC 4
// PARAM -ffddm_partitioner 0 -ffddm_gmres_restart 1000

/* example of wave guiding with gradient-index lenses */

include "ffddm.idp"
load "bem"

complex k = 10; // wavenumber

real lambda = 2*pi/real(k);

real nloc = 10./lambda;
real n = nloc*2*pi;

int[int] nsl(10); // number of lenses
nsl = n;

real theta = 10*2*pi/360; // angular shift between lenses

real[int] nsx(nsl.n), nsy(nsl.n);
nsx[0] = 3.2;
nsy[0] = 0;

for (int i=1; i< nsl.n; i++) {
  nsx[i] = nsx[i-1] + 2.01*cos(i*theta);
  nsy[i] = nsy[i-1] + 2.01*sin(i*theta);
}

int interface = 1;
int waveguide = 2;

border circle(t=0, 2*pi; i){x=cos(t)+nsx[i]; y=sin(t)+nsy[i]; z=0; label=interface;}

real L = 2;
real l = 0.3;
real dd = 0.1;

func finc = exp(-100*((x+L-0.1)^2+y^2+z^2)); // source for waveguide excitation

// interface around the waveguide :
border a(t=-L-dd, L+dd){x=t; y=-l-dd; z=0; label=interface;}
border b(t=-l-dd, l+dd){x=L+dd; y=t; z=0; label=interface;}
border c(t=L+dd, -L-dd){x=t; y=l+dd; z=0; label=interface;}
border d(t=l+dd, -l-dd){x=-L-dd; y=t; z=0; label=interface;}

// waveguide :
border ga(t=-L, L){x=t; y=-l; z=0; label=waveguide;}
border gc(t=L, -L){x=t; y=l; z=0; label=waveguide;}
border gd(t=l, -l){x=-L; y=t; z=0; label=waveguide;}

mesh Thwaveguide = buildmesh(a(nloc*(2*L+2*dd))+b(nloc*(2*l+2*dd))+c(nloc*(2*L+2*dd))+d(nloc*(2*l+2*dd))
                         +ga(nloc*2*L)+gc(nloc*2*L)+gd(nloc*2*l));

mesh Thlenses = buildmesh(circle(nsl));
// we do not want any region with label 0 for the lenses, as 0 is the waveguide region in Thwaveguide :
Thlenses = change(Thlenses,fregion=region+1);

mesh Th = Thwaveguide + Thlenses; // FEM mesh

meshL ThL = extract(Th);
ThL = trunc(ThL,region==interface);
ThL = OrientNormal(ThL,unbounded=1); // BEM mesh

plot(ThL);

// change region labels in the fem mesh to define gradient index in lenses
fespace Ph(Th,P0);
Ph reg = region;

int[int] regs(2*nsl.n+2);
for (int i=0; i< nsl.n; i++){
  int regi = reg(nsx[i], nsy[i]); // the lenses
  regs[2*i] = regi;
  regs[2*i+1] = i+1; 
}
regs[2*nsl.n] = reg(0,0); // the waveguide
regs[2*nsl.n+1] = -1;
Th = change(Th,region=regs);

func ind = region == -1 ? 1 : 2./(1+((x-nsx[region-1])^2+(y-nsy[region-1])^2)); // gradient index in lenses

// DD mesh partitioning
fespace FPh(Th,P0);
FPh fpart;
if (mpirank == 0) {
  int[int] lpart(FPh.ndof);
  if (mpisize > 1)
    scotch(lpart,Th,mpisize);
  else
    lpart = 0;
  for(int i=0;i<lpart.n;++i)
    fpart[][i]=lpart[i];
}
broadcast(processor(0), fpart[]);

macro Fsimple(PhGlobal, part, comm) part[] = fpart[]; // EOM

macro Bsimple(PhGlobal, part, comm) part = fpart; // EOM

macro Fdim 2//
ffddmbuildDmesh(F, Th, mpiCommWorld)

macro Bdim 3L//
ffddmbuildDmesh(B, ThL, mpiCommWorld)

macro def(i)i// EOM                         // scalar field definition
macro init(i)i// EOM                        // scalar field initialization
ffddmbuildDfespace(F, F, complex, def, init, P1)

ffddmbuildDfespace(B, B, complex, def, init, P1)

macro VarfSL(varfName, meshName, PhName)
  varf varfName(u,v) = int1dx1d(meshName)(meshName)(BEM(BemKernel("SL",k=k),u,v));
//

macro VarfTDL(varfName, meshName, PhName)
varf varfName(u,v) = int1dx1d(meshName)(meshName)(BEM(BemKernel("TDL",k=k),u,v)) + int1d(meshName)(0.5*u*v);  

// assemble BEM blocks :
ffddmprecond = "ras";
ffddmsetupOperatorBem(B, B, VarfSL)

ffddmsetupPrecond(B, VarfSL)

ffddmsetupOperatorBem(TDL, B, VarfTDL)

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

BVhi<complex> ubi, bbi;

bbi[] = 0;

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

macro VarfF(varfName, meshName, PhName)
    varf varfName(u,v) =
  int2d(meshName)(-ind*k^2*u*v+Grad(u)'*Grad(v))
  + on(waveguide,u=0); // EOM

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

macro VarfFrhs(varfName, meshName, PhName)
    varf varfName(u,v) =
                int2d(meshName)(finc*v)
                + on(waveguide,u=0); // EOM

macro VarfFmass(varfName, meshName, PhName)
    varf varfName(u,v) =
  int1d(meshName,interface)(u*v); // EOM

// assemble FEM blocks :
ffddmprecond = "oras";
ffddmsetup(F, F, VarfF, VarfFOpt)

ffddmsetupOperator(Fmass, F, VarfFmass)

FVhi<complex> ufi, bfi;

ffddmbuildrhs(F, VarfFrhs, bfi[])

int Cn = FVhi.ndof + BVhi.ndof;

// rhs for the coupled problem :
complex[int] bC = [bfi[], bbi[]];

// restriction matrix from local fem space to local bem space :
matrix<complex> RFtoB;
{
  matrix tmp = interpolate(BVhi,FVhi);
  RFtoB = tmp;
}

/* the coupled problem we want to solve is :
(  F     TDL ) (ufem) = (bfi)
(  mass  -SL ) (ubem)   ( 0 ) */

// define custom operators for the coupled problem :
macro Cmyscalprod
func complex Cscalprod(complex[int] va, complex[int] vb)
{
  complex resf = Fscalprod(va(0:FVhi.ndof-1),vb(0:FVhi.ndof-1));
  complex resb = Bscalprod(va(FVhi.ndof:Cn-1),vb(FVhi.ndof:Cn-1));
  return resf+resb;
}
//

ffddminitDfespacef(C, B, complex, def, init, P1, def, init, P1)

macro CmyA
func complex[int] CA(complex[int] &u) 
{
  complex[int] res(u.n);
  res(0:FVhi.ndof-1) = FA(u(0:FVhi.ndof-1));
  res(FVhi.ndof:Cn-1) = (-1)*BA(u(FVhi.ndof:Cn-1));
  
  complex[int] uFm = FmassA(u(0:FVhi.ndof-1));
  complex[int] tmp = RFtoB*uFm;
  Bupdate(tmp,true);
  res(FVhi.ndof:Cn-1) += tmp;
  
  complex[int] utdl = TDLA(u(FVhi.ndof:Cn-1));
  complex[int] tmp2 = RFtoB'*utdl;
  Fupdate(tmp2,true);
  res(0:FVhi.ndof-1) += tmp2;

  return res;
}
//

// block diagonal preconditioner :
macro CmyPREC1
func complex[int] CPREC1(complex[int] &u) 
{
  complex[int] res(u.n);
  res(0:FVhi.ndof-1) = FPREC1(u(0:FVhi.ndof-1));
  res(FVhi.ndof:Cn-1) = (-1)*BPREC1(u(FVhi.ndof:Cn-1));
  return res;
}
//

macro Cwithhpddmkrylov()1//

ffddmsetupinit(C,C);

complex[int] uC(Cn);

uC = CfGMRES(uC, bC, 1.e-4, 2000, "right");

Fwritesummary

Bwritesummary

Cwritesummary

ufi[] = uC(0:FVhi.ndof-1); // FEM local solution

ubi[] = uC(FVhi.ndof:Cn-1); // BEM local solution

FVhglob<complex> ufem;
FfromVhi(ufi[],FVhglob,ufem[]) // FEM global solution

ffddmplot(F,real(ufi),"FEM solution");

BVhglob<complex> ubem;

BfromVhi(ubi[], BVhglob, ubem[]) // BEM global solution

// output mesh for visualization
int np = 200/2;
real R = 20;
real rr = 20;

border b1(t=-rr, R){x=t; y=-rr;}
border b2(t=-rr, rr){x=R; y=t;}
border b3(t=R, -rr){x=t; y=rr;}
border b4(t=rr, -rr){x=-rr; y=t;}

nsl = -nsl;
// exterior mesh
mesh ThOut = buildmesh(b1(np*R/rr)+b2(np)+b3(np*R/rr)+b4(np)+circle(nsl)
              +a(-nloc*(2*L+2*dd))+b(-nloc*(2*l+2*dd))+c(-nloc*(2*L+2*dd))+d(-nloc*(2*l+2*dd)));
broadcast(processor(0),ThOut);

fespace UhOut(ThOut,P1);

varf vp(u,v)=int1d(ThL)(POT(BemPotential("SL",k=k),u,v));
HMatrix<complex> B = vp(BVhglob,UhOut);
if (mpirank == 0) cout << B.infos << endl;
//display(B);

UhOut<complex> uext;

uext[] = B*ubem[];

plot(ufem, uext, dim=2, fill=1, value=1, nbiso=40);
