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

load "PETSc"                        // PETSc plugin
macro removeZeros()true// EOM
macro trueRestrict()true// EOM
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
func Pk = P1;                       // finite element space

int s = getARGV("-split", 1);                                   // refinement factor

mesh Th = square(1, 1);
fespace Wh(Th, Pk);           // local finite element space
int[int][int] intersection;   // local-to-neighbors renumbering
real[int] D;                  // partition of unity
{
    Th = square(getARGV("-global", 40), getARGV("-global", 40));    // global mesh
    build(Th, s, intersection, D, Pk, mpiCommWorld)
}

varf vPb(u, v) = int2d(Th)(grad(u)' * grad(v)) + int2d(Th)(v) + on(1, u = 0.0);

matrix Loc = vPb(Wh, Wh, tgv = -2);
Mat A(Loc, intersection, D);
Wh<real> def(u);    // local solution

set(A, sparams = "-pc_type lu -ksp_type gmres");
real[int] rhs = vPb(0, Wh, tgv = -2);
real[int] rhsPETSc;
ChangeNumbering(A, rhs, rhsPETSc);
func real[int] funcRes(real[int]& inPETSc) {
    real[int] outPETSc(inPETSc.n);
    MatMult(A, inPETSc, outPETSc);
    outPETSc -= rhsPETSc;
    return outPETSc;
}
func int funcJ(real[int]& inPETSc) {
    return 0;
}
Wh xu;
xu = sqrt((x - 0.5)^2 + (y - 0.5)^2) < 0.4 ? 0.1 : 1.0;
plotMPI(Th, xu, Pk, def, real, cmm = "Upper bound");
real[int] xPETSc, xuPETSc;
u[] = 0;
ChangeNumbering(A, u[], xPETSc);
ChangeNumbering(A, xu[], xuPETSc);
SNESSolve(A, funcJ, funcRes, xPETSc, sparams = "-snes_monitor -snes_view -snes_vi_monitor -snes_type vinewtonrsls", xu = xuPETSc);
ChangeNumbering(A, u[], xPETSc, inverse = true, exchange = true);
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
