//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4
// PARAM

load "bem"
load "PETSc-complex"
include "macro_ddm.idp"

if(!HasType("Mat", "Htool")) exit(0);

complex k     = getARGV("-k", 10.0);

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

border circle(t=0, 2*pi) { x=cos(t); y=sin(t); z=0; }
meshL ThG, ThL = buildmeshL(circle(n));
ThL = OrientNormal(ThL, unbounded = 1);
ThG = ThL;
int[int] n2o;
macro ThLN2O()n2o//
Mat<complex> A;
{
    macro dimension()3// EOM
    macro meshtype()L// EOM
    MatCreate(ThL, A, P1);
    Mat<complex> T;
    MatCreate(ThL, T, P0);
    CoherentGlobalMesh(T, A, ThL, ThG);
}

fespace Uh(ThL, P1);
fespace UhG(ThG, P1);
varf vk(u, v) = int1dx1d(ThG)(ThG)(BEM(BemKernel("SL", k=k), u, v));
A = vk(UhG, UhG);

varf vRHS(u, v) = -int1d(ThL)(finc*v);
Uh<complex> b, u;
b[] = vRHS(0, Uh);
Mat<complex> E, D;
MatConvert(A, E);
MatConvert(E, D);
set(D, sparams = "-ksp_view -ksp_monitor -pc_type cholesky");
ObjectView(D, format = "info");
complex[int] uP, bP;
ChangeNumbering(A, b[], bP);
uP.resize(bP.n);
KSPSolve(D, bP, uP);
ChangeNumbering(A, u[], uP, inverse = true, exchange = true);
// visualization
if(!NoGraphicWindow) {
    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; }

    mesh ThOutG, ThOutL = buildmesh(b1(np) + b2(np) + b3(np) + b4(np) + circle(-np*pi/R));
    fespace UhOut(ThOutL, P1);
    int[int] n2o;//
    macro ThOutLN2O()n2o//
    Mat<complex> B, T;
    ThOutG = ThOutL;
    MatCreate(ThOutL, B, P1);
    MatCreate(ThOutL, T, P0);
    CoherentGlobalMesh(T, B, ThOutL, ThOutG);
    MatDestroy(T);
    varf vp(u, v) = int1d(ThL)(POT(BemPotential("SL", k=k), u, v));
    Mat<complex> HPot(B, A);
    fespace UhOutG(ThOutG, P1);
    HPot = vp(UhG, UhOutG);
    UhOut<complex> v, vinc;
    vinc = finc;
    v[] = HPot*u[];
    v[] += vinc[];
    UhOut vr = abs(v);
    plotMPI(ThOutL, vr, P1, def, real, cmm = "Global first kind solution");
}
