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

load "hpddm"                        // HPDDM plugin
macro dimension 3// 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), dz(u)]// EOM    // three-dimensional gradient
func Pk = P1;                               // finite element space

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

int[int] LL = [2,2, 1,2, 2,2];
mesh3 Th = cube(1, 1, 1, [x, y, z]);
fespace Wh(Th, Pk);           // local finite element space
int[int][int] intersection;   // local-to-neighbors renumbering
real[int] D;                  // partition of unity
{
    Th = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10), [x, y, z], label = LL);      // global mesh
    build(Th, s, intersection, D, Pk, mpiCommWorld)
}

varf vPb(u, v) = int3d(Th)(grad(u)' * grad(v)) + int3d(Th)(v) + on(1, u = 1.0);
matrix<real> Mat = vPb(Wh, Wh, tgv = -1);
real[int] rhs = vPb(0, Wh, tgv = -1);

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

int solver = getOption("schwarz_method");
if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) { // two-level Schwarz methods
    matrix<real> noPen = vPb(Wh, Wh, sym = 1);
    AttachCoarseOperator(mpiCommWorld, A, A = noPen);
}

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

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");
