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

include "getARGV.idp"
include "cube.idp"
load "metis"
load "parmetis"
load "medit"
int[int] Nxyz = [getARGV("-global", 50), getARGV("-global", 50), getARGV("-global", 50)];
real[int,int] Bxyz = [[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]];
int[int,int] Lxyz = [[1, 2], [2, 2], [2, 2]];
mesh3 Th = Cube(Nxyz, Bxyz, Lxyz);
fespace Ph(Th, P0);
{
    mpiBarrier(mpiCommWorld);
    real time = mpiWtime();
    Ph part;
    if(mpirank == 0)
        metisdual(part[], Th, getARGV("-lpart", mpisize));
    broadcast(processor(0, mpiCommWorld), part[]);
    mpiBarrier(mpiCommWorld);
    if(mpirank == 0) {
        cout << "METIS: " << mpiWtime() - time << endl;
        medit("METIS", Th, part);
    }
}
{
    mpiBarrier(mpiCommWorld);
    real time = mpiWtime();
    Ph part;
    parmetis(part[], Th, getARGV("-lpart", mpisize), communicator = mpiCommWorld, worker = getARGV("-worker", mpisize));
    mpiBarrier(mpiCommWorld);
    if(mpirank == 0) {
        cout << "ParMETIS: " << mpiWtime() - time << endl;
        medit("ParMETIS", Th, part);
    }
}
