//  stokes_3d.edp
//
//  Discussion:
//
//    Solve a Stokes problem in a 3D cube.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/freefem++/stokes_3d.edp
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    05 December 2014
//
//  Reference:
//
//    Frederic Hecht,
//    Freefem++,
//    Third Edition, version 3.22
//

load "msh3" load "medit"  // dynamics load tools for 3d.

int nn = 8;
mesh Th2=square(nn,nn);
fespace Vh2(Th2,P2);  Vh2 ux,uz,p2;
//
//  Create the 3D mesh.
//
int[int] rup=[0,2];
int[int] rdown=[0,1];
int[int] rmid=[1,1,2,1,3,1,4,1];
real zmin=0;
real zmax=1;

mesh3 Th = buildlayers(Th2,nn,
  zbound=[zmin,zmax],  reffacemid=rmid, 
  reffaceup = rup,     reffacelow = rdown);
//
//  Display the 3D mesh using MEDIT.
//
medit ( "c10x10x10", Th );
//
//  Define VVh, the finite element space for 3D velocity + Pressure.
//
fespace VVh(Th,[P2,P2,P2,P1]);

macro Grad(u) [dx(u),dy(u),dz(u)] // EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3))  //EOM
//
//  Define (U,P) and (V,Q), elements of the VVH function space.
//
VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];

solve vStokes([u1,u2,u3,p],[v1,v2,v3,q]) 
  = int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) 
                       + Grad(u2)'*Grad(v2) 
                       +  Grad(u3)'*Grad(v3)
                       - div(u1,u2,u3)*q 
                       - div(v1,v2,v3)*p 
                       + 1e-10*q*p ) 
  + on(2,u1=1.,u2=0,u3=0) 
  + on(1,u1=0,u2=0,u3=0);

 plot(p,wait=1, nbiso=5);  // a 3d plot of iso  pressure. in progress... march 2009
 //  to see the 10 cut plan in 2d 

for(int i=1;i<10;i++)
{
 real yy=i/10.; // compute yy.
  // do 3d -> 2d interpolation.
 ux= u1(x,yy,y); uz= u3(x,yy,y);  p2= p(x,yy,y);
 plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1);
}