//ff-mpirun -np 4 Helmholtz_circle_BEM_geneo.edp -ffddm_geneo_nu 20 -ffddm_geneo_threshold 100 -ns -wg

// for the make check:
// NBPROC 4
// PARAM -ffddm_geneo_nu 20 -ffddm_geneo_threshold 100

load "bem"

include "getARGV.idp"

macro dimension 3L//
include "ffddm.idp"

complex k = 10;

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

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

Th = OrientNormal(Th,unbounded=1);

ffddmbuildDmeshAug(H, Th, mpiCommWorld)

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

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

NewMacro norm()Hm1 EndMacro

macro VarfC(varfName, meshName, PhName)
IFMACRO(!norm)
  varf varfName(u,v) = int1d(meshName)(u*v);
ENDIFMACRO
IFMACRO(norm,Hm1)
  varf varfName(u,v) = int1d(meshName)(u*v + dx(u)*dx(v)+dy(u)*dy(v));
ENDIFMACRO
//

macro Varfrhs(varfName, meshName, PhName)
varf varfName(u,v) =
  -int1d(meshName)(finc*v);
//

ffddmsetupOperatorBem(H, H, Varf)

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

ffddmsetupPrecond(H, Varf)

real htoolEtabackup = htoolEta;

htoolEta = -1; // comment to compress E
ffddmextendedharmonicgeneoBEMsetup(H,Varf,VarfC)

htoolEta = htoolEtabackup;

/*
if (mpirank == 0) cout << HE.infos << endl;
display(HE);
*/

HVhi<complex> ui, bi;

ffddmbuildrhs(H, Varfrhs, bi[])

ui[] = HfGMRES(ui[], bi[], 1.e-6, 200, "right");

Hwritesummary

fespace Uh(Th,P1);
Uh<complex> u;

HfromVhi(ui[], Uh, u[])

// Mesh output
int np = 200;
int R = 4;

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));

fespace UhOut(ThOut,P1);

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

UhOut<complex> v, vinc;

vinc = finc;

v[] = B*u[];

v[] += vinc[];

UhOut vr = abs(v);

plot(vr, dim=2, fill=1, value=1, nbiso=20);
