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

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

real ttgv=1e10;
//string ssparams="nprow=1, npcol="+mpisize;

int nn=getARGV("-n",15);
mesh Th2=square(nn,nn);
fespace Vh2(Th2,P2);  Vh2 ux,uz,p2;
int[int] rup=[0,2],  rdown=[0,1], rmid=[1,1,2,1,3,1,4,1];
real zmin=0,zmax=1;
mesh3 Th=buildlayers(Th2,nn,
  zbound=[zmin,zmax],  labelmid=rmid, 
  reffaceup = rup,     reffacelow = rdown);
  
//medit("c10x10x10",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;

Uh u1,u2,u3;
Uh v1,v2,v3;
Ph p,q;

real timeI=clock();
real time1=clock();

//VVh [uu1,uu2,uu3,up];
varf vlaplace([u1,u2,u3],[v1,v2,v3]) = 
  int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) +  Grad(u2)'*Grad(v2) +  Grad(u3)'*Grad(v3) ) 
  + on(2,u1=fup,u2=0,u3=0) + on(1,u1=0,u2=0,u3=0) ; //'

varf vdiv([u1,u2,u3],[q])=int3d(Th,qforder=3)(- div(u1,u2,u3)*q );
varf vdivT([q],[u1,u2,u3])=int3d(Th,qforder=3)(- div(u1,u2,u3)*q );
varf vmass(p,q) = int3d(Th,qforder=3)( 1e-10*p*q );

matrix M=vmass(Ph,Ph);
matrix L=vlaplace(UUh,UUh,tgv=ttgv);
matrix B=vdiv(UUh,Ph);
matrix BT=vdivT(Ph,UUh);

matrix A = [ [ L, BT ],     //'
	     [ B, M  ]];

time1=clock()-time1;

real timeF=clock();
set(A,solver=sparsesolver);
timeF=clock()-timeF;

real time2=clock();
real[int] b(VVh.ndof);
real[int] xx(VVh.ndof);
{ 
  real[int] b1 = vlaplace(0,UUh);
  
  for(int ii=0; ii<UUh.ndof; ii++)
    {
      b[ii]=b1[ii];
    }
  for(int ii=0; ii<Ph.ndof; ii++)
    {
      b[ii+UUh.ndof]=0;
    }
}
time2=clock()-time2;

real time3=clock();
xx = A^-1*b;
time3=clock()-time3;

for(int ii=0; ii<Uh.ndof; ii++)
    {
      u1[][ii]=xx[3*ii];
      u2[][ii]=xx[3*ii+1];
      u3[][ii]=xx[3*ii+2];
    }
  for(int ii=0; ii<Ph.ndof; ii++)
    {
      p[][ii]=xx[ii+UUh.ndof];
    }

timeI=clock()-timeI;

cout << "============= CPU TIME ============" << 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) medit("UV2 PV2",Th,[u1,u2,u3],p);

