# Writing to file appears to get wrong values,

**URL:** <https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224>\
**Category:** General Discussion\
**Created:** [January 16, 2023, 5:45pm UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224 "2023-01-16T17:45:29Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [January 16, 2023, 5:45pm UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224/1 "2023-01-16T17:45:29Z")

</div>

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](https://doc.freefem.org/models/navier-stokes-equations.html)

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

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [January 17, 2023, 3:09pm UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224/2 "2023-01-17T15:09:39Z")

</div>

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](https://community.freefem.org/uploads/short-url/stYNVEqFwf203RKnRhbMydIa6Bh.edp) (4.6 KB)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [January 17, 2023, 4:29pm UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224/3 "2023-01-17T16:29:30Z")

</div>

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 🙂  
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.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [January 18, 2023, 1:01am UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224/4 "2023-01-18T01:01:07Z")

</div>

ok, it looks like most of the “stuff” works and I can worry more about  
math and physics 🙂 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 🙂

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.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [January 18, 2023, 11:12am UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224/5 "2023-01-18T11:12:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [January 18, 2023, 8:59pm UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224/6 "2023-01-18T20:59:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [January 19, 2023, 12:10am UTC](https://community.freefem.org/t/writing-to-file-appears-to-get-wrong-values/2224/7 "2023-01-19T00:10:33Z")

</div>

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 🙂