[ff3d-users] help!

[Augusto Almeida]

I need run a enviroment (a quadratic tube), with a viscous fluid in there, 
the aim is steady state and export the vector components (discretization 
70x70x70) in a archive with this format :

vx vy vz

The archive contais 343000 registers.

Follow the code (sorry some parts is in portuguese!) :

Thanks too much!

double limite_inf        =  0.0;
double limite_inf_r      =  0.0;
double limite_sup        =  1.0;
double limite_sup_r      =  1.0;
double dimensao          = 20.0;
double dimensao_refinada = 70.0;

// Definindo a dimensao do solido (deve ser menor para evitar calculos
// desnecessarios pois dentro do solido a velocidade e zero)
vector normal   = (dimensao, dimensao, dimensao);

// Definindo a dimensao da malha onde o fluido percorre
vector refinado = (dimensao_refinada, dimensao_refinada, dimensao_refinada);

// Definindo os limites superiores e inferiores da "caixa" onde ocorre a
// simulacao
vertex ini   = (limite_inf  , limite_inf  , limite_inf);
vertex ini_r = (limite_inf_r, limite_inf_r, limite_inf_r);
vertex fin   = (limite_sup  , limite_sup  , limite_sup);
vertex fin_r = (limite_sup_r, limite_sup_r, limite_sup_r);

// Criando as malhas
mesh MALHA          = structured(normal  , ini, fin);
mesh MALHA_REFINADA = structured(refinado, ini_r, fin_r);

// Carregando o arquivo com as geometrias usadas no experimento
scene placas = pov("tubo_quadrado.pov");

// Definindo o dominio fora (outside) dos solidos
domain DOMINIO = domain(placas, outside(<0,0,1>));

// Definindo solido com refinamento inferior ao da malha
mesh obstacle = surface(DOMINIO, MALHA);

// Funcao de penalizacao para nao realizar calculo dentro do solido
function P = 1E3 * (one(<0,0,1>)) + 1;

// Calculo da area
double area = int[MALHA](1);

function          uin = 0.001;
function           mu = 0.1;
femfunction u0(MALHA) = 0;
femfunction v0(MALHA) = 0;
femfunction w0(MALHA) = 0;
femfunction p0(MALHA_REFINADA) = 0;

double dt = 0.01;
double pp;
double qq;

double i = 1;

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

do{
   cout << "========== Tubo, Navier-Stokes passo " << i << "\n";

   // Calculo da componente x
   solve(u) in MALHA
   memory(matrix=none),
   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 MALHA xmin;
        u = 0   on MALHA ymin;
        u = 0   on MALHA ymax;
        u = 0   on MALHA zmin;
        u = 0   on MALHA zmax;
   };

   // Calculo da componente y
   solve(v) in MALHA
   memory(matrix=none),
   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 MALHA xmin;
        v = 0 on MALHA ymin;
        v = 0 on MALHA ymax;
        v = 0 on MALHA zmin;
        v = 0 on MALHA zmax;
   };

   // Calculo da componente z
   solve(w) in MALHA
   memory(matrix=none),
   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 MALHA xmin;
        w = 0 on MALHA ymin;
        w = 0 on MALHA ymax;
        w = 0 on MALHA zmin;
        w = 0 on MALHA zmax;
   };

   qq = int[MALHA](dx(u) + dy(v) + dz(w)) / area;
   solve(q) in MALHA_REFINADA
     krylov(type=cg,precond=ichol)
   {
        pde(q)
           - div(dt * grad(q)) = (dx(u) + dy(v) + dz(w) - qq);

        q = 0 on MALHA_REFINADA xmax;
   };

   p0 = p0 - q;
   pp = int[MALHA](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 <= 100);

// transforma o objeto em uma figura tetraédrica
mesh T = tetrahedrize(DOMINIO, MALHA_REFINADA);

// salva a malha
//save(vtk, "ns_placas",[u,v,w],MALHA_REFINADA);

//save(raw,"veloctuboquadrado00.vec", [u, v], MALHA_REFINADA);
//save(raw,"veloctuboquadrado003d.vec", [u, v, w], MALHA_REFINADA);
//save(medit,"tubo_quadrado_vel",T);
//save(medit,"tubo_quadrado_vel",[u,v,w],T);
//save(medit,"tubo_quadrado_pressao",T);
//save(medit,"tubo_quadrado_pressao",p0,T);


.pov
union{
     box {<0.0,0.0,0.0> , <1.0,1.0,0.2>}
     box {<0.0,0.0,0.2> , <1.0,0.2,0.8>}
     box {<0.0,0.8,0.2> , <1.0,1.0,0.8>}
     box {<0.0,0.0,0.8> , <1.0,1.0,1.0>}
     pigment{ color rgb <0,0,1>}
     }

_________________________________________________________________
Mande torpedos SMS do seu messenger para o celular dos seus amigos 
http://mobile.msn.com/