Trisurf Monte Carlo simulator
Samo Penic
2020-04-09 4f5ffcfd5ba38b6b7dd6d3b15de8bb8677537b9f
src/io.c
@@ -17,7 +17,7 @@
#include <dirent.h>
#include <errno.h>
#include <snapshot.h>
#include <b64zlib_compression.h>
ts_bool parse_args(int argc, char **argv){
    int c, retval;
@@ -245,6 +245,7 @@
   if(vesicle->filament_list!=NULL){
      if(vesicle->filament_list->poly[0]!=NULL){
      filno=vesicle->filament_list->n;
      fprintf(stderr,"WAS HERE AND I SHOULDNT BE\n");
      fonono=vesicle->filament_list->poly[0]->vlist->n;
      fil=1;
      }
@@ -255,16 +256,21 @@
   fprintf(fh, " <UnstructuredGrid>\n");
    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n);
    fprintf(fh,"<PointData Scalars=\"vertices_idx\">\n<DataArray type=\"Int64\" Name=\"vertices_idx\" format=\"binary\">");
   int *int_vector=(int *)malloc((vlist->n+monono*polyno+fonono*filno)*sizeof(ts_uint));
   long *int_vector=(long *)malloc((vlist->n+monono*polyno+fonono*filno)*sizeof(long));
   ts_double *float_vector=(double *)malloc((vlist->n+monono*polyno+fonono*filno)*sizeof(double));
   ts_double *float_cell_vector=(double *)malloc((blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n)*sizeof(double));
   ts_double *coordinate_vector=(double *)malloc(3*(vlist->n+monono*polyno+fonono*filno)*sizeof(double));
// sizes are as expected 4 bytes for int and 1 byte for char
//   printf("size of ts_uint %ld, of int %ld, of char %ld",sizeof(ts_uint), sizeof(int), sizeof(char));
   int offset=0;
      for(i=0;i<vlist->n;i++){
   //   fprintf(fh,"%u ",vtx[i]->idx);
      int_vector[i+offset]=vtx[i]->idx;
    }
   offset=offset+i;
   //polymeres
   if(poly){
      poly_idx=vlist->n;
      offset=offset+i;
      for(i=0;i<vesicle->poly_list->n;i++){
         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
            //fprintf(fh,"%u ", poly_idx);
@@ -285,26 +291,31 @@
         offset=offset+j;
      }
   }
   char *printout=ts_compress_intlist(int_vector,(vlist->n+monono*polyno+fonono*filno)*sizeof(ts_uint));
   printf("Offset in the end is %d, should be %d",offset,(vlist->n+monono*polyno+fonono*filno) );
   char *printout=ts_compress((char *)int_vector,offset*sizeof(long));
   fprintf(fh,"%s",printout);
   free(printout);
       fprintf(fh,"</DataArray>\n");
   if(cstlist!=NULL){
      fprintf(fh,"<DataArray type=\"Int64\" Name=\"vertices_in_cluster\" format=\"ascii\">");
      for(i=0;i<vlist->n;i++){
      fprintf(fh,"<DataArray type=\"Int64\" Name=\"vertices_in_cluster\" format=\"binary\">");
      offset=0;
      for(i=0;i<vlist->n;i++,offset++){
         if(vtx[i]->cluster!=NULL){
            fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
//            fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
            int_vector[offset]=vtx[i]->cluster->nvtx;
         } else {
            fprintf(fh,"-1 ");
//            fprintf(fh,"-1 ");
            int_vector[offset]=-1;
         }
          }
      //polymeres
      if(poly){
         poly_idx=vlist->n;
         for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
               fprintf(fh,"-1 ");
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
               //fprintf(fh,"-1 ");
               int_vector[offset]=-1;
            }
         }
      }
@@ -312,12 +323,15 @@
      if(fil){
         poly_idx=vlist->n+monono*polyno;
         for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
      //   fprintf(stderr,"was here\n");
               fprintf(fh,"-1 ");
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
   //               fprintf(fh,"-1 ");
               int_vector[offset]=-1;
            }
         }
      }
      char *printout=ts_compress((char *)int_vector,offset*sizeof(long));
      fprintf(fh,"%s",printout);
      free(printout);
      fprintf(fh,"</DataArray>\n");
@@ -325,16 +339,19 @@
   }
   //here comes additional data as needed. Currently only spontaneous curvature
   fprintf(fh,"<DataArray type=\"Float64\" Name=\"spontaneous_curvature\" format=\"ascii\">");
   for(i=0;i<vlist->n;i++){
      fprintf(fh,"%.17e ",vtx[i]->c);
   fprintf(fh,"<DataArray type=\"Float64\" Name=\"spontaneous_curvature\" format=\"binary\">");
   offset=0;
   for(i=0;i<vlist->n;i++, offset++){
//      fprintf(fh,"%.17e ",vtx[i]->c);
      float_vector[offset]=vtx[i]->c;
   }
      //polymeres
      if(poly){
         poly_idx=vlist->n;
         for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
               fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->c);
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
               //fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->c);
               float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->c;
            }
         }
      }
@@ -342,25 +359,34 @@
      if(fil){
         poly_idx=vlist->n+monono*polyno;
         for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
      //   fprintf(stderr,"was here\n");
               fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->c);
               //fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->c);
               float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->c;
            }
         }
      }
      printout=ts_compress((char *)float_vector,offset*sizeof(double));
      fprintf(fh,"%s",printout);
      free(printout);
    fprintf(fh,"</DataArray>\n");
   //here comes additional data. Energy!
   fprintf(fh,"<DataArray type=\"Float64\" Name=\"bending_energy\" format=\"ascii\">");
   for(i=0;i<vlist->n;i++){
      fprintf(fh,"%.17e ",vtx[i]->energy*vtx[i]->xk);
   fprintf(fh,"<DataArray type=\"Float64\" Name=\"bending_energy\" format=\"binary\">");
   offset=0;
   for(i=0;i<vlist->n;i++,offset++){
//      fprintf(fh,"%.17e ",vtx[i]->energy*vtx[i]->xk);
      float_vector[offset]=vtx[i]->energy*vtx[i]->xk;
   }
      //polymeres
      if(poly){
         poly_idx=vlist->n;
         for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
               fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k);
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
               //fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k);
               float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k;
            }
         }
      }
@@ -368,68 +394,108 @@
      if(fil){
         poly_idx=vlist->n+monono*polyno;
         for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
      //   fprintf(stderr,"was here\n");
               fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k);
//               fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k);
               float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k;
            }
         }
      }
   printout=ts_compress((char *)float_vector,offset*sizeof(double));
      fprintf(fh,"%s",printout);
      free(printout);
    fprintf(fh,"</DataArray>\n");
   
   fprintf(fh,"</PointData>\n<CellData>\n");
   if(vesicle->tape->stretchswitch==1){
      fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"ascii\">");
      for(i=0;i<blist->n;i++){
         fprintf(fh, "0.0 ");
      fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"binary\">");
      offset=0;
      for(i=0;i<blist->n;i++,offset++){
         //fprintf(fh, "0.0 ");
         float_cell_vector[offset]=0.0;
      }
      for(i=0;i<monono*polyno+filno*(fonono-1);i++){
         fprintf(fh,"0.0 ");
      for(i=0;i<monono*polyno+filno*(fonono-1);i++,offset++){
         float_cell_vector[offset]=0.0;
         //fprintf(fh,"0.0 ");
      }
      for(i=0;i<vesicle->tlist->n;i++){
         fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
      for(i=0;i<vesicle->tlist->n;i++,offset++){
         float_cell_vector[offset]=vesicle->tlist->tria[i]->energy;
         //fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
      }
      printout=ts_compress((char *)float_cell_vector,offset*sizeof(double));
      fprintf(fh,"%s",printout);
      free(printout);
       fprintf(fh,"</DataArray>\n");
   }
   fprintf(fh,"</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
   for(i=0;i<vlist->n;i++){
      fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
   fprintf(fh,"</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"binary\">");
   offset=0;
   for(i=0;i<vlist->n;i++,offset+=3){
      //fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
      coordinate_vector[offset]=vtx[i]->x;
      coordinate_vector[offset+1]=vtx[i]->y;
      coordinate_vector[offset+2]=vtx[i]->z;
   }
   //polymeres
   if(poly){
      for(i=0;i<vesicle->poly_list->n;i++){
         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
            fprintf(fh,"%.17e %.17e %.17e\n", vesicle->poly_list->poly[i]->vlist->vtx[j]->x,vesicle->poly_list->poly[i]->vlist->vtx[j]->y, vesicle->poly_list->poly[i]->vlist->vtx[j]->z );
         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,offset+=3){
            //fprintf(fh,"%.17e %.17e %.17e\n", vesicle->poly_list->poly[i]->vlist->vtx[j]->x,vesicle->poly_list->poly[i]->vlist->vtx[j]->y, vesicle->poly_list->poly[i]->vlist->vtx[j]->z );
      coordinate_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->x;
      coordinate_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[j]->y;
      coordinate_vector[offset+2]=vesicle->poly_list->poly[i]->vlist->vtx[j]->z;
         }
      }
   }
   //filaments
   if(fil){
      for(i=0;i<vesicle->filament_list->n;i++){
      for(i=0;i<vesicle->filament_list->n;i++,offset+=3){
         for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            fprintf(fh,"%.17e %.17e %.17e\n", vesicle->filament_list->poly[i]->vlist->vtx[j]->x,vesicle->filament_list->poly[i]->vlist->vtx[j]->y, vesicle->filament_list->poly[i]->vlist->vtx[j]->z );
   //      fprintf(fh,"%.17e %.17e %.17e\n", vesicle->filament_list->poly[i]->vlist->vtx[j]->x,vesicle->filament_list->poly[i]->vlist->vtx[j]->y, vesicle->filament_list->poly[i]->vlist->vtx[j]->z );
      coordinate_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->x;
      coordinate_vector[offset+1]=vesicle->filament_list->poly[i]->vlist->vtx[j]->y;
      coordinate_vector[offset+2]=vesicle->filament_list->poly[i]->vlist->vtx[j]->z;
         }
      }
   }
   printout=ts_compress((char *)coordinate_vector,offset*sizeof(double));
      fprintf(fh,"%s",printout);
      free(printout);
    fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
   for(i=0;i<blist->n;i++){
         fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
    fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"binary\">");
   free(int_vector);
   int_vector=(long *)malloc(((blist->n+monono*polyno+filno*(fonono-1))*2+vesicle->tlist->n*3)*sizeof(long));
   offset=0;
   for(i=0;i<blist->n;i++,offset+=2){
//         fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
         int_vector[offset]=blist->bond[i]->vtx1->idx;
         int_vector[offset+1]=blist->bond[i]->vtx2->idx;
   }
   //polymeres
   if(poly){
      poly_idx=vlist->n;
      for(i=0;i<vesicle->poly_list->n;i++){
         for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
//            fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx);
            fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono);
         for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++,offset+=2){
            //fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono);
            int_vector[offset]= vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono;
            int_vector[offset+1]=vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono;
         }
   //grafted bonds
      fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->grafted_vtx->idx, vesicle->poly_list->poly[i]->vlist->vtx[0]->idx+vlist->n+i*monono);
      //fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->grafted_vtx->idx, vesicle->poly_list->poly[i]->vlist->vtx[0]->idx+vlist->n+i*monono);
      int_vector[offset]=vesicle->poly_list->poly[i]->grafted_vtx->idx;
      int_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[0]->idx+vlist->n+i*monono;
      offset+=2;
      }
   }
@@ -438,32 +504,62 @@
   if(fil){
      poly_idx=vlist->n+monono*polyno;
      for(i=0;i<vesicle->filament_list->n;i++){
         for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
            fprintf(fh,"%u %u\n", vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono);
         for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++,offset+=2){
//            fprintf(fh,"%u %u\n", vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono);
            int_vector[offset]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono;
            int_vector[offset+1]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono;
//      fprintf(stderr,"was here\n");
         
         }
      }
   }
   for(i=0;i<vesicle->tlist->n;i++){
      fprintf(fh,"%u %u %u\n", vesicle->tlist->tria[i]->vertex[0]->idx, vesicle->tlist->tria[i]->vertex[1]->idx, vesicle->tlist->tria[i]->vertex[2]->idx);
   for(i=0;i<vesicle->tlist->n;i++,offset+=3){
//      fprintf(fh,"%u %u %u\n", vesicle->tlist->tria[i]->vertex[0]->idx, vesicle->tlist->tria[i]->vertex[1]->idx, vesicle->tlist->tria[i]->vertex[2]->idx);
      int_vector[offset]=vesicle->tlist->tria[i]->vertex[0]->idx;
      int_vector[offset+1]=vesicle->tlist->tria[i]->vertex[1]->idx;
      int_vector[offset+2]=vesicle->tlist->tria[i]->vertex[2]->idx;
   }
    fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
    for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2){
    fprintf(fh,"%u ",i);
      printout=ts_compress((char *)int_vector,offset*sizeof(long));
      fprintf(fh,"%s",printout);
      free(printout);
    fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"binary\">");
    offset=0;
    for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2,offset++){
//    fprintf(fh,"%u ",i);
   int_vector[offset]=i;
    }
   for(j=i+1;j<i+3*(vesicle->tlist->n);j+=3){ //let's continue counting from where we left of
      fprintf(fh,"%u ", j);
   for(j=i+1;j<i+3*(vesicle->tlist->n);j+=3,offset++){ //let's continue counting from where we left of
//      fprintf(fh,"%u ", j);
   int_vector[offset]=j;
   }
    fprintf(fh,"\n");
    fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
     for (i=0;i<blist->n+monono*polyno+(fonono-1)*filno;i++){
        fprintf(fh,"3 ");
//    fprintf(fh,"\n");
   printout=ts_compress((char *)int_vector,offset*sizeof(long));
   fprintf(fh,"%s",printout);
   free(printout);
    fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"binary\">");
   free(int_vector);
   free(float_vector);
   free(float_cell_vector);
   free(coordinate_vector);
   unsigned char *char_vector=(unsigned char *)malloc((blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n)*sizeof(char));
   offset=0;
     for (i=0;i<blist->n+monono*polyno+(fonono-1)*filno;i++,offset++){
        //fprintf(fh,"3 ");
   char_vector[offset]=3;
    }
   for(i=0;i<vesicle->tlist->n;i++){
      fprintf(fh,"5 ");
   for(i=0;i<vesicle->tlist->n;i++,offset++){
   //   fprintf(fh,"5 ");
      char_vector[offset]=5;
   }
   printout=ts_compress((char *)char_vector,offset*sizeof(char));
   fprintf(fh,"%s",printout);
   free(printout);
   free(char_vector);
    fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
    fclose(fh);
    return TS_SUCCESS;