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

/*
to change locatation of daynamic lib. set   freefem++.pref   file like :
more freefem++.pref
loadpath = "../../plugin/seq/"
loadpath += "./"
*/

//load "medit"
load "MUMPS_mpi"
include "getARGV.idp"

verbosity=0;
real ttgv=1e10;



int nn=getARGV("-n",15);
int pplot=getARGV("-plot",0);


int[int] ll=[1,1,1,1,1,2];
mesh3 Th;
if(mpirank==0)
{
	Th=cube(nn,nn,nn,label=ll);
	Th=change(Th,fregion=nuTriangle%mpisize);
}
broadcast(processor(0),Th);
fespace VVh(Th,[P2,P2,P2,P1]);
fespace UUh(Th,[P2,P2,P2]);
fespace Uh(Th,P2);
fespace Ph(Th,P1);

macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM

func fup = (1-x)*(x)*y*(1-y)*16;


VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];

real timeI=mpiWtime();
real time1=mpiWtime();
varf vStokes([u1,u2,u3,p],[v1,v2,v3,q]) =
  int3d(Th,qforder=3,mpirank)( Grad(u1)'*Grad(v1) +  Grad(u2)'*Grad(v2) +  Grad(u3)'*Grad(v3)   //' for emacs
                  - div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
  + on(2,u1=fup,u2=0,u3=0) + on(1,u1=0,u2=0,u3=0) ;


matrix MStokes=vStokes(VVh,VVh,tgv=ttgv);
if(mpirank==0)
cout << "size of matrix " << MStokes.n << " x " << MStokes.m << " nn nzero coef  = " << MStokes.nbcoef << endl;
time1=mpiWtime()-time1;

real timeF=mpiWtime();
set(MStokes,solver=sparsesolver,master=-1);
timeF=mpiWtime()-timeF;

real time2=mpiWtime();
real[int] b=vStokes(0,VVh);
time2=mpiWtime()-time2;

real time3=mpiWtime();
u1[] = MStokes^-1*b;
time3=mpiWtime()-time3;

timeI=mpiWtime()-timeI;

if(mpirank==0)
{
cout << "============= CPU TIME ============" << endl;
cout << "size of matrix " << MStokes.n << " x " << MStokes.m << " nn nzero coef  = " << MStokes.nbcoef << endl;
cout << " matrix                  " <<  time1 << endl;
cout << " Fact                    " <<  timeF << endl;
cout << " second member           " <<  time2 << endl;
cout << " solve                   " <<  time3 << endl;
cout << "                          ------------" << endl;
cout << " all                     " <<  timeI << endl;
cout << "============= CPU TIME ============" << endl;
}
//if(mpirank==0 && pplot) medit("UV2 PV2",Th,[u1,u2,u3],p);
