Bug in saving vector field to txt?

General Discussion
1

Hi,

I’m trying to save the 2D vector field from Stokes problem. However, no matter how I do it u2 and u1 are exactly the same. Weirdly, plot(u1) and plot(u2) still show correctly different values. Here is my script:

include "./freefem_matlab_octave_plot/release-v2.0/demos/ffmatlib.idp"

macro grad(u) [dx(u),dy(u)]//
macro Grad(u1,u2) [ grad(u1), grad(u2)]//
macro div(u1,u2) (dx(u1)+dy(u2))//


int nn = 30; // number of edge in each direction
// mesh Th=square(nn,nn,[2*pi*x,2*pi*y],flags=3);

real r=0.05;
real sep=1.12;
real L=3.0;
real rcx=L/2.0-r*sep;
real rcy=L/2.0;
real rcx2=L/2.0+r*sep;
real rcy2=rcy;

int[int] labs = [1,2,3,4];

// border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
border C01(t=0, 1){x=L*t; y=0; label=labs[0];}
border C02(t=0, 1){x=L; y=L*t; label=labs[1];}
border C03(t=0, 1){x=L*(1-t); y=L; label=labs[2];}
border C04(t=0, 1){x=0; y=L*(1-t); label=labs[3];}

border b(t=0, 2*pi){x=rcx+r*cos(t); y=rcy+r*sin(t); label=5;}
border b2(t=0, 2*pi){x=rcx2+r*cos(t); y=rcy2+r*sin(t); label=6;}

real factormsh = 3;
mesh Th = buildmesh(C01(nn) + C02(nn) + C03(nn) + C04(nn) 
    + b(-nn*factormsh) + b2(-nn*factormsh) );

mesh ThU=trunc(Th,1,split=2);  // Velocity mesh


fespace Uh(ThU,[P1,P1],periodic=[[2, y], [4, y], [1, x], [3, x]]);
fespace Ph(Th,P1,periodic=[[2, y], [4, y], [1, x], [3, x]]);
Uh [u1,u2], [v1,v2];
Ph p,q;
solve Stokes(<[u1,u2],[p]>,<[v1,v2],[q]>) =
int2d(ThU)( (Grad(u1,u2):Grad(v1,v2)) )
+ int2d(ThU)( - div(u1,u2)*q - div(v1,v2)*p )
+ int2d(Th)( -1e-10*p*q )
+ on(5,u1=-1,u2=0)
+ on(6,u1=1,u2=0);
cout << "u1 dof=" << u1[].n << endl;
for (int i=0; i<4; i++){
    cout << u1[][i] << endl;
}

cout << "u2 dof=" << u2[].n << endl;

for (int i=0; i<4; i++){
    cout << u2[][i] << endl;
}
plot(u1, cmm="u1" , fill=1,value=1);
plot(u2, cmm="u2" , fill=1,value=1);
string fnamedata="u.txt";
ffSaveData2(u1,u2,fnamedata)

So the cout lines give exactly same outputs. And ffSaveData2(u1,u2,fnamedata) also gives exactly same values. the txt file looks like:

1 -1
-9.376217058e-32 -9.376217058e-32
-1 -1
-1.202971747e-31 -1.202971747e-31
-1 -1
...
2

Hi,
when you declare
fespace Uh(ThU,[P1,P1],periodic=[[2, y], [4, y], [1, x], [3, x]]); Uh [u1,u2], [v1,v2];
You have to know that [u1,u2] is a whole object in itself containing all the components of u1 and u2 together.
When you write u1[] it gives the list of all components of the vector [u1,u2].
u2[] gives the same i.e. u2[]=u1[] (formally it represents [u1,u2][]!)
For the plot there is no problem, when you write u1 FreemFem understands to plot the first component of [u1,u2].
If you want to get the components of u1 alone you have to add

fespace UUh(ThU,[P1],periodic=[[2, y], [4, y], [1, x], [3, x]]);
UUh ux,uy;
ux=u1;
uy=u2;
cout << "ux dof=" << ux[].n << endl;
for (int i=0; i<2; i++){
    cout << ux[][i] << endl;
}
cout << "uy dof=" << uy[].n << endl;
for (int i=0; i<2; i++){
    cout << uy[][i] << endl;
}

The result is

u1 dof=53092
-1
-5.92836e-31
-1
-5.2406e-31
u2 dof=53092
-1
-5.92836e-31
-1
-5.2406e-31
ux dof=26546
-1
-1
uy dof=26546
-5.92836e-31
-5.2406e-31

You see that u1[] contains both ux[] and uy[] in an interlaced form.

1 Like
3

Please do not post issues both here and on GitHub.

4

Thank you! That solved my issue.

5

Sorry! Thanks for the reply.