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


// other
load "medit"
load "MUMPS_mpi"

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


int nn=5;
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);

real ccc = mpisize/real(Th.nt) ;
Th=change(Th,fregion= min(mpisize-1,int(nuTriangle* ccc+1e-10)));
if(mpirank==0)
for(int i=0;i<=mpisize;++i)
	 cout<< " Volume region  " <<  i << " " << int3d(Th,i)(1.) << endl;
//if(mpirank==0) 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;

//VVh [uu1,uu2,uu3,up];
real timeI=clock();
real time1=clock();
varf vlaplace1(u1,v1) = 
  int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) ) + on(2,u1=fup) + on(1,u1=0); //' for emacs 

varf vlaplacenull(u2,v2) = 
  int3d(Th,qforder=3)( Grad(u2)'*Grad(v2) ) + on(2,u2=0) + on(1,u2=0); //' for emacs 

varf vdx(u1,q) = int3d(Th,qforder=3)( -dx(u1)*q );
varf vdy(u2,q) = int3d(Th,qforder=3)( -dy(u2)*q ); 
varf vdz(u3,q) = int3d(Th,qforder=3)( -dz(u3)*q );

varf vmass(p,q) = int3d(Th,qforder=3)( 1e-10*p*q );

matrix M     = vmass(Ph,Ph,tgv=ttgv);
matrix L1    = vlaplace1(Uh,Uh,tgv=ttgv);
matrix Lnull = vlaplacenull(Uh,Uh);
matrix Bx = vdx(Uh,Ph);
matrix By = vdy(Uh,Ph);
matrix Bz = vdz(Uh,Ph);


matrix A = [ [ L1,     0,     0, Bx' ],     
	         [  0, Lnull,     0, By' ],
             [  0,     0, Lnull, Bz' ],
             [ Bx,    By,    Bz,  M  ] ];    //' for emacs


time1=clock()-time1;

real time2=clock();
real[int] b(VVh.ndof);
real[int] xx(VVh.ndof);
{ 
  real[int] b1 = vlaplace1(0,Uh,tgv=ttgv);
  real[int] b2 = vlaplacenull(0,Uh,tgv=ttgv);
  
  for(int ii=0; ii<Uh.ndof; ii++)
    {
      b[ii]=b1[ii];
      b[ii+Uh.ndof]=b2[ii];
      b[ii+2*Uh.ndof]=b2[ii];
    }
  for(int ii=UUh.ndof; ii<UUh.ndof+Ph.ndof; ii++)
    {
      b[ii]=0;
    }
}
time2=clock()-time2;

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

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

for(int ii=0; ii<Uh.ndof; ii++)
    {
      u1[][ii]=xx[ii];
      u2[][ii]=xx[ii+Uh.ndof];
      u3[][ii]=xx[ii+2*Uh.ndof];
    }
  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);

