// File Stokes3D.edp
//
// 3D STOKES EQUATIONS
//

load "msh3"
		
load "medit"
		
//
// NUMERICAL VALUES TO PROVIDE
//

verbosity = 10;

int NSigma=100, NBase=99, NTop=101, Nin=200, Nout=201;
int N1=9, N2=3, Nxy2=5, Nxy1=9, Nz=9;
int imesh=0;

real xnu=.1, coefin=1.;
real errmesh=1.e-2;

//
// AUXILIAR VARIABLES
//

real pis2=.5*pi, errold;

//
// THE 2D TRIANGULATION
//

border C1(s=0., 3.){x=s; y=0.; label=NSigma;}
border C2(t=-pis2, 0.){x=3.+2.*cos(t); y=2.+2.*sin(t); label=NSigma;}
border C3(s=2., 5.){x=5.; y=s; label=NSigma;}
border C4(s=5., 4.){x=s; y=5.; label=Nout;}
border C5(s=5., 2.){x=4.; y=s; label=NSigma;}
border C6(t=0., -pis2){x=3.+cos(t); y=2.+sin(t); label=NSigma;}
border C7(s=3., 0.){x=s; y=1.; label=NSigma;}
border C8(s=1., 0.){x=0.; y=s; label=Nin;}

plot(C1(N1)+C2(Nxy1)+C3(N1)+C4(N2)+C5(N1)+C6(Nxy2)+C7(N1)+C8(N2),
    wait=true);

mesh Th2d=buildmesh(C1(N1)+C2(Nxy1)+C3(N1)+C4(N2)+C5(N1)+C6(Nxy2)+C7(N1)+C8(N2));
plot(Th2d, wait=true);
//savemesh(Th2d,"E-Mesh2D.msh");

/*
errold=errmesh;
for (int i=0;i< 4;i++)
{
	imesh=i+1;
	Th2d=adaptmesh(Th2d,f2d0,err=errold);
	plot(Th2d, cmm="2D MESH NO. "+imesh, wait=true);
	errold = .5*errold;
};
*/

//
// THE 3D TRIANGULATION
//

int refint = Th2d(1.,.5).region;
		
int[int] rup=[refint,NTop,NSigma,NSigma,Nin,Nin,Nout,Nout],    // upper face --> NTop
		 rdown=[refint,NBase,NSigma,NSigma,Nin,Nin,Nout,Nout], // lower face --> NBase
		 rmid=[NSigma,NSigma,Nin,Nin,Nout,Nout];               // lateral faces, NSigma--> NSigma

mesh3 Th3d = buildlayers(Th2d, Nz, zbound=[0.,3.],
		 reffaceup=rup, reffacelow=rdown, reffacemid=rmid);
plot(Th3d,wait=true);
//plot(Th3d,wait=1,ps="mixedTh3d.eps");
//savemesh(Th3d,"E-Mesh3D.mesh");

medit("3D Mesh", Th3d);

//
// THE FUNCTIONS
//

func uin1=coefin*y*(1.-y)*z*(3.-z);
func uin2=0.;
func uin3=0.;

func fseg1=0.;
func fseg2=0.;
func fseg3=-1.;


//
// 	THE SPACES
//

fespace Xh(Th3d,P2);
Xh u, v, w, u1, u2, u3, v1, v2, v3, f1=fseg1, f2=fseg2, f3=fseg3;

fespace Mh(Th3d,P1);
Mh pp, q;

//
// THE PROBLEM
//

macro Prod(u,v) dx(u)*dx(v)+dy(u)*dy(v)+dz(u)*dz(v)  // EOM
macro Grad(q) [dx(q),dy(q),dz(q)]  // EOM
macro Div(u,v,w) (dx(u)+dy(v)+dz(w))  // EOM

problem Stokes([u1,u2,u3,pp],[v1,v2,v3,q]) = int3d(Th3d)(
        xnu*(Prod(u1,v1)+Prod(u2,v2)+Prod(u3,v3)) + [v1,v2,v3]'*Grad(pp)
        + Div(u1,u2,u3)*q + 1.e-5*pp*q)
		- int3d(Th3d)(f1*v1+f2*v2+f3*v3)
        + on(Nin, u1=uin1,u2=uin2,u3=uin3)
        + on(NTop,NBase,NSigma, u1=0.,u2=0.,u3=0.)
//        + on(NBase,NSigma,NTop, u=gseg)
;

// SOLVING THE PROBLEM

Stokes;

//
// DISPLAYING, WRITING AND PRINTING
//

// 3D VIEWS

plot(Th3d, u1, wait=true, value=true, fill=false, cmm="THE SOLUTION");
plot(Th3d, u2, wait=true, value=true, fill=false, cmm="THE SOLUTION");
plot(Th3d, u3, wait=true, value=true, fill=false, cmm="THE SOLUTION");

plot(Th3d, pp, wait=true, value=true, fill=false, cmm="THE SOLUTION");

// ADDITIONAL (HORIZONTAL) 2D VIEWS

plot(Th2d, wait=true);

fespace Vh2d(Th2d,P2); Vh2d uxy1, uxy2;

real hz = .3, zz;
for(int i=0;i<=10;i++)
{
	zz = i*hz;
	uxy1 = u1(x,y,zz);
	uxy2 = u2(x,y,zz);
    plot(uxy1, cmm="CUT AT z ="+zz, value=true, fill=true, wait=true);
    plot(uxy2, cmm="CUT AT z ="+zz, value=true, fill=true, wait=true);
    plot([uxy1,uxy2], cmm="CUT AT z ="+zz, value=true, fill=true, wait=true);
}

//{ // save solution
//ofstream file("Solution3D.txt");
//file << u[] << endl;
//} // close the file (end block)

// { // read
// ifstream file("f.txt");
// file >> g[] ;
// } // close reading file (end block)