// -*- c++ -*-

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2003  Laboratoire J.-L. Lions UPMC Paris

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: stokes.ff,v 1.5 2005/11/28 00:54:32 delpinux Exp $

// Stokes equations

vector n = (10,10,10);
vertex A = (0,0,0);
vertex B = (1,1,1);

mesh M = structured(n,A,B);

function teta0 = 0;
double eps = 1e-2;

function u1n=0;
function u2n=0;
function u3n=0;

function u1=0;
function u2=0;
function u3=0;
function p =0;

double i=0;
do {
   
solve(u1,u2,u3,p) in M
    krylov(precond=diagonal,type=bicg),
    cg(epsilon=1e-5) {

	test(U1,U2,U3,P)
	  int(grad(u1)*grad(U1) - p*dx(U1)
	      +grad(u2)*grad(U2) - p*dy(U2)
	      +grad(u3)*grad(U3) - p*dz(U3)
	      +u1n*dx(u1)*U1+u2n*dx(u2)*U2+u3n*dx(u3)*U3
	      -eps*grad(p)*grad(P) - eps*p*P
	      -dx(u1)*P - dy(u2)*P - dz(u3)*P) = 0;

	u1 = 0 on M xmin;
	u1 = 0 on M xmax;
	u1 = 0 on M ymin;
	u1 = 1 on M ymax;
	u1 = 0 on M zmin;
	u1 = 0 on M zmax;

	u2 = 0 on M;

	u3 = 0 on M;
    };

   u1n=u1;
   u2n=u2;
   u3n=u3;
   
   i++;
} while (i<=3);

// Post-processing
mesh T = tetrahedrize(M);
save(medit,"velocity",T);
save(medit,"velocity",[u1, u2, u3],T);
save(medit,"pressure",T);
save(medit,"pressure",p,T);

save(vtk,"ns",{[u1,u2,u3],p},T);