Trisurf Monte Carlo simulator
Samo Penic
2014-03-11 3f5c83bc26e98dd41e93fc9023a610cd0ec1d363
commit | author | age
f74313 1 #include "general.h"
d7639a 2 #include<stdio.h>
SP 3 #include "io.h"
f74313 4 #include <confuse.h>
SP 5 #include "vertex.h"
6 #include "bond.h"
d7639a 7 #include<string.h>
a6b1b5 8 #include<stdlib.h>
d7639a 9 #include <sys/types.h>
SP 10 #include <dirent.h>
b14a8d 11 #include "initial_distribution.h"
a2db52 12 #include "poly.h"
d7639a 13
083e03 14 #include <getopt.h>
SP 15 #include <sys/stat.h>
16 #include <sys/types.h>
17 #include <dirent.h>
18 #include <errno.h>
e9eab4 19 /** DUMP STATE TO DISK DRIVE **/
SP 20
1ab449 21 ts_bool dump_state(ts_vesicle *vesicle, ts_uint iteration){
e9eab4 22
SP 23     /* save current state with wrong pointers. Will fix that later */
24     ts_uint i,j,k;
25     FILE *fh=fopen("dump.bin","wb");
26
27     /* dump vesicle */
86f5e7 28     fwrite(vesicle, sizeof(ts_vesicle),1,fh);
e9eab4 29     /* dump vertex list */
SP 30     fwrite(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
31     /* dump bond list */
32     fwrite(vesicle->blist, sizeof(ts_bond_list),1,fh);
33     /* dump triangle list */
34     fwrite(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
35     /* dump cell list */
36     fwrite(vesicle->clist, sizeof(ts_cell_list),1,fh);
37     /* dump poly list */
38     fwrite(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
39     /* level 1 complete */
40
41     /*dump vertices*/
42     for(i=0;i<vesicle->vlist->n;i++){
43         fwrite(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
44         /* dump pointer offsets for:
45                     neigh
46                     bond
47                     tria
48                     cell is ignored
49         */
50         for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
3c1ac1 51             fwrite(&vesicle->vlist->vtx[i]->neigh[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 52         }
SP 53         for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
3c1ac1 54             fwrite(&vesicle->vlist->vtx[i]->bond[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 55         }
SP 56         for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
3c1ac1 57             fwrite(&vesicle->vlist->vtx[i]->tristar[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 58         }
SP 59     }
60
61     /*dump bonds*/
62     for(i=0;i<vesicle->blist->n;i++){
63         fwrite(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
64         /* dump pointer offsets for vtx1 and vtx2 */
3c1ac1 65         //off=(ts_ulong)(vesicle->blist->bond[i]->vtx1-vesicle->vlist->vtx[0]);
SP 66         fwrite(&vesicle->blist->bond[i]->vtx1->idx,sizeof(ts_uint),1,fh); 
67         //off=(ts_ulong)(vesicle->blist->bond[i]->vtx2-vesicle->vlist->vtx[0]);
68         fwrite(&vesicle->blist->bond[i]->vtx2->idx,sizeof(ts_uint),1,fh); 
e9eab4 69     }
SP 70
71     /*dump triangles*/
72     for(i=0;i<vesicle->tlist->n;i++){
73         fwrite(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
74         /* dump pointer offsets for vertex */
3c1ac1 75         fwrite(&vesicle->tlist->tria[i]->vertex[0]->idx,sizeof(ts_uint),1,fh); 
SP 76         fwrite(&vesicle->tlist->tria[i]->vertex[1]->idx,sizeof(ts_uint),1,fh); 
77         fwrite(&vesicle->tlist->tria[i]->vertex[2]->idx,sizeof(ts_uint),1,fh); 
e9eab4 78         /* dump pointer offsets for neigh */
SP 79         for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
3c1ac1 80             fwrite(&vesicle->tlist->tria[i]->neigh[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 81         }
SP 82     }
83
84
85     /*dump polymeres */
86     for(i=0;i<vesicle->poly_list->n;i++){
87         fwrite(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
3c1ac1 88         fwrite(vesicle->poly_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
SP 89         fwrite(vesicle->poly_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
e9eab4 90     } 
SP 91      
92     /* dump poly vertex(monomer) list*/
93     for(i=0;i<vesicle->poly_list->n;i++){
94         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
95             fwrite(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
96             /* dump offset for neigh and bond */
97             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
3c1ac1 98                // off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]-vesicle->poly_list->poly[i]->vlist->vtx[0]);
SP 99                 fwrite(&vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]->idx,sizeof(ts_uint),1,fh); 
e9eab4 100             }
SP 101             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
3c1ac1 102                 //off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]-vesicle->poly_list->poly[i]->blist->bond[0]);
SP 103                 fwrite(&vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]->idx,sizeof(ts_uint),1,fh); 
e9eab4 104             }
SP 105         }
3c1ac1 106     // grafted vtx on vesicle data dump
SP 107         fwrite(&vesicle->poly_list->poly[i]->grafted_vtx->idx, sizeof(ts_uint),1,fh);
e9eab4 108     }
SP 109     /* dump poly bonds between monomers list*/
110     for(i=0;i<vesicle->poly_list->n;i++){
111         for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
112             fwrite(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
113             /* dump vtx1 and vtx2 offsets */
3c1ac1 114             //off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx1-vesicle->poly_list->poly[i]->vlist->vtx[0]);
SP 115             fwrite(&vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx,sizeof(ts_uint),1,fh); 
116 //            off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]);
117             fwrite(&vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx,sizeof(ts_uint),1,fh); 
e9eab4 118         }
SP 119     }
120
121 /* pointer offsets for fixing the restored pointers */
122 /* need pointers for 
123     vlist->vtx
124     blist->bond
125     tlist->tria
126     clist->cell
127     poly_list->poly
128     and for each poly:
129         poly_list->poly->vtx
130         poly_list->poly->bond
131 */
132
063289 133     fwrite(vesicle->clist, sizeof(ts_cell_list),1,  fh);
1ab449 134     
SP 135     fwrite(&iteration, sizeof(ts_uint),1,fh);
e9eab4 136     fclose(fh);
SP 137     return TS_SUCCESS;
138 }
139
140
141 /** RESTORE DUMP FROM DISK **/
1ab449 142 ts_vesicle *restore_state(ts_uint *iteration){
e9eab4 143     ts_uint i,j,k;
SP 144     FILE *fh=fopen("dump.bin","rb");
3f5c83 145
SP 146     struct stat sb;
147     if (stat("dump.bin", &sb) == -1) {
148         //dump file does not exist.
149         return NULL;
150     }
151
152     //check if it is regular file
153     if((sb.st_mode & S_IFMT) != S_IFREG) {
154         //dump file is not a regular file.
155         ts_fprintf(stderr,"Dump file is not a regular file!\n");
156         return NULL;
157     }
158
e9eab4 159     ts_uint retval;
3c1ac1 160     ts_uint idx;
e9eab4 161
SP 162 /* we restore all the data from the dump */
163     /* restore vesicle */
164     ts_vesicle *vesicle=(ts_vesicle *)calloc(1,sizeof(ts_vesicle));
86f5e7 165     retval=fread(vesicle, sizeof(ts_vesicle),1,fh);
SP 166 //    fprintf(stderr,"was here! %e\n",vesicle->dmax);
167
e9eab4 168     /* restore vertex list */
SP 169     vesicle->vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
170     retval=fread(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
171     /* restore bond list */
172     vesicle->blist=(ts_bond_list *)malloc(sizeof(ts_bond_list));
173     retval=fread(vesicle->blist, sizeof(ts_bond_list),1,fh);
174     /* restore triangle list */
175     vesicle->tlist=(ts_triangle_list *)malloc(sizeof(ts_triangle_list));
176     retval=fread(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
177     /* restore cell list */
178     vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
179     retval=fread(vesicle->clist, sizeof(ts_cell_list),1,fh);
180     /* restore poly list */
181     vesicle->poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
182     retval=fread(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
183     /* level 1 complete */
184
3c1ac1 185 /* prerequisity. Bonds must be malloced before vertexes are recreated */
SP 186   vesicle->blist->bond=(ts_bond **)calloc(vesicle->blist->n,sizeof(ts_bond *));
187     for(i=0;i<vesicle->blist->n;i++){
188         vesicle->blist->bond[i]=(ts_bond *)malloc(sizeof(ts_bond));
189     }
190 /* prerequisity. Triangles must be malloced before vertexes are recreated */
191   vesicle->tlist->tria=(ts_triangle **)calloc(vesicle->tlist->n,sizeof(ts_triangle *));
192     for(i=0;i<vesicle->tlist->n;i++){
193         vesicle->tlist->tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle));
194 }
a37ba2 195 /* prerequisity. Vertices must be malloced before vertexes are recreated */
e9eab4 196     vesicle->vlist->vtx=(ts_vertex **)calloc(vesicle->vlist->n,sizeof(ts_vertex *));
a37ba2 197  for(i=0;i<vesicle->vlist->n;i++){
e9eab4 198         vesicle->vlist->vtx[i]=(ts_vertex *)malloc(sizeof(ts_vertex));
a37ba2 199  }
SP 200   /*restore vertices*/
201     for(i=0;i<vesicle->vlist->n;i++){
e9eab4 202         retval=fread(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
SP 203         /*restore neigh, bond, tristar. Ignoring cell */
204         vesicle->vlist->vtx[i]->neigh=(ts_vertex **)calloc(vesicle->vlist->vtx[i]->neigh_no, sizeof(ts_vertex *));
205         for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
3c1ac1 206             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 207             vesicle->vlist->vtx[i]->neigh[j]=vesicle->vlist->vtx[idx];
e9eab4 208         }
SP 209         vesicle->vlist->vtx[i]->bond=(ts_bond **)calloc(vesicle->vlist->vtx[i]->bond_no, sizeof(ts_bond *));
210         for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
3c1ac1 211             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 212 /* pointer can be assigned only when list of bonds is fully initialized in memory. Thus bondlist popularization must be done before vertex can reference to it */
213             vesicle->vlist->vtx[i]->bond[j]=vesicle->blist->bond[idx];    
e9eab4 214         }
SP 215
216         vesicle->vlist->vtx[i]->tristar=(ts_triangle **)calloc(vesicle->vlist->vtx[i]->tristar_no, sizeof(ts_triangle *));
217         for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
3c1ac1 218             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 219 /* same comment as above */
220             vesicle->vlist->vtx[i]->tristar[j]=vesicle->tlist->tria[idx];
e9eab4 221         }
SP 222
223     }
224
225     /*restore bonds*/
3c1ac1 226    // vesicle->blist->bond=(ts_bond **)calloc(vesicle->blist->n,sizeof(ts_bond *)); // done before.
e9eab4 227     for(i=0;i<vesicle->blist->n;i++){
3c1ac1 228      //   vesicle->blist->bond[i]=(ts_bond *)malloc(sizeof(ts_bond)); //done before.
e9eab4 229         retval=fread(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
SP 230         /* restore vtx1 and vtx2 */
3c1ac1 231         retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 232         vesicle->blist->bond[i]->vtx1=vesicle->vlist->vtx[idx];
233         retval=fread(&idx,sizeof(ts_uint),1,fh);
234         vesicle->blist->bond[i]->vtx2=vesicle->vlist->vtx[idx];
e9eab4 235     }
SP 236
237     /*restore triangles*/
3c1ac1 238 //    vesicle->tlist->tria=(ts_triangle **)calloc(vesicle->tlist->n,sizeof(ts_triangle *)); // done before
e9eab4 239     for(i=0;i<vesicle->tlist->n;i++){
3c1ac1 240  //       vesicle->tlist->tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle)); // done before
e9eab4 241         retval=fread(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
SP 242         /* restore pointers for vertices */
3c1ac1 243         retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 244         vesicle->tlist->tria[i]->vertex[0]=vesicle->vlist->vtx[idx];
245         retval=fread(&idx,sizeof(ts_uint),1,fh);
246         vesicle->tlist->tria[i]->vertex[1]=vesicle->vlist->vtx[idx];
247         retval=fread(&idx,sizeof(ts_uint),1,fh);
248         vesicle->tlist->tria[i]->vertex[2]=vesicle->vlist->vtx[idx];
e9eab4 249         /* restore pointers for neigh */
3c1ac1 250      vesicle->tlist->tria[i]->neigh=(ts_triangle **)malloc(vesicle->tlist->tria[i]->neigh_no*sizeof(ts_triangle *));
e9eab4 251         for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
3c1ac1 252             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 253             vesicle->tlist->tria[i]->neigh[j]=vesicle->tlist->tria[idx];
e9eab4 254         }
SP 255
256     }
257    
258     /*restore cells */
259 /*TODO: do we need to recalculate cells here? */
260 /*    vesicle->clist->cell=(ts_cell **)malloc(vesicle->clist->cellno*sizeof(ts_cell *));
261     for(i=0;i<vesicle->clist->cellno;i++){
262         vesicle->clist->cell[i]=(ts_cell *)malloc(sizeof(ts_cell));
263         retval=fread(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
264     }
265 */
266     /*restore polymeres */
267     vesicle->poly_list->poly = (ts_poly **)calloc(vesicle->poly_list->n,sizeof(ts_poly *));
268     for(i=0;i<vesicle->poly_list->n;i++){
269         vesicle->poly_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly));
270         retval=fread(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
3c1ac1 271         vesicle->poly_list->poly[i]->vlist=(ts_vertex_list *)calloc(1,sizeof(ts_vertex_list));
SP 272         retval=fread(vesicle->poly_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
273         vesicle->poly_list->poly[i]->blist=(ts_bond_list *)calloc(1,sizeof(ts_bond_list));
274         retval=fread(vesicle->poly_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
275     /* initialize adress space for pointers that will hold specific vertices (monomers) and bonds */
276         vesicle->poly_list->poly[i]->vlist->vtx=(ts_vertex **)calloc(vesicle->poly_list->poly[i]->vlist->n,sizeof(ts_vertex *));
277         vesicle->poly_list->poly[i]->blist->bond=(ts_bond **)calloc(vesicle->poly_list->poly[i]->blist->n,sizeof(ts_bond *));
278      for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
279             vesicle->poly_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex));
280     }
281     for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
282             vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
283     }
284
e9eab4 285     } 
3c1ac1 286
e9eab4 287      
SP 288     /* restore poly vertex(monomer) list*/
289     for(i=0;i<vesicle->poly_list->n;i++){
290         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
291             retval=fread(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
3c1ac1 292                 
e9eab4 293             /* restore neigh and bonds */
SP 294             vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh=(ts_vertex **)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no, sizeof(ts_vertex *));
295             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
3c1ac1 296                 retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 297                 vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->poly_list->poly[i]->vlist->vtx[idx];
e9eab4 298             }
SP 299             vesicle->poly_list->poly[i]->vlist->vtx[j]->bond=(ts_bond **)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no, sizeof(ts_bond *));
300             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
3c1ac1 301                 retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 302                 vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->poly_list->poly[i]->blist->bond[idx];
e9eab4 303             }
SP 304
305         }
3c1ac1 306     /* restore grafted vtx on vesicle and grafted_poly */
SP 307                 retval=fread(&idx,sizeof(ts_uint),1,fh);
308         vesicle->vlist->vtx[idx]->grafted_poly=vesicle->poly_list->poly[i];
309         vesicle->poly_list->poly[i]->grafted_vtx=vesicle->vlist->vtx[idx];    
e9eab4 310     }
3c1ac1 311
e9eab4 312     /* restore poly bonds between monomers list*/
SP 313     for(i=0;i<vesicle->poly_list->n;i++){
314         for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
3c1ac1 315        //     vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
e9eab4 316             retval=fread(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
SP 317             /* restore vtx1 and vtx2 */
3c1ac1 318                 retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 319                 vesicle->poly_list->poly[i]->blist->bond[j]->vtx1=vesicle->poly_list->poly[i]->vlist->vtx[idx];
320                 retval=fread(&idx,sizeof(ts_uint),1,fh);
321                 vesicle->poly_list->poly[i]->blist->bond[j]->vtx2=vesicle->poly_list->poly[i]->vlist->vtx[idx];
e9eab4 322         }
SP 323     }
324
063289 325 // recreating space for cells // 
SP 326     vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
327     retval=fread(vesicle->clist, sizeof(ts_cell_list), 1,fh); 
328     vesicle->clist->cell=(ts_cell **)malloc(sizeof(ts_cell *)*vesicle->clist->ncmax[0]*vesicle->clist->ncmax[1]*vesicle->clist->ncmax[2]);
329     for(i=0;i<vesicle->clist->ncmax[0]*vesicle->clist->ncmax[1]*vesicle->clist->ncmax[2];i++){
330             vesicle->clist->cell[i]=(ts_cell *)calloc(1,sizeof(ts_cell));
331             vesicle->clist->cell[i]->idx=i+1; // We enumerate cells! Probably never required!
332         }
e9eab4 333
1ab449 334     retval=fread(iteration,sizeof(ts_uint),1,fh);
e9eab4 335     if(retval); 
SP 336     fclose(fh);
337     return vesicle;
338 }
339
340
341
083e03 342 ts_bool parse_args(int argc, char **argv){
3f5c83 343     int c, retval;
SP 344     struct stat sb;
345     sprintf(command_line_args.path, "./"); //clear string;
346     sprintf(command_line_args.output_fullfilename,"./output.pvd");
347     sprintf(command_line_args.dump_fullfilename,"./dump.bin");
348     sprintf(command_line_args.tape_fullfilename,"./tape");
349             FILE *file;
350     
083e03 351 while (1)
SP 352      {
353        static struct option long_options[] =
354          {
8a6614 355            {"force-from-tape", no_argument,       &(command_line_args.force_from_tape), 1},
SP 356        {"reset-iteration-count", no_argument, &(command_line_args.reset_iteration_count), 1},
083e03 357            {"tape",     no_argument,       0, 't'},
SP 358            {"output-file",  required_argument, 0, 'o'},
359            {"directory",  required_argument, 0, 'd'},
3f5c83 360            {"dump-filename", required_argument,0, 'f'},
083e03 361            {0, 0, 0, 0}
SP 362          };
363        /* getopt_long stores the option index here. */
364        int option_index = 0;
365
3f5c83 366        c = getopt_long (argc, argv, "d:f:o:t:",
083e03 367                         long_options, &option_index);
SP 368
369        /* Detect the end of the options. */
370        if (c == -1)
371          break;
372
373        switch (c)
374          {
375          case 0:
376            /* If this option set a flag, do nothing else now. */
377            if (long_options[option_index].flag != 0)
378              break;
379            printf ("option %s", long_options[option_index].name);
380            if (optarg)
381              printf (" with arg %s", optarg);
382            printf ("\n");
383            break;
384
385          case 't':
386             //check if tape exists. If not, fail immediately.
3f5c83 387             if (stat(optarg, &sb) == -1) {
SP 388                 ts_fprintf(stderr,"Tape '%s' does not exist!\n",optarg);
389                 fatal("Please select correct tape",1);
390             } else {
391                 strcpy(command_line_args.tape_fullfilename,optarg);
392             }
083e03 393            break;
SP 394
395          case 'o':
396             //set filename of master output file
3f5c83 397             if ((file = fopen(optarg, "w")) == NULL) {
SP 398                 fprintf(stderr,"Could not create output file!\n");
399                 fatal("Please specify correct output file",1);
400                 
401             } else {
402                 fclose(file);
403             }
404             strcpy(command_line_args.output_fullfilename, optarg);
405             break;
083e03 406
SP 407          case 'd':
408             //check if directory exists. If not create one. If creation is
409             //successful, set directory for output files.
410             //printf ("option -d with value `%s'\n", optarg);
3f5c83 411             if (stat(optarg, &sb) == -1) {
SP 412                 //directory does not exist
083e03 413                 retval=mkdir(optarg, 0700);
SP 414                 if(retval){
3f5c83 415                     fatal("Could not create requested directory. Check if you have permissions",1);
083e03 416                 }
SP 417             }
3f5c83 418             //check if is a proper directory
SP 419             else if((sb.st_mode & S_IFMT) != S_IFDIR) {
420                 //it is not a directory. fatal error.
421                 ts_fprintf(stderr,"%s is not a directory!\n",optarg);
422                 fatal("Cannot continue",1);
083e03 423             }
3f5c83 424             strcpy(command_line_args.path, optarg);
083e03 425            break;
SP 426
427         case 'f':
428             //check if dump file specified exists. Defaults to dump.bin
3f5c83 429             if ((file = fopen(optarg, "w")) == NULL) {
SP 430                 fprintf(stderr,"Could not create dump file!\n");
431                 fatal("Please specify correct dump file",1);
432                 
433             } else {
434                 fclose(file);
435             }
436             strcpy(command_line_args.dump_fullfilename, optarg);
083e03 437             break;
SP 438
439          case '?':
440            /* getopt_long already printed an error message. */
441
442             ts_fprintf(stderr,"\n\nhere comes the help.\n\n");
443             fatal("Ooops, read help first",1);
444            break;
445
446          default:
447            exit (1);
448          }
449      }
450
451     return TS_SUCCESS;
452
453 }
454
455
456
457
d7639a 458 ts_bool print_vertex_list(ts_vertex_list *vlist){
SP 459     ts_uint i;
460     printf("Number of vertices: %u\n",vlist->n);
461     for(i=0;i<vlist->n;i++){
a6b1b5 462         printf("%u: %f %f %f\n",
8f6a69 463 vlist->vtx[i]->idx,vlist->vtx[i]->x, vlist->vtx[i]->y, vlist->vtx[i]->z);
d7639a 464     }
SP 465     return TS_SUCCESS;
466 }
467
468 ts_bool print_vertex_neighbours(ts_vertex_list *vlist){
469     ts_uint i,j;
a6b1b5 470     ts_vertex **vtx=vlist->vtx;
d7639a 471     printf("Vertex id(neigh no): (neighvertex coord) (neighvertex coord) ...\n");
SP 472     for(i=0;i<vlist->n;i++){
8f6a69 473         printf("%u(%u): ",vtx[i]->idx,vtx[i]->neigh_no);
SP 474         for(j=0;j<vtx[i]->neigh_no;j++){
475             printf("(%f,%f,%f)",vtx[i]->neigh[j]->x,
476 vtx[i]->neigh[j]->y,vtx[i]->neigh[j]->z);
d7639a 477         }
SP 478         printf("\n");
479     }
480
481 return TS_SUCCESS;
482 }
483
484 ts_bool write_vertex_fcompat_file(ts_vertex_list *vlist,ts_char *filename){
a6b1b5 485     ts_vertex **vtx=vlist->vtx;
d7639a 486     ts_uint i;
SP 487     FILE *fh;
488     
489     fh=fopen(filename, "w");
490     if(fh==NULL){
491         err("Cannot open file %s for writing");
492         return TS_FAIL;
493     }
494     for(i=0;i<vlist->n;i++)
8f6a69 495         fprintf(fh," %E\t%E\t%E\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a 496
SP 497     fclose(fh);
498 return TS_SUCCESS;
499 }
500
501
502 ts_bool fprint_vertex_list(FILE *fh,ts_vertex_list *vlist){
503     ts_uint i,j;
504     for(i=0;i<vlist->n;i++){
8f6a69 505         fprintf(fh," %.17E\t%.17E\t%.17E\t%u\n",vlist->vtx[i]->x,
SP 506             vlist->vtx[i]->y, vlist->vtx[i]->z,
507             vlist->vtx[i]->neigh_no);
508         for(j=0;j<vlist->vtx[i]->neigh_no;j++){
509             fprintf(fh,"\t%u",(ts_uint)(vlist->vtx[i]->neigh[j]->idx));
a6b1b5 510         //-vlist->vtx+1));
d7639a 511         }
SP 512         fprintf(fh,"\n");
513     }
514     return TS_SUCCESS;
515 }
516
517 ts_bool fprint_tristar(FILE *fh, ts_vesicle *vesicle){
518     ts_uint i,j;
a6b1b5 519     for(i=0;i<vesicle->vlist->n;i++){
8f6a69 520         fprintf(fh,"\t%u",vesicle->vlist->vtx[i]->tristar_no);
SP 521         for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
522             fprintf(fh,"\t%u",(ts_uint)(vesicle->vlist->vtx[i]->tristar[j]->idx));//-vesicle->tlist->tria+1));
d7639a 523         }
SP 524         fprintf(fh,"\n");
525     }
526     return TS_SUCCESS;
527 }
528
529 ts_bool fprint_triangle_list(FILE *fh, ts_vesicle *vesicle){
a6b1b5 530         ts_triangle_list *tlist=vesicle->tlist;
d7639a 531       ts_uint i,j;
SP 532     for(i=0;i<tlist->n;i++){
41a035 533         fprintf(fh,"\t%u",tlist->tria[i]->neigh_no);
SP 534         for(j=0;j<tlist->tria[i]->neigh_no;j++){
535             fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->neigh[j]->idx));//-tlist->tria+1)); 
d7639a 536         }
SP 537         fprintf(fh,"\n");
538             for(j=0;j<3;j++){
41a035 539             fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->vertex[j]->idx));//-vesicle->vlist->vtx+1)); 
d7639a 540             }
SP 541         fprintf(fh,"\n");
41a035 542         fprintf(fh,"%.17E\t%.17E\t%.17E\n",tlist->tria[i]->xnorm,
SP 543 tlist->tria[i]->ynorm,tlist->tria[i]->znorm);
d7639a 544         fprintf(fh,"0.00000000000000000\n0.00000000000000000\n");
SP 545     }
546     return TS_SUCCESS;
547 }
548
549 ts_bool fprint_vertex_data(FILE *fh,ts_vertex_list *vlist){
550     ts_uint i,j;
551     for(i=0;i<vlist->n;i++){
552         fprintf(fh," %.17E\t%.17E\t%.17E\t%.17E\t%.17E\t%u\n",
8f6a69 553         vlist->vtx[i]->xk,vlist->vtx[i]->c,vlist->vtx[i]->energy,
SP 554         vlist->vtx[i]->energy_h, vlist->vtx[i]->curvature, 0);
cef6ca 555         for(j=0;j<vlist->vtx[i]->bond_no;j++){
8f6a69 556             fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length_dual);
d7639a 557         }
SP 558             fprintf(fh,"\n");
cef6ca 559         for(j=0;j<vlist->vtx[i]->bond_no;j++){
8f6a69 560             fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length);
d7639a 561         }
SP 562             fprintf(fh,"\n");
563     }
564     return TS_SUCCESS;
565 }
566
567 ts_bool fprint_bonds(FILE *fh,ts_vesicle *vesicle){
568     ts_uint i;
a6b1b5 569     for(i=0;i<vesicle->blist->n;i++){
e016c4 570         fprintf(fh,"\t%u\t%u\n",(ts_uint)(vesicle->blist->bond[i]->vtx1->idx),
a6b1b5 571 //-vesicle->vlist->vtx+1),
e016c4 572         (ts_uint)(vesicle->blist->bond[i]->vtx2->idx));
a6b1b5 573     //-vesicle->vlist.vtx+1));
d7639a 574     }
SP 575     return TS_SUCCESS;
576 }
577
578
579 ts_bool write_dout_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
580     FILE *fh;
581     fh=fopen(filename, "w");
582     if(fh==NULL){
583         err("Cannot open file %s for writing");
584         return TS_FAIL;
585     }
586     fprintf(fh,"%.17E\n%.17E\n",vesicle->stepsize,vesicle->dmax);
a6b1b5 587     fprint_vertex_list(fh,vesicle->vlist);
d7639a 588     fprint_tristar(fh,vesicle);
SP 589     fprint_triangle_list(fh,vesicle);
a6b1b5 590     fprint_vertex_data(fh,vesicle->vlist);
d7639a 591     fprint_bonds(fh,vesicle);
SP 592     fclose(fh);    
593     return TS_SUCCESS;
594 }
595
596 ts_bool read_tape_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
597     FILE *fh;
598     char line[255];
599     fh=fopen(filename, "r");
600         if(fh==NULL){
601                 err("Cannot open file for reading... Nonexistant file?");
602                 return TS_FAIL;
603         }
a6b1b5 604     ts_uint retval=1;
d7639a 605     while(retval!=EOF){
SP 606         retval=fscanf(fh,"%s",line);
607         
608         fprintf(stderr,"%s",line);
609     }    
610     fclose(fh);    
611     return TS_SUCCESS;
612 }
613
614 ts_bool write_master_xml_file(ts_char *filename){
615      FILE *fh;
616     ts_char *i,*j;
617     ts_uint tstep;
618     ts_char *number;
619         fh=fopen(filename, "w");
620         if(fh==NULL){
621                 err("Cannot open file %s for writing");
622                 return TS_FAIL;
623         }
624
625     fprintf(fh,"<?xml version=\"1.0\"?>\n<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n<Collection>");
626     DIR *dir = opendir(".");
627     if(dir){
628         struct dirent *ent;
629         tstep=0;
630         while((ent = readdir(dir)) != NULL)
631         {
632             i=rindex(ent->d_name,'.');
633             if(i==NULL) continue;
634             if(strcmp(i+1,"vtu")==0){
635                     j=rindex(ent->d_name,'_');
636                     if(j==NULL) continue;
637                     number=strndup(j+1,j-i); 
a6b1b5 638                 fprintf(fh,"<DataSet timestep=\"%u\" group=\"\" part=\"0\" file=\"%s\"/>\n",atoi(number),ent->d_name);
d7639a 639                 tstep++;
SP 640                     free(number);
641             }  
642         }
643     }
f74313 644     free(dir);
d7639a 645     fprintf(fh,"</Collection>\n</VTKFile>\n");
SP 646     fclose(fh);
647     return TS_SUCCESS;
648 }
649
650 ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno){
a6b1b5 651     ts_vertex_list *vlist=vesicle->vlist;
SP 652     ts_bond_list *blist=vesicle->blist;
653     ts_vertex **vtx=vlist->vtx;
40aa5b 654     ts_uint i,j;
d7639a 655         char filename[255];
SP 656     FILE *fh;
657
658         sprintf(filename,"timestep_%.6u.vtu",timestepno);
659     fh=fopen(filename, "w");
660     if(fh==NULL){
661         err("Cannot open file %s for writing");
662         return TS_FAIL;
663     }
664     /* Here comes header of the file */
40aa5b 665
SP 666     //find number of extra vtxs and bonds of polymeres
3c1ac1 667     ts_uint monono=0, polyno=0, poly_idx=0;
40aa5b 668     ts_bool poly=0;
SP 669     if(vesicle->poly_list!=NULL){
670         if(vesicle->poly_list->poly[0]!=NULL){
671         polyno=vesicle->poly_list->n;
672         monono=vesicle->poly_list->poly[0]->vlist->n;
673         poly=1;
674         }
675     }
d7639a 676     fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n <UnstructuredGrid>\n");
40aa5b 677     fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno, blist->n+monono*polyno);
d7639a 678     fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
SP 679        for(i=0;i<vlist->n;i++){
a6b1b5 680         fprintf(fh,"%u ",vtx[i]->idx);
d7639a 681     }
40aa5b 682     //polymeres
SP 683     if(poly){
3c1ac1 684         poly_idx=vlist->n;
40aa5b 685         for(i=0;i<vesicle->poly_list->n;i++){
3c1ac1 686             for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
SP 687                 fprintf(fh,"%u ", poly_idx);
40aa5b 688             }
SP 689         }
690     }
d7639a 691
SP 692     fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
693     for(i=0;i<vlist->n;i++){
8f6a69 694         fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
40aa5b 695     }
SP 696     //polymeres
697     if(poly){
698         for(i=0;i<vesicle->poly_list->n;i++){
699             for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
700                 fprintf(fh,"%e %e %e\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 );
701             }
702         }
d7639a 703     }
SP 704
705     fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
706     for(i=0;i<blist->n;i++){
e016c4 707             fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a 708     }
40aa5b 709     //polymeres
SP 710     if(poly){
3c1ac1 711         poly_idx=vlist->n;
40aa5b 712         for(i=0;i<vesicle->poly_list->n;i++){
SP 713             for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
3c1ac1 714 //                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);
SP 715                 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);
40aa5b 716             }
SP 717     //grafted bonds
3c1ac1 718         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);
40aa5b 719         }
SP 720
721     }
722     
723
d7639a 724     fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
40aa5b 725     for (i=2;i<(blist->n+monono*polyno)*2+1;i+=2){
d7639a 726     fprintf(fh,"%u ",i);
SP 727     }
728     fprintf(fh,"\n");
729     fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
40aa5b 730      for (i=0;i<blist->n+monono*polyno;i++){
d7639a 731         fprintf(fh,"3 ");
SP 732     }
733
734     fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
735     fclose(fh);
736     return TS_SUCCESS;
737
738 }
739
740 ts_bool write_vertex_vtk_file(ts_vesicle *vesicle,ts_char *filename, ts_char *text){
a6b1b5 741     ts_vertex_list *vlist=vesicle->vlist;
SP 742     ts_bond_list *blist=vesicle->blist;
743     ts_vertex **vtx=vlist->vtx;
744     ts_uint i;
d7639a 745     FILE *fh;
SP 746     
747     fh=fopen(filename, "w");
748     if(fh==NULL){
749         err("Cannot open file %s for writing");
750         return TS_FAIL;
751     }
752     /* Here comes header of the file */
753 //    fprintf(stderr,"NSHELL=%u\n",nshell);
754     fprintf(fh, "# vtk DataFile Version 2.0\n");
755     /* TODO: Do a sanity check on text. Max 255 char, must not me \n terminated */ 
756     fprintf(fh, "%s\n", text);
757     fprintf(fh,"ASCII\n");
758     fprintf(fh,"DATASET UNSTRUCTURED_GRID\n");
759     fprintf(fh,"POINTS %u double\n", vlist->n);
760     for(i=0;i<vlist->n;i++){
8f6a69 761         fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a 762     }
SP 763     
764     fprintf(fh,"CELLS %u %u\n",blist->n,3*blist->n);
765     for(i=0;i<blist->n;i++){
e016c4 766             fprintf(fh,"2 %u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a 767     }
SP 768     fprintf(fh,"CELL_TYPES %u\n",blist->n);
769     for(i=0;i<blist->n;i++)
770         fprintf(fh,"3\n");
771
772     fprintf(fh,"POINT_DATA %u\n", vlist->n);
773     fprintf(fh,"SCALARS scalars long 1\n");
774     fprintf(fh,"LOOKUP_TABLE default\n");
775
776     for(i=0;i<vlist->n;i++)
8f6a69 777         fprintf(fh,"%u\n",vtx[i]->idx);
d7639a 778
SP 779     fclose(fh);
780     return TS_SUCCESS;
781 }
782
783
784
1ab449 785 ts_tape *parsetape(char *filename){
SP 786   //  long int nshell=17,ncxmax=60, ncymax=60, nczmax=60, npoly=10, nmono=20, pswitch=0;  // THIS IS DUE TO CONFUSE BUG!
787     ts_tape *tape=(ts_tape *)calloc(1,sizeof(ts_tape));
788     tape->multiprocessing=calloc(255,sizeof(char));
789   /*  long int brezveze0=1;
b14a8d 790     long int brezveze1=1;
SP 791     long int brezveze2=1;
e86357 792     ts_double xk0=25.0, dmax=1.67,stepsize=0.15,kspring=800.0,pressure=0.0;
d7a113 793     long int iter=1000, init=1000, mcsw=1000;
1ab449 794 */    
40aa5b 795
d7639a 796     cfg_opt_t opts[] = {
1ab449 797         CFG_SIMPLE_INT("nshell", &tape->nshell),
SP 798         CFG_SIMPLE_INT("npoly", &tape->npoly),
799         CFG_SIMPLE_INT("nmono", &tape->nmono),
800         CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
801         CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
802     CFG_SIMPLE_INT("pswitch",&tape->pswitch),
803     CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
804     CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
805         CFG_SIMPLE_FLOAT("stepsize",&tape->stepsize),
806         CFG_SIMPLE_INT("nxmax", &tape->ncxmax),
807         CFG_SIMPLE_INT("nymax", &tape->ncymax),
808         CFG_SIMPLE_INT("nzmax", &tape->nczmax),
809         CFG_SIMPLE_INT("iterations",&tape->iterations),
810     CFG_SIMPLE_INT("mcsweeps",&tape->mcsweeps),
811     CFG_SIMPLE_INT("inititer", &tape->inititer),
812         CFG_SIMPLE_BOOL("quiet",&tape->quiet),
813         CFG_SIMPLE_STR("multiprocessing",tape->multiprocessing),
814         CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
815         CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
816         CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
d7639a 817         CFG_END()
SP 818     };
819     cfg_t *cfg;    
820     ts_uint retval;
821     cfg = cfg_init(opts, 0);
1ab449 822     retval=cfg_parse(cfg, filename);
d7639a 823     if(retval==CFG_FILE_ERROR){
a6b1b5 824     fatal("No tape file.",100);
d7639a 825     }
SP 826     else if(retval==CFG_PARSE_ERROR){
827     fatal("Invalid tape!",100);
828     }
a2db52 829
d7639a 830     cfg_free(cfg);
40aa5b 831
SP 832
1ab449 833     /* global variables are set automatically */
SP 834     quiet=tape->quiet;
835     return tape;
836 }
d7639a 837
1ab449 838 ts_bool tape_free(ts_tape *tape){
SP 839     free(tape->multiprocessing);
840     free(tape);
841     return TS_SUCCESS;
d7639a 842 }
f74313 843
SP 844
845 ts_bool read_geometry_file(char *fname, ts_vesicle *vesicle){
846     FILE *fh;
847     ts_uint i, nvtx,nedges,ntria;
848     ts_uint vtxi1,vtxi2;
849     float x,y,z;
850     ts_vertex_list *vlist;
851     fh=fopen(fname, "r");
852         if(fh==NULL){
853                 err("Cannot open file for reading... Nonexistant file?");
854                 return TS_FAIL;
855         }
856     ts_uint retval;
857     retval=fscanf(fh,"%u %u %u",&nvtx, &nedges, &ntria);
858     vesicle->vlist=init_vertex_list(nvtx);
859     vlist=vesicle->vlist;
860     for(i=0;i<nvtx;i++){
8f6a69 861    //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->x,&vlist->vtx[i]->y,&vlist->vtx[i]->z);
f74313 862        retval=fscanf(fh,"%F %F %F",&x,&y,&z);
8f6a69 863         vlist->vtx[i]->x=x;
SP 864         vlist->vtx[i]->y=y;
865         vlist->vtx[i]->z=z;
f74313 866     }
SP 867     for(i=0;i<nedges;i++){
868         retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
869         bond_add(vesicle->blist,vesicle->vlist->vtx[vtxi1-1],vesicle->vlist->vtx[vtxi2-1]);
870     }
871     //TODO: neighbours from bonds,
872     //TODO: triangles from neigbours
873
874 //    Don't need to read triangles. Already have enough data
875     /*
876     for(i=0;i<ntria;i++){
877         retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
8f6a69 878         vtxi1=vesicle->blist->vertex1->idx;
SP 879         vtxi2=vesicle->blist->vertex1->idx;
f74313 880         
SP 881     }
882     */
41a035 883     if(retval);
f74313 884     fclose(fh);    
SP 885
886
887
888     return TS_SUCCESS;
889 }