//ff-mpirun  -np 4 diffusion-3d-minimal-direct.edp -wg -raspart -global 20

// If you have openmpi you may need to add the option --oversubscribe to allow more processes than the number of cores available on your computer

// Abstract "Lap: Direct solver (MUMPS) :" << endl;//direct parallel solver 

// for the make check:
// NBPROC 4
// PARAM -raspart -global 20


macro dimension 3// EOM            // 2D or 3D

include "ffddm.idp"

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

int[int] LL = [2,2, 1,2, 2,2];
mesh3 ThGlobal = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10),
    [x, y, z], label = LL);      // global mesh

macro Varf(varfName, meshName, PhName)
    varf varfName(u,v) = int3d(meshName)(grad(u)' * grad(v)) + int3d(meshName)(v) + on(1, u = 1.0); // EOM
       
// Domain decomposition
ffddmbuildDmesh( Lap , ThGlobal , mpiCommWorld )
ffddmbuildDfespace( Lap , Lap , real , def , init , P1 )
ffddmsetupOperator(Lap ,Lap , Varf)

// Distributed Direct solve
real[int] rhs(LapVhi.ndof);//rhs(1) works as well 
ffddmbuildrhs(Lap , Varf , rhs )
LapVhi def(udirectsolver);

//Direct solve
udirectsolver[] = Lapdirectsolve(rhs);
Lapwritesummary//process 0 prints convergence history
ffddmplot(Lap,udirectsolver, "Lap Global solution with direct solver");


