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

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 def2(i)[i, i#B]// EOM
macro init2(i)[i, i]// EOM
func Pk = P1;                               // finite element space
func Zk = [P1, P2];                         // finite element space

mesh3 Th = cube(1, 1, 1);
fespace Wh(Th, Pk);            // local finite element space
fespace Qh(Th, Zk);            // local finite element space
int[int][int] intersectionP;   // local-to-neighbors renumbering
int[int][int] intersectionZ;   // local-to-neighbors renumbering
real[int] DP;                  // partition of unity
real[int] DZ;                  // partition of unity
{
    int[int] l = [1, 1, 1, 1];
    mesh3 ThBorder;
    Th = cube(getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5));    // global mesh
    fespace Ph(Th, P0);
    real[int] part(Ph.ndof);
    if(mpirank == 0)
        partitionerSeq(part, Th, mpisize);
    partitionerPar(part, Th, mpiCommWorld, mpisize);
    buildOverlapWithPartitioning(Th, ThBorder, part, 11111, 1, 1, intersectionZ, DP, Pk, mpiCommWorld, false);
    Th = cube(getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5)); // need to reload the global mesh because it is destroyed in ''build*''
    buildOverlapEdgeWithPartitioning(Th, ThBorder, part, 22222, 1, 1, intersectionP, DZ, Zk, mpiCommWorld, false, Zk, def2, init2);
    ThBorder = cube(1, 1, 1);
}
