Save fields under binary format

General Discussion
1

Hi every one,

I would like to export my data under binary format. I use

{ofstream export(“BF_binary.btxt”, binary);
export << Ux << endl;
export << Rd1 << endl;
}

but the resulting file BF_binary.btxt is not a binary file.

To overcome this, I tried to use medit using

savesol(“Medit.solb”,Th, [Ux,Uy,P], Rd1);

and then

uu = readsol(“Medit.sol”);

but the solution I read with readsol does not match the one I saved earlier. I see it by a plotting it.

Any solution please ?

Regards,

Adrien

2

Hi !

I have tried your first approach, namely

and it seems that indeed it saves the date in ordinary .txt format.

Regarding your second approach with “savesol”, you need to be more specific with your choice of “fespace” for “[Ux,Uy,P]”, what is “uu” (from which fespace) and what is Rd1 (integer or fe function) ? Here is also the reference for using “savesol” - Developers - I couldn’d really find the “readsol” in documentation (it says “Todo” in External libraries, I guess the documentation on this topic isn’t completed yet).

From Developers it appears that “order” in “savesol” has to be 0 or 1 (and no other) - are your finite element spaces all P1 or P0 ? This might be the first issue.

  1. I assume “int Rd1=1”, i.e. “Rd1” is “order” in your code above. The code below works well for me (only one vector finite element space)

     mesh Th = square(10,10);
 fespace Vh(Th,[P1,P1,P1]);
 Vh [vx,vy,p]=[x,y,x+y];

 plot([vx,vy],wait=1,cmm="velocity before");
 plot(p,fill=1,wait=1,cmm="pressure before");

 savesol("vec.solb",Th,[vx,vy,p],order=1);
 [vx,vy,p]=[0,0,0];

 vx[] = readsol("vec.solb");
 plot([vx,vy],wait=1,cmm="velocity after");
 plot(p,fill=1,wait=1,cmm="pressure after");

  2. In case your “Rd1” is finite element function from scalar fespace (P1). In this case you have two different fespaces of order P1, one vector fespace and one scalar. The following code seems to work well for me:

     mesh Th = square(10,10);
 fespace Uh(Th,P1);
 Uh uh=x+y; plot(uh,wait=1,fill=1,value=1,cmm="uh initial");

 fespace Vh(Th,[P1,P1]);
 Vh [vhx,vhy]=[x,y]; plot([vhx,vhy],wait=1,fill=1,value=1,cmm="[vhx,vhy] initial");

 savesol("mix.solb",Th,uh,[vhx,vhy],order=1);
 uh=0; [vhx,vhy]=[0,0];

 real[int] sols=readsol("mix.solb");
 for(int k=0; k<Uh.ndof; ++k){
         uh[][k]=sols[3*k];
         vhx[][2*k]=sols[3*k+1];
         vhx[][2*k+1]=sols[3*k+2];
 }

 plot(uh,wait=1,fill=1,value=1,cmm="uh after reading");
 plot([vhx,vhy],wait=1,fill=1,value=1,cmm="[vhx,vhy] after reading");

So, “savesol” saves your functions as

    uh[][k]   vhx[][2k]   vhx[][2k+1]              ... for k = 0 : Uh.ndof

note, Vh.ndof = 2*Uh.ndof. Also, you might want to look into how dofs are numbered for vector fespaces.

  1. As a general suggestion, its not a bad idea to create very coarse mesh and fespaces on it, and savesol in txt format and then you can personally investigate what is going on.

You can also check your options with vtk/vtu format (for visualizetion in Paraview) - External libraries. Documentation is very clear here.

if

3

Thank you very much for this precise answer!

In my case, [Ux, Uy, P] are fields respectively [P2, P2, P1] and Rd1 is an integer. The fact of using P2 fields is surely the source of my error then.

Adrien

4

Effectivelly,

for finite element [P2,P2,P1], to see which dof are asocial to field P2,P2 ou P1

   mesh Th= ...
   fespace Vh(Th,[P2,P2,P1]);
   fespace Uh(Th,P2);
   fespace Ph(Th,P1);
  Uh un;
  Ph pn;
  un[]=0:un[].n-1; //    value  0, 1, .... 
  pn[]=0:p[].n-1; // value  0, 1, .... 

  Vh [u1,u2,p]=[1,2,3];
  Vh [v1,v2,q]=[un,un,pn]; 

  cout << u1[] << endl; 
 if( u1[][dof] == 1) dof associated to u1  is v1[][dof] in P2 
 if( u1[][dof] == 2) dof associated to u2   is v1[][dof] in P2 
 if( u1[][dof] == 3) dof associated to u3  is  is u1[][dof] in P1