// -*- 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: navier-stokes.ff,v 1.6 2007/05/21 21:23:17 delpinux Exp $

// Navier-Stokes flow around a cylinder
// Uses a projection scheme over the velocity divergence free space.

// Fictitious domain method uses volume penalty method

// Only 5 iterations are performed

vector n2 = (25, 5, 5);
vector n1 = 2*n2;

vertex A = (0,0,0);
vertex B = (5,1,1);
mesh M = structured(n1,A,B);
mesh M2= structured(n2,A,B);
scene S = pov("navier-stokes.pov");
domain Omega = domain(S,outside(<1,0,0>));

mesh obstacle = surface(Omega,M);

function P = 1E3*one(<1,0,0>)+1;

// Computes the area of Omega
double area = int[M](1);

function uin = 12.5;
function mu = 0.1;
femfunction u0(M) = 0;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M2) = 0;

double dt = 0.01;
double pp;
double qq;

double i = 1;

function u = 0;
function v = 0;
function w = 0;

do{

    cout << "========== Navier-Stokes step " << i << "\n";

    solve(u) in M
      krylov(type=cg,precond=ichol)
    {
	pde(u)
	    (P/dt)*u - div(mu*grad(u))
	    = ((1/dt)*convect([u0,v0,w0],dt,u0) - dx(p0));
	u = uin on M xmin;
	u = 0 on M ymin;
	u = 0 on M ymax;
	u = 0 on M zmin;
	u = 0 on M zmax;
    };

    solve(v) in M
    krylov(type=cg,precond=ichol)
    {
      	pde(v)
	    (P/dt)*v - div(mu*grad(v))
	    = ((1/dt)*convect([u0,v0,w0],dt,v0) - dy(p0));
	v = 0 on M xmin;
	v = 0 on M ymin;
	v = 0 on M ymax;
	v = 0 on M zmin;
	v = 0 on M zmax;
    };

    solve(w) in M
      krylov(type=cg,precond=ichol)
    {
	pde(w)
	    (P/dt)*w - div(mu*grad(w))
	    = ((1/dt)*convect([u0,v0,w0],dt,w0) - dz(p0));
	w = 0 on M xmin;
	w = 0 on M ymin;
	w = 0 on M ymax;
	w = 0 on M zmin;
	w = 0 on M zmax;
    };

    qq = int[M](dx(u)+dy(v)+dz(w))/area;
    solve(q) in M2
      krylov(type=cg,precond=ichol)
    {

	pde(q)
	    - div(dt*grad(q)) = (dx(u)+dy(v)+dz(w) - qq);

	q = 0 on M xmax;
    };

    p0 = p0 - q;
    pp = int[M](p0)/area;
    p0 = p0 - pp;
    u0 = u + (dt*dx(q));
    v0 = v + (dt*dy(q));
    w0 = w + (dt*dz(q));

    i = i+1;

} while(i <= 10);

// Final results.

mesh T = tetrahedrize(Omega,M);

save(vtk,"ns",[u,v,w],M);

save(medit,"velocity",T);
save(medit,"velocity",[u,v,w],T);
save(medit,"pressure",T);
save(medit,"pressure",p0,T);