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

load "PETSc"                        // PETSc plugin
include "macro_ddm.idp"             // additional DDM functions

macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
int[int] labPeriodic = [2, 4, 1, 3];

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 Th = buildmesh(a(getARGV("-global", 40)) + b(getARGV("-global", 40)) + c(getARGV("-global", 40)) + d(getARGV("-global", 40)) + e(getARGV("-global", 40)));
Mat A;
int[int] n2o;
macro ThN2O()n2o//
macro ThPeriodic()[[labPeriodic[0],x+y], [labPeriodic[1],x+y], [labPeriodic[2],x-y], [labPeriodic[3],x-y]]//
macro Pk()P2, periodic = ThPeriodic// EOM
DmeshCreate(Th);
MatCreate(Th, A, Pk);

fespace Wh(Th, Pk);                 // local finite element space
func f = (y+x+1) * (y+x-1) * (y-x+1) * (y-x-1);
varf vPb(u, v) = int2d(Th)(grad(u)' * grad(v)) - int2d(Th)((0.39 - f) * v) + on(0, u = 0.0);
real[int] rhs = vPb(0, Wh);

set(A, sparams = "-ksp_view");
Wh<real> u;                         // local solution

A = vPb(Wh, Wh);
u[] = A^-1 * rhs;

real[int] err = A * u[];            // global matrix-vector product
exchange(A, rhs, scaled = true);
err -= rhs;

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

set(A, sparams = "-pc_type gamg -ksp_type gmres -pc_gamg_threshold -1.0 -ksp_max_it 200");
u[] = 0.0;
u[] = A^-1 * rhs;
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
if(!NoGraphicWindow) {
    real[int] sol;
    ChangeNumbering(A, u[], sol);
    ChangeNumbering(A, u[], sol, inverse = true);
    mesh ThPlt = buildmesh(a(getARGV("-global", 40)) + b(getARGV("-global", 40)) + c(getARGV("-global", 40)) + d(getARGV("-global", 40)) + e(getARGV("-global", 40)));
    fespace WhPlt(ThPlt, Pk);
    WhPlt plt;
    WhPlt reduce;
    int[int] rest = restrict(Wh, WhPlt, n2o);
    plt[](rest) = u[];
    mpiReduce(plt[], reduce[], processor(0, mpiCommWorld), mpiSUM);
    if(mpirank == 0)
        plot(reduce);
}
