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

load "hpddm"                        // HPDDM plugin
macro partitioner()metis// EOM      // metis, scotch, or parmetis
macro dimension()3// EOM            // 2D or 3D
load "Element_Mixte3d"
include "macro_ddm.idp"
include "cube.idp"

macro def(i)[i, i#y, i#z]// EOM             // vector field definition
macro init(i)[i, i, i]// EOM                // vector field initialization
macro defPart(u)u// EOM                     // partition of unity definition
macro initPart(u)u// EOM                    // partition of unity initialization
macro Curl(ux, uy, uz)[dy(uz)-dz(uy), dz(ux)-dx(uz), dx(uy)-dy(ux)]// EOM
macro CrossN(ux, uy, uz)[uy*N.z-uz*N.y, uz*N.x-ux*N.z, ux*N.y-uy*N.x]// EOM
func Pk = Edge03d;
func PkPart = Edge03ds0;

int Dirichlet = 1, Robin = 2;

int overlap = getARGV("-overlap", 2);                           // geometric overlap between subdomains
int fakeInterface = getARGV("-interface", 10);                  // interface between subdomains
int s = getARGV("-split", 1);                                   // refinement factor

if(verbosity > 0 && mpirank == 0) {
    cout << " --- " << mpirank << "/" << mpisize;
    cout << " - maxwell-3d.edp - input parameters: refinement factor = " << s << " - overlap = " << overlap << endl;
}

int [int,int] LL = [[1,2], [3,4], [5,6]];
real[int, int] BB = [[0,1], [0,1], [0,1]];
int[int] NN = [getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10)];
int[int] CC = [1, 1, 1];
mesh3 Th = Cube(NN, BB, LL);
func k = 6 * pi;
int[int] chlab = [1, Robin, 2, Robin, 3, Robin, 4, Robin, 5, Robin, 6, Robin];
Th = change(Th, refface = chlab);

fespace Wh(Th, Pk);

mesh3 ThBorder;
int[int][int] intersection;   // local-to-neighbors renumbering
real[int] D;
{
    buildOverlapEdge(Th, ThBorder, 10, s, overlap, intersection, D, Pk, mpiCommWorld, 0, PkPart, defPart, initPart)
}

matrix<complex> Mat;
matrix<complex> Opt;
complex[int] rhs(Wh.ndof);
{
    mesh3 ThAugmented = Th + ThBorder;
    ThBorder = Cube(CC, BB, LL);
    fespace WhAugmented(ThAugmented, Pk);
    varf vPb([Ex,Ey,Ez],[vx,vy,vz]) =
        int3d(ThAugmented)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
        - int3d(ThAugmented)(k^2*[vx,vy,vz]'*[Ex,Ey,Ez])
        - int2d(ThAugmented, Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
        + on(Dirichlet, Ex=0, Ey=0, Ez=0);
    Mat = vPb(WhAugmented, WhAugmented, solver=GMRES);
    func f = exp(-((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)*60);
    varf vPbRhs([Ex,Ey,Ez],[vx,vy,vz]) =
        - int3d(ThAugmented)([vx,vy,vz]'*[0,0,f])
        + on(Dirichlet, Ex=0,Ey=0,Ez=0);
    if(mpisize > 1) {
        complex[int] rhsFull(WhAugmented.ndof);
        rhsFull = vPbRhs(0, WhAugmented);
        matrix R = interpolate(Wh, WhAugmented);
        R.thresholding(1.0e-10);
        renumbering(Mat, R, rhsFull, rhs);
    }
    else
        rhs = vPbRhs(0, WhAugmented);
}

schwarz<complex> A(Mat, intersection, D, scaled = false, communicator = mpiCommWorld);
set(A, sparams = "-hpddm_schwarz_method oras");

complex[int] nEx(Wh.ndof);
nEx = 0.0;

int solver = getOption("schwarz_method");
if(solver == 1 || solver == 2 || solver == 4) { // optimized Schwarz methods
    varf vPbOpt([Ex,Ey,Ez],[vx,vy,vz]) =
        int3d(Th)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
        - int3d(Th)(k^2*[vx,vy,vz]'*[Ex,Ey,Ez])
        - int2d(Th, Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
        - int2d(Th, fakeInterface)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
        + on(Dirichlet, Ex=0, Ey=0, Ez=0);
    Opt = vPbOpt(Wh, Wh);
}
DDM(A, rhs, nEx, O = Opt);
Wh<complex> def(sol);
sol[] = nEx;
medit("Global solution", Th, [real(sol), real(soly), real(solz)]);
