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

load "PETSc"                        // PETSc plugin
load "mat_edgeP1"
load "Element_Mixte3d"
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"
include "cube.idp"

macro Curl(ux, uy, uz)[dy(uz)-dz(uy), dz(ux)-dx(uz), dx(uy)-dy(ux)]// EOM
macro CrossN(ux, uy, uz)[uy*N.z-uz*N.y, uz*N.x-ux*N.z, ux*N.y-uy*N.x]// EOM
func Pk = Edge03d;
func PkNode = P1;

int Dirichlet = 1;

mesh3 Th = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10));
int[int] chlab = [1, Dirichlet, 2, Dirichlet, 3, Dirichlet, 4, Dirichlet, 5, Dirichlet, 6, Dirichlet];
Th = change(Th, refface = chlab);

fespace Wh(Th, Pk);

Mat A;
int[int] n2o;
macro ThN2O()n2o// EOM
DmeshCreate(Th);
{
    macro ThPkPart()Edge03ds0// EOM special FE for the domain partitioning
    MatCreate(Th, A, Pk);
}

real[int] rhs(Wh.ndof);
varf vPb([Ex,Ey,Ez],[vx,vy,vz]) =
    int3d(Th)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
    + int3d(Th)([vx,vy,vz]'*[Ex,Ey,Ez])
    + on(Dirichlet, Ex=0, Ey=0, Ez=0);
A = vPb(Wh, Wh);
func f = exp(-2.0*((x-0.5)^2+(y-0.5)^2+(z-0.5)^2));
varf vPbRhs([Ex,Ey,Ez],[vx,vy,vz]) =
    - int3d(Th)([vx,vy,vz]'*[0,0,f])
    + on(Dirichlet, Ex=0,Ey=0,Ez=0);
rhs = vPbRhs(0, Wh);

{
    Mat B;
    MatCreate(Th, B, P1);
    fespace VhNode(Th, PkNode);
    matrix GLoc;
    MatrixEdgeP1(GLoc, Th);
    Mat G(A, B, GLoc);
    VhNode coord = x;
    real[int] tmp;
    ChangeNumbering(B, coord[], tmp);
    real[int, int] coordPETSc(3, tmp.n);
    for(int j = 0; j < 3; ++j) {
        for[i, ui : tmp] coordPETSc(j, i) = ui;
        if(j == 0)
            coord = y;
        else if(j == 1)
            coord = z;
        ChangeNumbering(B, coord[], tmp);
    }
    set(A, sparams = "-pc_type hypre -pc_hypre_type ams -ksp_monitor -ksp_view -ksp_max_it 40", coordinates = coordPETSc, gradient = G);
}

Wh [solx, soly, solz];
solx[] = A^-1 * rhs;
if(!NoGraphicWindow) {
    real[int] sol;
    ChangeNumbering(A, solx[], sol);
    ChangeNumbering(A, solx[], sol, inverse = true);
    mesh3 ThPlt = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10));
    fespace WhPlt(ThPlt, Pk);
    WhPlt [pltx, plty, pltz];
    WhPlt [reducex, reducey, reducez];
    int[int] rest = restrict(Wh, WhPlt, n2o);
    for[i, v : rest] pltx[][v] = solx[][i];
    mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
    if(mpirank == 0)
        medit("Global solution", ThPlt, [real(reducex), real(reducey), real(reducez)]);
}
