So I realized that Stokes does not work and I needed vStokes to solve

my 2 cubes in the wind example ( drafting a semi truck lol ). The plot looks

plausible and I wanted to explore in R. However, all the velocity and pressure

values are identical on each point AFAICT. What is wrong where I write

to p.txt? Thanks.

///home/ubuntu/dev/freefem/FreeFem-sources-master/examples/3d$ vi Poisson-cube-ballon.edp

//marchywka@happy:/home/ubuntu/dev/freefem/FreeFem-sources-master/examples/3d$

// Navier-Stokes equations

load “msh3”

load “tetgen”

load “medit”

load “mshmet”

load “mmg”

//

real volumetet; // use in tetg.

meshS Thx0,Thx1, Thy0, Thy1, Thz0, Thz1;

int nx=10,ny=10,nz=10;

func meshS mycube ( real x0, real x1, real y0,real y1, real z0, real z1, int orient, int lb)

{

meshS ThHex;

// first build the 6 faces of the hex.

//real x0=-1,x1=1;

//real y0=-1.1,y1=1.1;

//real z0=-1.2,z1=1.2;

// a volume of on tet.

volumetet= (x1-x0)*(y1-y0)*(z1-z0)/ (nx*ny*nz) /6.;

mesh Thx = square(ny,nz,[y0+(y1-y0)*x,z0+(z1-z0)*y]);

mesh Thy = square(nx,nz,[x0+(x1-x0)*x,z0+(z1-z0)*y]);

mesh Thz = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

int ii=lb-1;

int[int] refz=[0,5+ii]; // bas

int[int] refZ=[0,6+ii]; // haut

int[int] refy=[0,3+ii]; // devant

int[int] refY=[0,4+ii]; // derriere

int[int] refx=[0,1+ii]; // gauche

int[int] refX=[0,2+ii]; // droite

int f1=-1;

int f2=1;

if (orient==1) { f1=1; f2=-1; }

int[int] l0=[0,lb+0];

//meshS

Thx0 = movemesh23(Thx,transfo=[x0,x,y],label=l0,orientation=f1,region=refx,removeduplicate=false);

++l0[1];

//meshS

Thx1 = movemesh23(Thx,transfo=[x1,x,y],label=l0,orientation=f2,region=refX,removeduplicate=false);

++l0[1];

//meshS

Thy0 = movemesh23(Thy,transfo=[x,y0,y],label=l0,orientation=f2,region=refy,removeduplicate=false);

++l0[1];

//meshS

Thy1 = movemesh23(Thy,transfo=[x,y1,y],label=l0,orientation=f1,region=refY,removeduplicate=false);

++l0[1];

//meshS

Thz0 = movemesh23(Thz,transfo=[x,y,z0],label=l0,orientation=f1,region=refz,removeduplicate=false);

++l0[1];

//meshS

Thz1 = movemesh23(Thz,transfo=[x,y,z1],label=l0,orientation=f2,region=refZ,removeduplicate=false);

cout<<“l0=”<<l0<<endl;

//medit(" — ", Thx0,Thx1,Thy0,Thy1,Thz0,Thz1);

ThHex = Thx0+Thx1+Thy0+Thy1+Thz0+Thz1;

return ThHex;

}

real x0=-1,x1=1; real yy0=-1.1,yy1=.1; real z0=-1.2,z1=.2;

real xa0=-.5; real ya0=-.8; real za0=-1;

real xb0=-.1; real yb0=-.8; real zb0=-1;

real sz=.25;

meshS ThHex = mycube(x0,x1,yy0,yy1,z0,z1,1,1);

real vol=volumetet;

meshS Thsrc=Thx0;

nx=10; ny=10; nz=10;

meshS Tha = mycube(xa0,xa0+sz,ya0,ya0+sz,za0,za0+sz,1,100);

meshS Thb = mycube(xb0,xb0+sz,yb0,yb0+sz,zb0,zb0+sz,1,200);

meshS Thc= Tha+Thb+ThHex;

//medit(" asdf ",Thc);

real[int] domaine = [-.9,-1,-1.1,.5,vol];

mesh3 Th = tetg(Thc,switch=“pqaAAYYQ”,nbofregions=1,regionlist=domaine);

for( int i=0; i<1; ++i)

{

cout<<" ++++++++++++++++++++++++++++++++++++ "<<i<<endl;

real uxc=100;

// vStokes example a lot odifferent doh

// example uses 1,2,3 instead of x,y,z

fespace VVh(Th, [P2, P2, P2, P1]);

VVh [ux, uy, uz, p];

VVh [vx, vy, vz, q];

fespace Xh(Th,P1);

// this does not fing work.

//savemesh(Th,“./ux.txt”,[x,y,z,ux]);

{ofstream ff(“dump.txt”); dumptable(ff); }

macro grad(u) [dx(u),dy(u),dz(u)]// def of grad operator

macro Grad(u) [dx(u), dy(u), dz(u)] //

macro div(ux,uy,uz) (dx(ux) + dy(uy) + dz(uz)) //

solve vStokes ([ux, uy, uz, p], [vx, vy, vz, q],solver=CG)

= int3d(Th, qforder=3)(

Grad(ux)’ * Grad(vx)

- Grad(uy)’ * Grad(vy)
- Grad(uz)’ * Grad(vz)

- div(ux, uy, uz) * q
- div(vx, vy, vz) * p

- 1e-10 * q * p

)

+on(1,ux=uxc)

+on(2,ux=uxc)

//+on(3,ux=uxc) +on(4,ux=uxc) +on(5,ux=uxc) +on(6,ux=uxc)

+on(1,p=1000)

+on(2,p=2000)

+on(100,101,102,103,104,105,200,201,202,203,204,205,ux=0, uy=0,uz=0)

;

cout<<" ++++++++++++ done solve +++++++++++++++++++ "<<i<<endl;

//real[int] foo(Th.nv);

{

// wtf?

//{ofstream ff(“u.txt”); ff<<ux<<uy<<uz; }

//{ofstream ff(“p.txt”); ff<<p; }

{ofstream ff(“p.txt”);

// probably this is number of triangles or tets

for (int i = 0; i < Th.nt; i++)

{

for (int j = 0; j < 1+(3*0); j++)

ff << Th[i][j].x

<< " "<< Th[i][j].y

<< " "<< Th[i][j].z

<< " " << p[VVh(i,j)]

<< " " << ux[VVh(i,j)]

<< " " << uy[VVh(i,j)]

<< " " << uz[VVh(i,j)]

<< endl;

}

}

} // scoping

//real[int] viso=[0,.01,.05,.1,.2,.3];

//plot(ux, viso=viso, value=true, wait=1);

//plot([ux,uy,uz],p, value=true, wait=1);

plot(p, value=true, wait=1);

plot(ux, value=true, wait=1);

Xh hh=p; real lerr=.001;

real[int] hx=mshmet(Th,(dx(p)*dx(p)+dy(p)*dy(p)+dz(p)*dz(p)),normalization=1,aniso=0,nbregul=1,hmin=Th.measure/80,hmax=Th.measure/2,err=lerr);

cout<<" +++++++++ RECO ++++++++++++++++++++ "<<i<<endl;

// zero regions?

//Th=tetgreconstruction(Th,switch=“raAQ”,nbofregions=1,regionlist=domaine,sizeofvolume=hh*hh*hh/6.);

Th=mmg3d(Th,metric=hx);

medit(“tetg”,Th,wait=1);

} //i