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

// for the make check:
// NBPROC 2
// PARAM -frequency 2.e+8

load "bem"

include "getARGV.idp"

real c0= 299792458;
real f = getARGV("-frequency",6e+8);
real k = 2*pi*f/c0;
real lambda = c0/f;

cout << "lambda=" << lambda << endl;
// incident wave
real thetaI =  2*pi*0./360.;
real[int] dir=[cos(thetaI),0,sin(thetaI)];

// plane wave polarisation y
func fincx = 0;
func fincy = exp(1i*k*(x*dir[0]+y*dir[1]+z*dir[2]));
func fincz = 0;

int nloc = 4./lambda;
mesh3 Th3 = cube(nloc,nloc,nloc,[x-0.5,y-0.5,z-0.5]);

meshS ThS = extract(Th3);

plot(ThS);

macro def(u)[u,u#2,u#3]//

fespace Uh(ThS,RT0S);
Uh<complex> def(u), def(b);

varf vM([u1,u2,u3],[v1,v2,v3]) = int2dx2d(ThS)(ThS)(BEM(BemKernel("MA_SL",k=k),[u1,u2,u3],[v1,v2,v3]));  
HMatrix<complex> H = vM(Uh,Uh,eta=10,eps=1e-3);

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

varf vrhs(def(u),def(v)) = int2d(ThS)(def(v)'*[-fincx,-fincy,-fincz]);

b[] = vrhs(0,Uh);

// resolution of the courant J
u[] = H^-1*b[];

// computation of the potential to reconstruct E
nloc *= 20;
nloc -= nloc%8;
meshS ThOut = square3(nloc,nloc,[4*(x-0.5),4*(y-0.5),0]);
ThOut = trunc(ThOut, max(abs(x),abs(y)) > 0.5 + 4./nloc);

fespace UhOut(ThOut,[P1,P1,P1]);
 
varf vP([u1,u2,u3],[v1,v2,v3])=int2d(ThS)(POT(BemPotential("MA_SL",k=k),[u1,u2,u3],[v1,v2,v3])) ;
HMatrix<complex> B = vP(Uh,UhOut,eta=10,eps=1e-3);

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

UhOut<complex> [Ex,Ey,Ez];

Ex[] = B*u[];

fespace UhOutP1(ThOut,P1);
UhOutP1 Etotalr = sqrt(real(Ex+fincx)^2+ real(Ey+fincy)^2 + real(Ez+fincz)^2);

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