Trisurf Monte Carlo simulator
Samo Penic
2019-10-18 842378276ab4c4a8dabd02324bb24b8da39c4384
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
f74313 2 #include "general.h"
d7639a 3 #include<stdio.h>
SP 4 #include "io.h"
f74313 5 #include "vertex.h"
SP 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"
055dd7 13 #include "cell.h"
083e03 14 #include <getopt.h>
SP 15 #include <sys/stat.h>
16 #include <sys/types.h>
17 #include <dirent.h>
18 #include <errno.h>
854cb6 19 #include <snapshot.h>
e9eab4 20
SP 21
083e03 22 ts_bool parse_args(int argc, char **argv){
3f5c83 23     int c, retval;
SP 24     struct stat sb;
25     sprintf(command_line_args.path, "./"); //clear string;
7b0c07 26     sprintf(command_line_args.output_fullfilename,"output.pvd");
SP 27     sprintf(command_line_args.dump_fullfilename,"dump.bin");
28     sprintf(command_line_args.tape_fullfilename,"tape");
1bf3c3 29     sprintf(command_line_args.tape_templatefull,"./tape");
3f5c83 30             FILE *file;
SP 31     
083e03 32 while (1)
SP 33      {
34        static struct option long_options[] =
35          {
8a6614 36            {"force-from-tape", no_argument,       &(command_line_args.force_from_tape), 1},
SP 37        {"reset-iteration-count", no_argument, &(command_line_args.reset_iteration_count), 1},
083e03 38            {"tape",     no_argument,       0, 't'},
fd8126 39        {"version", no_argument, 0, 'v'},
083e03 40            {"output-file",  required_argument, 0, 'o'},
SP 41            {"directory",  required_argument, 0, 'd'},
3f5c83 42            {"dump-filename", required_argument,0, 'f'},
7b0c07 43            {"tape-options",required_argument,0,'c'},
1bf3c3 44            {"tape-template", required_argument,0,0},
ab798b 45             {"restore-from-vtk",required_argument,0,0},
083e03 46            {0, 0, 0, 0}
SP 47          };
48        /* getopt_long stores the option index here. */
49        int option_index = 0;
50
fd8126 51        c = getopt_long (argc, argv, "d:f:o:t:c:v",
083e03 52                         long_options, &option_index);
SP 53
54        /* Detect the end of the options. */
55        if (c == -1)
56          break;
57
58        switch (c)
59          {
60          case 0:
61            /* If this option set a flag, do nothing else now. */
62            if (long_options[option_index].flag != 0)
63              break;
1bf3c3 64 /*           printf ("option %s", long_options[option_index].name);
083e03 65            if (optarg)
1bf3c3 66              printf (" with arg %s", optarg); 
SP 67            printf ("\n"); */
68             //TODO: find a better way.
69             if(strcmp(long_options[option_index].name,"tape-template")==0){
70                 strcpy(command_line_args.tape_templatefull,optarg);
71             }
ab798b 72             if(strcmp(long_options[option_index].name,"restore-from-vtk")==0){
ee84bd 73                 strcpy(command_line_args.dump_from_vtk,optarg);
SP 74             }
083e03 75            break;
fd8126 76      case 'v':
SP 77         fprintf(stdout,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);
78             fprintf(stdout,"Programming done by: Samo Penic and Miha Fosnaric\n");
79             fprintf(stdout,"Released under terms of GPLv3\n");
80         exit(0);
083e03 81
7b0c07 82          case 'c':
SP 83               strcpy(command_line_args.tape_opts,optarg);
84             break;
85          case 't': //tape
3f5c83 86                 strcpy(command_line_args.tape_fullfilename,optarg);
083e03 87            break;
SP 88
7b0c07 89          case 'o':  //set filename of master pvd output file
3f5c83 90             strcpy(command_line_args.output_fullfilename, optarg);
SP 91             break;
083e03 92
SP 93          case 'd':
94             //check if directory exists. If not create one. If creation is
95             //successful, set directory for output files.
96             //printf ("option -d with value `%s'\n", optarg);
3f5c83 97             if (stat(optarg, &sb) == -1) {
SP 98                 //directory does not exist
083e03 99                 retval=mkdir(optarg, 0700);
SP 100                 if(retval){
3f5c83 101                     fatal("Could not create requested directory. Check if you have permissions",1);
083e03 102                 }
SP 103             }
3f5c83 104             //check if is a proper directory
SP 105             else if((sb.st_mode & S_IFMT) != S_IFDIR) {
106                 //it is not a directory. fatal error.
107                 ts_fprintf(stderr,"%s is not a directory!\n",optarg);
108                 fatal("Cannot continue",1);
083e03 109             }
3f5c83 110             strcpy(command_line_args.path, optarg);
083e03 111            break;
SP 112
113         case 'f':
3f5c83 114             strcpy(command_line_args.dump_fullfilename, optarg);
083e03 115             break;
SP 116
117          case '?':
118            /* getopt_long already printed an error message. */
ba73ab 119             print_help(stdout);
SP 120 //ts_fprintf(stderr,"\n\nhere comes the help.\n\n");
083e03 121             fatal("Ooops, read help first",1);
SP 122            break;
123
124          default:
125            exit (1);
126          }
127      }
128
7b0c07 129 //Here we set correct values for full filenames!
SP 130     char *buffer=(char *)malloc(10000*sizeof(char));
131     //correct the path and add trailing /
132     if(command_line_args.path[strlen(command_line_args.path)-1]!='/') strcat(command_line_args.path,"/");
133    
134 /* master pvd output file */ 
135     strcpy(buffer,command_line_args.path);
136     strcat(buffer,command_line_args.output_fullfilename);
137     if ((file = fopen(buffer, "w")) == NULL) {
138                 fprintf(stderr,"Could not create output file %s!\n", buffer);
139                 fatal("Please specify correct output file or check permissions of the file",1);
1bf3c3 140                 //there is a tape template. make a copy into desired directory
7b0c07 141                 
SP 142             } else {
143                 fclose(file);
144                 strcpy(command_line_args.output_fullfilename,buffer);
145             }
146
147 /* tape file */
148     strcpy(buffer,command_line_args.path);
149     strcat(buffer,command_line_args.tape_fullfilename);
150     if (stat(buffer, &sb) == -1) {
1bf3c3 151
SP 152                 //tape does not exist. does tape template exist?
153                 if(stat(command_line_args.tape_templatefull, &sb)==-1){ 
154                     ts_fprintf(stderr,"Tape '%s' does not exist and no tape template was specified (or does not exist)!\n",buffer);
155                     fatal("Please select correct tape or check permissions of the file",1);
156                 } else {
157                     //tape template found
158                     fatal("Samo did not program template copy yet",1); 
159                 }
7b0c07 160             } else {
SP 161                 strcpy(command_line_args.tape_fullfilename,buffer);
162             }
163
164
165 /* dump file */
166             strcpy(buffer,command_line_args.path);
167             strcat(buffer,command_line_args.dump_fullfilename);
168             //check if dump file exist first.
169             if (stat(buffer, &sb) == -1) {
170                 //no dump file. check if we can create one.
171                 if ((file = fopen(buffer, "w")) == NULL) {
172                     fprintf(stderr,"Could not create dump file '%s'!\n",buffer);
173                     fatal("Please specify correct dump file or check permissions of the file",1);
174                 } else {
175                     fclose(file);
176                     //good, file is writeable. delete it for now.
177                     remove(buffer);
178                 }
179             }
180             strcpy(command_line_args.dump_fullfilename, buffer);
181
182
183     free(buffer);
083e03 184     return TS_SUCCESS;
SP 185
186 }
187
188
ba73ab 189 ts_bool print_help(FILE *fd){
SP 190     fprintf(fd,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);
191     fprintf(fd,"Programming done by: Samo Penic and Miha Fosnaric\n");
192     fprintf(fd,"Released under terms of GPLv3\n\n");
193
194     fprintf(fd, "Invoking trisurf-ng without any flags results in default operation. Program reads 'tape' file and 'dump.bin' binary representation of the vesicle from disk and continues the simulation where it was aborted (as specified in 'dump.bin').\n\n");
195     fprintf(fd,"If 'tape' has different values than binary dump, those are used (if possible -- some values in tape, such as nshell cannot be modified).\n\n\n");
196     fprintf(fd,"However, if dump.bin is not present, user is notified to specify --force-from-tape flag. The vesicle will be created from specifications in tape.\n\n\n");
197     fprintf(fd,"Flags:\n\n");
198     fprintf(fd,"--force-from-tape\t\t makes initial shape of the vesicle from tape. Ignores already existing binary dump and possible simulation results.\n");
199     fprintf(fd,"--restore-from-vtk\t\t VTK's file ending with '.vtu' are preferred way to make state snapshots for restoration. With this flag the restoration of the vesicle from vtk is possible. The simulation will continue if hidden '.status' file with last iteration done is present. Otherwise it will start simulation from timestep 0.\n");
c4178c 200     fprintf(fd,"--reset-iteration-count\t\t starts simulation from the beginning (using binary dump).\n");
ba73ab 201     fprintf(fd,"--tape (or -t)\t\t specifies tape filename. For --force-from-tape and restoring from binary dump. Defaults to 'tape'.\n");
SP 202     fprintf(fd,"--version (or -v)\t\t Prints version information.\n");
203     fprintf(fd,"--output-file (or -o)\t\t Specifies filename of .PVD file. Defaults to 'output.pvd'\n");
c4178c 204     fprintf(fd,"--dump-filename (or -f)\t\t specifies filename for binary dump&restore. Defaults to 'dump.bin'\n\n\n");
SP 205     fprintf(fd,"Examples:\n\n");
206     fprintf(fd,"trisurf --force-from-tape\n");
207     fprintf(fd,"trisurf --reset-iteration-count\n");
208     fprintf(fd,"trisurf --restore-from-vtk filename.vtu\n");
209     fprintf(fd,"\n\n");
210
d7639a 211     return TS_SUCCESS;
SP 212 }
213
0a2c81 214 ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno, ts_cluster_list *cstlist){
a6b1b5 215     ts_vertex_list *vlist=vesicle->vlist;
SP 216     ts_bond_list *blist=vesicle->blist;
217     ts_vertex **vtx=vlist->vtx;
40aa5b 218     ts_uint i,j;
13d445 219     //ts_double senergy=0.0;
267db5 220         char filename[10000];
SP 221         char just_name[255];
d7639a 222     FILE *fh;
267db5 223         strcpy(filename,command_line_args.path);
SP 224         sprintf(just_name,"timestep_%.6u.vtu",timestepno);
225         strcat(filename,just_name);
d7639a 226
SP 227     fh=fopen(filename, "w");
228     if(fh==NULL){
229         err("Cannot open file %s for writing");
230         return TS_FAIL;
231     }
232     /* Here comes header of the file */
40aa5b 233
SP 234     //find number of extra vtxs and bonds of polymeres
58230a 235     ts_uint monono=0, polyno=0, poly_idx=0, filno=0, fonono=0;
M 236     ts_bool poly=0, fil=0;
40aa5b 237     if(vesicle->poly_list!=NULL){
SP 238         if(vesicle->poly_list->poly[0]!=NULL){
239         polyno=vesicle->poly_list->n;
240         monono=vesicle->poly_list->poly[0]->vlist->n;
241         poly=1;
242         }
243     }
58230a 244
M 245     if(vesicle->filament_list!=NULL){
246         if(vesicle->filament_list->poly[0]!=NULL){
247         filno=vesicle->filament_list->n;
248         fonono=vesicle->filament_list->poly[0]->vlist->n;
249         fil=1;
250         }
251     }
252
854cb6 253     fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n");
SP 254     xml_trisurf_data(fh,vesicle);
255     fprintf(fh, " <UnstructuredGrid>\n");
8db203 256     fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n);
80ebbe 257     fprintf(fh,"<PointData Scalars=\"vertices_idx\">\n<DataArray type=\"Int64\" Name=\"vertices_idx\" format=\"binary\">");
a69203 258     long *int_vector=(long *)malloc((vlist->n+monono*polyno+fonono*filno)*sizeof(long));
842378 259     ts_double *float_vector=(double *)malloc((vlist->n+monono*polyno+fonono*filno)*sizeof(double));
SP 260     ts_double *float_cell_vector=(double *)malloc((blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n)*sizeof(double));
261     ts_double *coordinate_vector=(double *)malloc(3*(vlist->n+monono*polyno+fonono*filno)*sizeof(double));
a69203 262 // sizes are as expected 4 bytes for int and 1 byte for char
SP 263 //    printf("size of ts_uint %ld, of int %ld, of char %ld",sizeof(ts_uint), sizeof(int), sizeof(char));
80ebbe 264     int offset=0;
d7639a 265        for(i=0;i<vlist->n;i++){
80ebbe 266     //    fprintf(fh,"%u ",vtx[i]->idx);
SP 267         int_vector[i+offset]=vtx[i]->idx;
d7639a 268     }
a69203 269     offset=offset+i;
40aa5b 270     //polymeres
SP 271     if(poly){
3c1ac1 272         poly_idx=vlist->n;
40aa5b 273         for(i=0;i<vesicle->poly_list->n;i++){
3c1ac1 274             for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
80ebbe 275                 //fprintf(fh,"%u ", poly_idx);
SP 276                 int_vector[j+offset]=poly_idx;
58230a 277             }
80ebbe 278             offset=offset+j;
58230a 279         }
M 280     }
281     //filaments
282     if(fil){
283         poly_idx=vlist->n+monono*polyno;
284         for(i=0;i<vesicle->filament_list->n;i++){
285             for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
286     //    fprintf(stderr,"was here\n");
80ebbe 287                 //fprintf(fh,"%u ", poly_idx);
SP 288                 int_vector[j+offset]=poly_idx;
40aa5b 289             }
80ebbe 290             offset=offset+j;
40aa5b 291         }
SP 292     }
a69203 293     printf("Offset in the end is %d, should be %d",offset,(vlist->n+monono*polyno+fonono*filno) );
842378 294     char *printout=ts_compress((char *)int_vector,offset*sizeof(long));
80ebbe 295     fprintf(fh,"%s",printout);
SP 296     free(printout);
d7639a 297
45c708 298         fprintf(fh,"</DataArray>\n");
0a2c81 299     if(cstlist!=NULL){
842378 300         fprintf(fh,"<DataArray type=\"Int64\" Name=\"vertices_in_cluster\" format=\"binary\">");
SP 301         offset=0;
302         for(i=0;i<vlist->n;i++,offset++){
0a2c81 303             if(vtx[i]->cluster!=NULL){
842378 304 //                fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
SP 305                 int_vector[offset]=vtx[i]->cluster->nvtx;
0a2c81 306             } else {
842378 307 //                fprintf(fh,"-1 ");
SP 308                 int_vector[offset]=-1;
0a2c81 309             }
SP 310             }
311         //polymeres
312         if(poly){
313             poly_idx=vlist->n;
314             for(i=0;i<vesicle->poly_list->n;i++){
842378 315                 for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
SP 316                     //fprintf(fh,"-1 ");
317                     int_vector[offset]=-1;
0a2c81 318                 }
SP 319             }
320         }
321         //filaments
322         if(fil){
323             poly_idx=vlist->n+monono*polyno;
324             for(i=0;i<vesicle->filament_list->n;i++){
842378 325                 for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
SP 326     //                    fprintf(fh,"-1 ");
327                     int_vector[offset]=-1;
0a2c81 328                 }
SP 329             }
330         }
842378 331         char *printout=ts_compress((char *)int_vector,offset*sizeof(long));
SP 332         fprintf(fh,"%s",printout);
333         free(printout);
0a2c81 334
SP 335         fprintf(fh,"</DataArray>\n");
336
337
338     }
339
45c708 340     //here comes additional data as needed. Currently only spontaneous curvature
842378 341     fprintf(fh,"<DataArray type=\"Float64\" Name=\"spontaneous_curvature\" format=\"binary\">");
SP 342     offset=0;
343     for(i=0;i<vlist->n;i++, offset++){
344 //        fprintf(fh,"%.17e ",vtx[i]->c);
345         float_vector[offset]=vtx[i]->c;
45c708 346     }
SP 347         //polymeres
348         if(poly){
349             poly_idx=vlist->n;
350             for(i=0;i<vesicle->poly_list->n;i++){
842378 351                 for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
SP 352                     //fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->c);
353                     float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->c;
45c708 354                 }
SP 355             }
356         }
357         //filaments
358         if(fil){
359             poly_idx=vlist->n+monono*polyno;
360             for(i=0;i<vesicle->filament_list->n;i++){
842378 361                 for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
45c708 362         //    fprintf(stderr,"was here\n");
842378 363                     //fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->c);
SP 364                     float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->c;
45c708 365                 }
SP 366             }
367         }
842378 368         printout=ts_compress((char *)float_vector,offset*sizeof(double));
SP 369         fprintf(fh,"%s",printout);
370         free(printout);
371
372
45c708 373     fprintf(fh,"</DataArray>\n");
SP 374
5c64e2 375     //here comes additional data. Energy!
842378 376     fprintf(fh,"<DataArray type=\"Float64\" Name=\"bending_energy\" format=\"binary\">");
SP 377     offset=0;
378     for(i=0;i<vlist->n;i++,offset++){
379 //        fprintf(fh,"%.17e ",vtx[i]->energy*vtx[i]->xk);
380         float_vector[offset]=vtx[i]->energy*vtx[i]->xk;
5c64e2 381     }
SP 382         //polymeres
383         if(poly){
384             poly_idx=vlist->n;
385             for(i=0;i<vesicle->poly_list->n;i++){
842378 386                 for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
SP 387                     //fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k);
388                     float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k;
5c64e2 389                 }
SP 390             }
391         }
392         //filaments
393         if(fil){
394             poly_idx=vlist->n+monono*polyno;
395             for(i=0;i<vesicle->filament_list->n;i++){
842378 396                 for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
5c64e2 397         //    fprintf(stderr,"was here\n");
842378 398 //                    fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k);
SP 399                     float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k;
5c64e2 400                 }
SP 401             }
402         }
842378 403
SP 404     printout=ts_compress((char *)float_vector,offset*sizeof(double));
405         fprintf(fh,"%s",printout);
406         free(printout);
407
5c64e2 408     fprintf(fh,"</DataArray>\n");
SP 409
45c708 410     
13d445 411     fprintf(fh,"</PointData>\n<CellData>\n");
SP 412
7ec6fb 413     if(vesicle->tape->stretchswitch==1){
842378 414         fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"binary\">");
SP 415         offset=0;
416         for(i=0;i<blist->n;i++,offset++){
417             //fprintf(fh, "0.0 ");
418             float_cell_vector[offset]=0.0;
cf7aba 419         }
842378 420         for(i=0;i<monono*polyno+filno*(fonono-1);i++,offset++){
SP 421             float_cell_vector[offset]=0.0;
422             //fprintf(fh,"0.0 ");
7ec6fb 423         }
842378 424         for(i=0;i<vesicle->tlist->n;i++,offset++){
SP 425             float_cell_vector[offset]=vesicle->tlist->tria[i]->energy;
426             //fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
13d445 427         }
842378 428
SP 429         printout=ts_compress((char *)float_cell_vector,offset*sizeof(double));
430         fprintf(fh,"%s",printout);
431         free(printout);
432
7ec6fb 433         fprintf(fh,"</DataArray>\n");
SP 434     }
45c708 435
13d445 436
SP 437
842378 438     fprintf(fh,"</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"binary\">\n");
SP 439     offset=0;
440     for(i=0;i<vlist->n;i++,offset+=3){
441         //fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
442         coordinate_vector[offset]=vtx[i]->x;
443         coordinate_vector[offset+1]=vtx[i]->y;
444         coordinate_vector[offset+2]=vtx[i]->z;
40aa5b 445     }
SP 446     //polymeres
447     if(poly){
448         for(i=0;i<vesicle->poly_list->n;i++){
842378 449             for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,offset+=3){
SP 450                 //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 );
451         coordinate_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->x;
452         coordinate_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[j]->y;
453         coordinate_vector[offset+2]=vesicle->poly_list->poly[i]->vlist->vtx[j]->z;
58230a 454             }
M 455         }
456     }
457     //filaments
458     if(fil){
842378 459         for(i=0;i<vesicle->filament_list->n;i++,offset+=3){
58230a 460             for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
842378 461     //        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 );
SP 462         coordinate_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->x;
463         coordinate_vector[offset+1]=vesicle->filament_list->poly[i]->vlist->vtx[j]->y;
464         coordinate_vector[offset+2]=vesicle->filament_list->poly[i]->vlist->vtx[j]->z;
40aa5b 465             }
SP 466         }
d7639a 467     }
842378 468     printout=ts_compress((char *)coordinate_vector,offset*sizeof(double));
SP 469         fprintf(fh,"%s",printout);
470         free(printout);
d7639a 471
842378 472
SP 473
474     fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"binary\">");
475     free(int_vector);
476     int_vector=(long *)malloc(((blist->n+monono*polyno+filno*(fonono-1))*2+vesicle->tlist->n*3)*sizeof(long));
477     offset=0;
478     for(i=0;i<blist->n;i++,offset+=2){
479 //            fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
480             int_vector[offset]=blist->bond[i]->vtx1->idx;
481             int_vector[offset+1]=blist->bond[i]->vtx2->idx;
d7639a 482     }
40aa5b 483     //polymeres
SP 484     if(poly){
3c1ac1 485         poly_idx=vlist->n;
40aa5b 486         for(i=0;i<vesicle->poly_list->n;i++){
842378 487             for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++,offset+=2){
SP 488                 //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);
489                 int_vector[offset]= vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono;
490                 int_vector[offset+1]=vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono;
491                 
40aa5b 492             }
SP 493     //grafted bonds
842378 494         //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);
SP 495         int_vector[offset]=vesicle->poly_list->poly[i]->grafted_vtx->idx;
496         int_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[0]->idx+vlist->n+i*monono;
497         offset+=2;
40aa5b 498         }
SP 499
500     }
501     
58230a 502     //filaments
M 503     if(fil){
504         poly_idx=vlist->n+monono*polyno;
505         for(i=0;i<vesicle->filament_list->n;i++){
842378 506             for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++,offset+=2){
SP 507 //                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);
508                 int_vector[offset]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono;
509                 int_vector[offset+1]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono;
58230a 510 //        fprintf(stderr,"was here\n");
M 511             
512             }
513         }
514
515     }
842378 516     for(i=0;i<vesicle->tlist->n;i++,offset+=3){
SP 517 //        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);
518         int_vector[offset]=vesicle->tlist->tria[i]->vertex[0]->idx;
519         int_vector[offset+1]=vesicle->tlist->tria[i]->vertex[1]->idx;
520         int_vector[offset+2]=vesicle->tlist->tria[i]->vertex[2]->idx;
521         
8db203 522     }
842378 523         printout=ts_compress((char *)int_vector,offset*sizeof(long));
SP 524         fprintf(fh,"%s",printout);
525         free(printout);
526
527
528     fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"binary\">");
529     offset=0;
530     for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2,offset++){
531 //    fprintf(fh,"%u ",i);
532     int_vector[offset]=i;
d7639a 533     }
842378 534     for(j=i+1;j<i+3*(vesicle->tlist->n);j+=3,offset++){ //let's continue counting from where we left of
SP 535 //        fprintf(fh,"%u ", j);
536     int_vector[offset]=j;
8db203 537     }
842378 538 //    fprintf(fh,"\n");
SP 539     printout=ts_compress((char *)int_vector,offset*sizeof(long));
540     fprintf(fh,"%s",printout);
541     free(printout);
542     fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"binary\">\n");
543     free(int_vector);
544     free(float_vector);
545     free(float_cell_vector);
546     free(coordinate_vector);
547     unsigned char *char_vector=(unsigned char *)malloc((blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n)*sizeof(char));
548     offset=0;
549      for (i=0;i<blist->n+monono*polyno+(fonono-1)*filno;i++,offset++){
550         //fprintf(fh,"3 ");
551     char_vector[offset]=3;
d7639a 552     }
842378 553     for(i=0;i<vesicle->tlist->n;i++,offset++){
SP 554     //    fprintf(fh,"5 ");
555         char_vector[offset]=5;
8db203 556     }
842378 557     printout=ts_compress((char *)char_vector,offset*sizeof(char));
SP 558     fprintf(fh,"%s",printout);
559     free(printout);
560     free(char_vector);
561
d7639a 562     fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
SP 563     fclose(fh);
564     return TS_SUCCESS;
565
566 }
567
144784 568 ts_bool write_pov_file(ts_vesicle *vesicle, char *filename){
SP 569     FILE *fh;
570     ts_uint i;
571     
572     fh=fopen(filename, "w");
573     if(fh==NULL){
574         err("Cannot open file %s for writing");
575         return TS_FAIL;
576     }
577
578     for(i=0;i<vesicle->tlist->n;i++){
579     
580     fprintf(fh,"\ttriangle {");
581     fprintf(fh,"\t<%e,%e,%e> <%e,%e,%e> <%e,%e,%e> }\n", 
582     vesicle->tlist->tria[i]->vertex[0]->x,
583     vesicle->tlist->tria[i]->vertex[0]->y,
584     vesicle->tlist->tria[i]->vertex[0]->z,
585
586     vesicle->tlist->tria[i]->vertex[1]->x,
587     vesicle->tlist->tria[i]->vertex[1]->y,
588     vesicle->tlist->tria[i]->vertex[1]->z,
589
590     vesicle->tlist->tria[i]->vertex[2]->x,
591     vesicle->tlist->tria[i]->vertex[2]->y,
592     vesicle->tlist->tria[i]->vertex[2]->z
593     );
594     }
595         
596     fclose(fh);
597     return TS_SUCCESS;
598 }
599
600
1ab449 601 ts_tape *parsetape(char *filename){
698ae1 602     FILE *fd = fopen (filename, "r");
SP 603     long length;
604     size_t size;
605     fseek (fd, 0, SEEK_END);
606       length = ftell (fd);
607     fseek (fd, 0, SEEK_SET);
608     size=fread (tapetxt, 1, length, fd);
609     fclose(fd);
610     if(size);
611     ts_tape *tape=parsetapebuffer(tapetxt);
612     return tape;
613 }
614
615 ts_tape *parsetapebuffer(char *buffer){
1ab449 616     ts_tape *tape=(ts_tape *)calloc(1,sizeof(ts_tape));
SP 617     tape->multiprocessing=calloc(255,sizeof(char));
40aa5b 618
d7639a 619     cfg_opt_t opts[] = {
1ab449 620         CFG_SIMPLE_INT("nshell", &tape->nshell),
SP 621         CFG_SIMPLE_INT("npoly", &tape->npoly),
622         CFG_SIMPLE_INT("nmono", &tape->nmono),
58230a 623     CFG_SIMPLE_INT("nfil",&tape->nfil),
M 624     CFG_SIMPLE_INT("nfono",&tape->nfono),
e98482 625     CFG_SIMPLE_INT("internal_poly",&tape->internal_poly),
58230a 626     CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
37791b 627     CFG_SIMPLE_FLOAT("R_nucleusX",&tape->R_nucleusX),
SP 628     CFG_SIMPLE_FLOAT("R_nucleusY",&tape->R_nucleusY),
629     CFG_SIMPLE_FLOAT("R_nucleusZ",&tape->R_nucleusZ),
58230a 630     CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
ea1cce 631     CFG_SIMPLE_FLOAT("dmin_interspecies", &tape->dmin_interspecies),
1ab449 632         CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
SP 633     CFG_SIMPLE_INT("pswitch",&tape->pswitch),
9166cb 634     CFG_SIMPLE_INT("constvolswitch",&tape->constvolswitch),
c0ae90 635     CFG_SIMPLE_INT("constareaswitch",&tape->constareaswitch),
43c042 636     CFG_SIMPLE_FLOAT("constvolprecision",&tape->constvolprecision),
2ae815 637     CFG_SIMPLE_INT("stretchswitch",&tape->stretchswitch),
SP 638     CFG_SIMPLE_FLOAT("xkA0",&tape->xkA0),    
1ab449 639     CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
SP 640     CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
b30f45 641     CFG_SIMPLE_FLOAT("xi",&tape->xi),
1ab449 642         CFG_SIMPLE_FLOAT("stepsize",&tape->stepsize),
SP 643         CFG_SIMPLE_INT("nxmax", &tape->ncxmax),
644         CFG_SIMPLE_INT("nymax", &tape->ncymax),
645         CFG_SIMPLE_INT("nzmax", &tape->nczmax),
646         CFG_SIMPLE_INT("iterations",&tape->iterations),
647     CFG_SIMPLE_INT("mcsweeps",&tape->mcsweeps),
648     CFG_SIMPLE_INT("inititer", &tape->inititer),
819a09 649                 CFG_SIMPLE_BOOL("quiet",(cfg_bool_t *)&tape->quiet),
S 650         CFG_SIMPLE_STR("multiprocessing",&tape->multiprocessing),
1ab449 651         CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
SP 652         CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
653         CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
dc77e8 654     CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
e5858f 655     CFG_SIMPLE_INT("number_of_vertices_with_c0", &tape->number_of_vertices_with_c0),
SP 656     CFG_SIMPLE_FLOAT("c0",&tape->c0),
657     CFG_SIMPLE_FLOAT("w",&tape->w),
250de4 658     CFG_SIMPLE_FLOAT("F",&tape->F),
6c274d 659 /* Variables related to plane confinement */
514ebc 660     CFG_INT("plane_confinement_switch", 0, CFGF_NONE),
SP 661     CFG_FLOAT("plane_d", 15, CFGF_NONE),
662     CFG_FLOAT("plane_F", 1000, CFGF_NONE),
6c274d 663 /* Variables related to stretching */
36bc6d 664 //    CFG_FLOAT("stretchswitch", 0, CFGF_NONE),
SP 665 //    CFG_FLOAT("xkA0",0,CFGF_NONE),
d7639a 666         CFG_END()
SP 667     };
668     cfg_t *cfg;    
669     ts_uint retval;
670     cfg = cfg_init(opts, 0);
698ae1 671     retval=cfg_parse_buf(cfg, buffer);
514ebc 672     tape->plane_confinement_switch=cfg_getint(cfg,"plane_confinement_switch");
88bdd7 673     tape->plane_d=cfg_getfloat(cfg,"plane_d");
SP 674     tape->plane_F=cfg_getfloat(cfg,"plane_F");
514ebc 675
d7639a 676     if(retval==CFG_FILE_ERROR){
a6b1b5 677     fatal("No tape file.",100);
d7639a 678     }
SP 679     else if(retval==CFG_PARSE_ERROR){
680     fatal("Invalid tape!",100);
681     }
a2db52 682
7b0c07 683     /* here we override all values read from tape with values from commandline*/
SP 684     getcmdline_tape(cfg,command_line_args.tape_opts);
d7639a 685     cfg_free(cfg);
40aa5b 686
SP 687
1ab449 688     /* global variables are set automatically */
SP 689     quiet=tape->quiet;
690     return tape;
691 }
d7639a 692
1ab449 693 ts_bool tape_free(ts_tape *tape){
SP 694     free(tape->multiprocessing);
695     free(tape);
696     return TS_SUCCESS;
d7639a 697 }
f74313 698
SP 699
7b0c07 700
SP 701 ts_bool getcmdline_tape(cfg_t *cfg, char *opts){
702
703     char *commands, *backup, *saveptr, *saveopptr, *command, *operator[2];
07e3de 704     operator[0]=0;
SP 705     operator[1]=0;
7b0c07 706     ts_uint i,j;
SP 707     commands=(char *)malloc(10000*sizeof(char));
708     backup=commands; //since the pointer to commands will be lost, we acquire a pointer that will serve as backup.
709     strcpy(commands,opts);
710     for(i=0; ;i++, commands=NULL){
711         //breaks comma separated list of commands into specific commands.
712         command=strtok_r(commands,",",&saveptr);    
713         if(command==NULL) break;
714 //        fprintf(stdout,"Command %d: %s\n",i,command);    
715         //extracts name of command and value of command into operator[2] array.
716         for(j=0; j<2;j++,command=NULL){
717             operator[j]=strtok_r(command,"=",&saveopptr);
718             if(operator[j]==NULL) break;
719 //            fprintf(stdout," ---> Operator %d: %s\n",j,operator[j]);        
720         }
721         //1. check: must have 2 operators.
722         if(j!=2) fatal("Error. Command line tape options are not formatted properly",1);
723
724     //    cfg_setstr(cfg,operator[0],operator[1]);
725         cmdline_to_tape(cfg,operator[0],operator[1]);
726         //2. check: must be named properly.
727         //3. check: must be of right format (integer, double, string, ...)
728         
729     }
730     free(backup);
731     return TS_SUCCESS;
732 }
733
734
735 ts_bool cmdline_to_tape(cfg_t *cfg, char *key, char *val){
736
737     cfg_opt_t *cfg_opt=cfg_getopt(cfg,key);
738     if(cfg_opt==NULL) fatal("Commandline tape option not recognised",1); //return TS_FAIL; 
739     switch (cfg_opt->type){
740         case CFGT_INT:
741             cfg_setint(cfg,key,atol(val));
742             break;
743         case CFGT_FLOAT:
744             cfg_setfloat(cfg,key,atof(val));
745             break;
746 /*        case CFGT_BOOL:
747             cfg_setbool(cfg,operator[0],operator[1]);
748             break; */
749         case CFGT_STR:
750             cfg_setstr(cfg,key,val);
751             break;
752         default:
753             break;
754
755     }
f74313 756     return TS_SUCCESS;
SP 757 }