Writing to file appears to get wrong values,

General Discussion
1

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)/ (nxnynz) /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 :slight_smile:
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=hhhhhh/6.);

Th=mmg3d(Th,metric=hx);

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

2

I Have do some small correction ,

  1. remove mes mesh in small 2 small cube

  2. no BC on pressure,

  3. add BC on all cube

  4. in this case Stoke equation is not the correct modele, I think it is Navier Stokes (More difficult)
    and Do not do this.

bbbb.edp (4.6 KB)

3

Thanks for looking at it in details. Yes, I’ve been moving around things like BC’s to
get some idea of what is going on and intended to solve NS but it was not clear
exactly how specific the solver had to be and how much it could figure out :slight_smile:
I think I copied that from an example with the NS weak form calling vStokes.

I’ll probably look at the c++ source eventually. I had to specify CG which I guess is
conjugate gradient to get it to fit into memory.

However, the immediate problem is that p.txt appears inconsistent with my code.

I’m wiriting this,
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;

but all the values are identical on each line,

more p.txt | awk ‘{print $4-$5+$6-$7}’ | sort | uniq -c
19375 0

Thanks.

4

ok, it looks like most of the “stuff” works and I can worry more about
math and physics :slight_smile: First, thanks for looking at the BC’s although in this
reference frame the velocity on the road should still be the same
as the two vehicles are stationary and rhe road is moving so no boundary
layer except whatever the cars create. Getting the elements from the
solution still puzzles me but I have a working approach. I guess in general
I’m concerned about silent failure modes but rhis works now,

real x=Th[i][0].x;
real y=Th[i][0].y;
real z=Th[i][0].z;
ff<<x<<" “<<y<<” “<<z
<<” " <<p(x,y,z)
<<" " <<ux(x,y,z)
<<" " <<uy(x,y,z)
<<" " <<uz(x,y,z)
<< endl;

I’ll look now into taking the vStokes thing to NS and check the physics. I did find I
could integrae pressure on the front and back land saw a benefit from proximity
although some pressures were negative :slight_smile:

real p1f=int2d(Th,100)( p ) ;
real p1r=int2d(Th,101)(p) ;

real p2f=int2d(Th,200)( p ) ;
real p2r=int2d(Th,201)(p) ;

cout<<" vars del “<<del<<” uxc "<<uxc<< " pressures “<<p1f<<” “<<p1r<<” del1 “<<(p1f-p1r)<<” “<<p2f<<” “<<p2r<<” del2 "<<(p2f-p2r)<<endl;

Thanks.

5

Oh and thanks for removing the elements from the inner cube.
I thought the meshing was supposed to be confined between surfaces
in the direction indicated by “orientation” but that does not seem to be
the case. When I exported the mesh points I did find a few inside the cubes
although I did not look in detail as you had apparently removed them with
trunc which seems to work well.

6

Remember in Stokes with Dirichlet BC , the pressure is defined through a constant.

7

Thanks. I think in the equation I was looking at only grad(p) appeared but
in the weak form there was a p. But the other issues was just verifying that
the BC’s were being enforced where expected. I think rather than go full NS,
I’m going to try an electormagnetics problem. There seem to be a lot of Helmholtz
examples.
It looks like “trunc” is pretty versatile not just for removing vertices but also
refining that mesh. I guess regions can be defined with no consequence ( don’t
produce automatic changes to the mesh ) until trunc is called to either elminate
or split elements in the region.

Presumably I can find a lot of places like points and corners that need
refinement and do some of that a priori by hand.

I think I saw one helmholtz example but it claimed to only work for
real k. But if i can figure out how to use it I also noted there is a Sommerfield radiation
boundary conditions which may be nice- another tutortial I was going to try
was radiation pattern of some metal blobs that kind of look like an LC tank :slight_smile: