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

load "hpddm"                        // HPDDM plugin
macro dimension()2// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

macro def(i)i// EOM                 // scalar field definition
macro init(i)i// EOM                // scalar field initialization
macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
int[int] labPeriodic = [2, 4, 1, 3];
macro Pk() P2, periodic=[[labPeriodic[0],x+y], [labPeriodic[1],x+y], [labPeriodic[2],x-y], [labPeriodic[3],x-y]]// EOM

string deflation = getARGV("-deflation", "geneo");              // coarse space construction
int overlap = getARGV("-overlap", 1);                           // geometric overlap between subdomains
int fakeInterface = getARGV("-interface", 10);                  // interface between subdomains
int s = getARGV("-split", 1);                                   // refinement factor

mpiComm comm;
int p = getARGV("-hpddm_level_2_p", 1);
bool excluded = splitComm(mpiCommWorld, p, comm, topology = getARGV("-hpddm_level_2_topology", 0), exclude = (usedARGV("-hpddm_level_2_exclude") != -1));

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

real r = 0.25;
border a(t=0,1)    { x = -t+1; y = t;             label = 1; };
border b(t=0,1)    { x = -t; y = 1-t;             label = 2; };
border c(t=0,1)    { x = t-1; y = -t;             label = 3; };
border d(t=0,1)    { x = t; y = -1+t;             label = 4; };
border e(t=0,2*pi) { x = r*cos(t); y = -r*sin(t); label = 0; };

mesh ThBorder, Th = buildmesh(a(getARGV("-global", 40)) + b(getARGV("-global", 40)) + c(getARGV("-global", 40)) + d(getARGV("-global", 40)) + e(getARGV("-global", 40)));
fespace Wh(Th, Pk);           // local finite element space
int[int][int] intersection;   // local-to-neighbors renumbering
real[int] D;                  // partition of unity
{
    buildOverlapPeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, Pk, comm, excluded, labPeriodic)
}

func f = (y+x+1) * (y+x-1) * (y-x+1) * (y-x-1);
real[int] rhs;                              // local right-hand side
matrix<real> Mat;                           // local operator
{                                           // local weak form
    mesh ThAugmented = Th + ThBorder;
    varf vPb(u, v) = int2d(ThAugmented)(grad(u)' * grad(v)) - int2d(ThAugmented)((0.39 - f) * v) + on(0, u = 0.0);
    fespace WhAugmented(ThAugmented, Pk);
    Mat = vPb(WhAugmented, WhAugmented, tgv = -1);
    real[int] rhsFull = vPb(0, WhAugmented, tgv = -1);
    matrix R = interpolate(Wh, WhAugmented);
    renumbering(Mat, R, rhsFull, rhs);
}
ThBorder = square(1, 1);

schwarz A(Mat, intersection, D, scaled = false);
set(A, sparams = "-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction deflated -hpddm_geneo_nu 10");

matrix<real> Opt;                           // local operator with optimized boundary conditions
pair ret;
{
    int solver = getOption("schwarz_method");
    if(solver == 1 || solver == 2 || solver == 4) { // optimized Schwarz methods
        fespace Ph(Th, P0);
        real kZero = getARGV("-kZero", 10.0);
        Ph transmission = kZero;
        varf vOptimized(u, v) = int2d(Th)(grad(u)' * grad(v)) + int1d(Th, fakeInterface)(transmission * (u * v));
        Opt = vOptimized(Wh, Wh, tgv = -1);
    }
    if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) { // two-level Schwarz methods
        if(excluded)
            AttachCoarseOperator(mpiCommWorld, A/*, A = noPen, B = overlapRestriction, threshold = 2. * h[].max / diam*/);
        else {
            varf vPbNoPen(u, v) = int2d(Th)(grad(u)' * grad(v)) + on(1, u = 0.0);
            matrix<real> noPen = vPbNoPen(Wh, Wh, sym = 1);
            if(deflation == "geneo") // standard GenEO, no need for RHS -> deduced from LHS (Neumann matrix)
                AttachCoarseOperator(mpiCommWorld, A, A = noPen/*, threshold = 2. * h[].max / diam,*/, ret = ret);
            else if(deflation == "dtn") {
                varf vMass(def(u), def(v)) = int1d(Th, fakeInterface)(u * v);
                matrix<real> massMatrix = vMass(Wh, Wh, sym = 1);
                AttachCoarseOperator(mpiCommWorld, A, A = noPen, B = massMatrix, pattern = Opt/*, threshold = k,*/, ret = ret);
            }
            else if(deflation == "geneo-2") // GenEO-2 for optimized Schwarz methods, need for RHS (LHS is still Neumann matrix)
                AttachCoarseOperator(mpiCommWorld, A, A = noPen, B = Opt, pattern = Opt/*, threshold = 2. * h[].max / diam,*/, ret = ret);
        }
    }
}

Wh<real> def(u);    // local solution

if(Opt.n > 0)       // optimized Schwarz methods
    DDM(A, rhs, u[], excluded = excluded, ret = ret, O = Opt);
else
    u[] = A^-1 * rhs;

real[int] err(u[].n);
err = A * u[];      // global matrix-vector product
err -= rhs;

plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
u[] = err;
plotMPI(Th, u, Pk, def, real, cmm = "Global residual");
