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

macro def(i)i// EOM                 // scalar field definition
macro init(i)i// EOM                // scalar field initialization
include "getARGV.idp"
int overlap = getARGV("-overlap", 1);
int s = getARGV("-split", 2);
int[int][int] intersection;
{
    macro dimension()2// EOM            // 2D or 3D
    include "macro_ddm.idp"             // additional DDM functions

    func Pk = P2;                       // finite element space
    int level = 4;
    mesh[int] ThTab(level);
    mesh[int] ThBorderTab(level);
    real[int][int] DTab(level);
    ThTab[level - 1] = square(getARGV("-global", 10), getARGV("-global", 10));
    matrix[int] P(level - 1);
    buildOverlapEdgeRecursive(ThTab, ThBorderTab, 10, s, overlap, level, P, intersection, DTab, Pk, mpiCommWorld, false, Pk, def, init);
    for(int i = 0; i < level; ++i) {
        fespace Vh(ThTab[i], Pk);
        Vh u;
        u[] = DTab[i];
        plotMPI(ThTab[i], u, Pk, def, real, cmm = "Global recursive partition of unity at level " + (i + 1));
    }
}
{
    macro dimension()3// EOM            // 2D or 3D
    include "macro_ddm.idp"             // additional DDM functions

    func Pk = P1;                       // finite element space
    int level = 3;
    mesh3[int] ThTab(level);
    mesh3[int] ThBorderTab(level);
    real[int][int] DTab(level);
    ThTab[level - 1] = cube(getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [x, y, z]);
    matrix[int] P(level - 1);
    buildOverlapEdgeRecursive(ThTab, ThBorderTab, -111111, s, overlap, level, P, intersection, DTab, Pk, mpiCommWorld, false, Pk, def, init);
    for(int i = 0; i < level; ++i) {
        fespace Vh(ThTab[i], Pk);
        Vh u;
        u[] = DTab[i];
        plotMPI(ThTab[i], u, Pk, def, real, cmm = "Global recursive partition of unity at level " + (i + 1));
    }
}
