Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2016-05-24 cc95bd1ca5e69ce13766ca0acea11560531ed620
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 /** DUMP STATE TO DISK DRIVE **/
SP 21
1ab449 22 ts_bool dump_state(ts_vesicle *vesicle, ts_uint iteration){
e9eab4 23
SP 24     /* save current state with wrong pointers. Will fix that later */
25     ts_uint i,j,k;
267db5 26     FILE *fh=fopen(command_line_args.dump_fullfilename,"wb");
e9eab4 27
SP 28     /* dump vesicle */
d0a1d3 29     fwrite(vesicle, sizeof(ts_vesicle)-sizeof(ts_double),1,fh);
e9eab4 30     /* dump vertex list */
SP 31     fwrite(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
32     /* dump bond list */
33     fwrite(vesicle->blist, sizeof(ts_bond_list),1,fh);
34     /* dump triangle list */
35     fwrite(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
36     /* dump cell list */
37     fwrite(vesicle->clist, sizeof(ts_cell_list),1,fh);
38     /* dump poly list */
39     fwrite(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
bcf455 40     /* dump filament list */
M 41     fwrite(vesicle->filament_list, sizeof(ts_poly_list),1,fh);
e9eab4 42     /* level 1 complete */
SP 43
44     /*dump vertices*/
45     for(i=0;i<vesicle->vlist->n;i++){
46         fwrite(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
47         /* dump pointer offsets for:
48                     neigh
49                     bond
50                     tria
51                     cell is ignored
52         */
53         for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
3c1ac1 54             fwrite(&vesicle->vlist->vtx[i]->neigh[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 55         }
SP 56         for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
3c1ac1 57             fwrite(&vesicle->vlist->vtx[i]->bond[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 58         }
SP 59         for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
3c1ac1 60             fwrite(&vesicle->vlist->vtx[i]->tristar[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 61         }
SP 62     }
63
64     /*dump bonds*/
65     for(i=0;i<vesicle->blist->n;i++){
66         fwrite(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
67         /* dump pointer offsets for vtx1 and vtx2 */
3c1ac1 68         //off=(ts_ulong)(vesicle->blist->bond[i]->vtx1-vesicle->vlist->vtx[0]);
SP 69         fwrite(&vesicle->blist->bond[i]->vtx1->idx,sizeof(ts_uint),1,fh); 
70         //off=(ts_ulong)(vesicle->blist->bond[i]->vtx2-vesicle->vlist->vtx[0]);
71         fwrite(&vesicle->blist->bond[i]->vtx2->idx,sizeof(ts_uint),1,fh); 
e9eab4 72     }
SP 73
74     /*dump triangles*/
75     for(i=0;i<vesicle->tlist->n;i++){
76         fwrite(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
77         /* dump pointer offsets for vertex */
3c1ac1 78         fwrite(&vesicle->tlist->tria[i]->vertex[0]->idx,sizeof(ts_uint),1,fh); 
SP 79         fwrite(&vesicle->tlist->tria[i]->vertex[1]->idx,sizeof(ts_uint),1,fh); 
80         fwrite(&vesicle->tlist->tria[i]->vertex[2]->idx,sizeof(ts_uint),1,fh); 
e9eab4 81         /* dump pointer offsets for neigh */
SP 82         for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
3c1ac1 83             fwrite(&vesicle->tlist->tria[i]->neigh[j]->idx,sizeof(ts_uint),1,fh); 
e9eab4 84         }
SP 85     }
86
87
88     /*dump polymeres */
89     for(i=0;i<vesicle->poly_list->n;i++){
90         fwrite(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
3c1ac1 91         fwrite(vesicle->poly_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
SP 92         fwrite(vesicle->poly_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
e9eab4 93     } 
SP 94      
95     /* dump poly vertex(monomer) list*/
96     for(i=0;i<vesicle->poly_list->n;i++){
97         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
98             fwrite(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
99             /* dump offset for neigh and bond */
100             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
3c1ac1 101                // off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]-vesicle->poly_list->poly[i]->vlist->vtx[0]);
SP 102                 fwrite(&vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]->idx,sizeof(ts_uint),1,fh); 
e9eab4 103             }
SP 104             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
3c1ac1 105                 //off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]-vesicle->poly_list->poly[i]->blist->bond[0]);
SP 106                 fwrite(&vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]->idx,sizeof(ts_uint),1,fh); 
e9eab4 107             }
SP 108         }
3c1ac1 109     // grafted vtx on vesicle data dump
SP 110         fwrite(&vesicle->poly_list->poly[i]->grafted_vtx->idx, sizeof(ts_uint),1,fh);
e9eab4 111     }
SP 112     /* dump poly bonds between monomers list*/
113     for(i=0;i<vesicle->poly_list->n;i++){
114         for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
115             fwrite(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
116             /* dump vtx1 and vtx2 offsets */
3c1ac1 117             //off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx1-vesicle->poly_list->poly[i]->vlist->vtx[0]);
SP 118             fwrite(&vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx,sizeof(ts_uint),1,fh); 
119 //            off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]);
120             fwrite(&vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx,sizeof(ts_uint),1,fh); 
e9eab4 121         }
SP 122     }
123
bcf455 124
M 125   /*dump filamentes grandes svinjas */
126     for(i=0;i<vesicle->filament_list->n;i++){
127         fwrite(vesicle->filament_list->poly[i],sizeof(ts_poly),1,fh);
128         fwrite(vesicle->filament_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
129         fwrite(vesicle->filament_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
130     } 
131      
132     /* dump filamentes vertex(monomer) list*/
133     for(i=0;i<vesicle->filament_list->n;i++){
134         for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
135             fwrite(vesicle->filament_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
136             /* dump offset for neigh and bond */
137             for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
138                // off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]-vesicle->poly_list->poly[i]->vlist->vtx[0]);
139                 fwrite(&vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh[k]->idx,sizeof(ts_uint),1,fh); 
140             }
141             for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no;k++){
142                 //off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]-vesicle->poly_list->poly[i]->blist->bond[0]);
143                 fwrite(&vesicle->filament_list->poly[i]->vlist->vtx[j]->bond[k]->idx,sizeof(ts_uint),1,fh); 
144             }
145         }
146     }
147     /* dump poly bonds between monomers list*/
148     for(i=0;i<vesicle->filament_list->n;i++){
149         for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
150             fwrite(vesicle->filament_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
151             /* dump vtx1 and vtx2 offsets */
152             //off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx1-vesicle->poly_list->poly[i]->vlist->vtx[0]);
153             fwrite(&vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx,sizeof(ts_uint),1,fh); 
154 //            off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]);
155             fwrite(&vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx,sizeof(ts_uint),1,fh); 
156         }
157     }
158
159
160
e9eab4 161 /* pointer offsets for fixing the restored pointers */
SP 162 /* need pointers for 
163     vlist->vtx
164     blist->bond
165     tlist->tria
166     clist->cell
167     poly_list->poly
168     and for each poly:
169         poly_list->poly->vtx
170         poly_list->poly->bond
171 */
172
bd2210 173 //    fwrite(vesicle->clist, sizeof(ts_cell_list),1,  fh);
9166cb 174 /* write tape information on vesicle */
bd2210 175 //    fwrite(vesicle->tape,sizeof(ts_tape),1,fh);
1ab449 176     fwrite(&iteration, sizeof(ts_uint),1,fh);
e9eab4 177     fclose(fh);
SP 178     return TS_SUCCESS;
179 }
180
181
182 /** RESTORE DUMP FROM DISK **/
1ab449 183 ts_vesicle *restore_state(ts_uint *iteration){
e9eab4 184     ts_uint i,j,k;
267db5 185     FILE *fh=fopen(command_line_args.dump_fullfilename,"rb");
3f5c83 186
SP 187     struct stat sb;
267db5 188     if (stat(command_line_args.dump_fullfilename, &sb) == -1) {
3f5c83 189         //dump file does not exist.
SP 190         return NULL;
191     }
192
193     //check if it is regular file
194     if((sb.st_mode & S_IFMT) != S_IFREG) {
195         //dump file is not a regular file.
196         ts_fprintf(stderr,"Dump file is not a regular file!\n");
197         return NULL;
198     }
199
e9eab4 200     ts_uint retval;
3c1ac1 201     ts_uint idx;
e9eab4 202
SP 203 /* we restore all the data from the dump */
204     /* restore vesicle */
205     ts_vesicle *vesicle=(ts_vesicle *)calloc(1,sizeof(ts_vesicle));
d0a1d3 206     retval=fread(vesicle, sizeof(ts_vesicle)-sizeof(ts_double),1,fh);
86f5e7 207 //    fprintf(stderr,"was here! %e\n",vesicle->dmax);
SP 208
e9eab4 209     /* restore vertex list */
SP 210     vesicle->vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
211     retval=fread(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
212     /* restore bond list */
213     vesicle->blist=(ts_bond_list *)malloc(sizeof(ts_bond_list));
214     retval=fread(vesicle->blist, sizeof(ts_bond_list),1,fh);
215     /* restore triangle list */
216     vesicle->tlist=(ts_triangle_list *)malloc(sizeof(ts_triangle_list));
217     retval=fread(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
218     /* restore cell list */
219     vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
220     retval=fread(vesicle->clist, sizeof(ts_cell_list),1,fh);
221     /* restore poly list */
222     vesicle->poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
223     retval=fread(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
bcf455 224     /* restore filament list */
M 225     vesicle->filament_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
226     retval=fread(vesicle->filament_list, sizeof(ts_poly_list),1,fh);
e9eab4 227     /* level 1 complete */
SP 228
3c1ac1 229 /* prerequisity. Bonds must be malloced before vertexes are recreated */
SP 230   vesicle->blist->bond=(ts_bond **)calloc(vesicle->blist->n,sizeof(ts_bond *));
231     for(i=0;i<vesicle->blist->n;i++){
232         vesicle->blist->bond[i]=(ts_bond *)malloc(sizeof(ts_bond));
233     }
234 /* prerequisity. Triangles must be malloced before vertexes are recreated */
235   vesicle->tlist->tria=(ts_triangle **)calloc(vesicle->tlist->n,sizeof(ts_triangle *));
236     for(i=0;i<vesicle->tlist->n;i++){
237         vesicle->tlist->tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle));
238 }
a37ba2 239 /* prerequisity. Vertices must be malloced before vertexes are recreated */
e9eab4 240     vesicle->vlist->vtx=(ts_vertex **)calloc(vesicle->vlist->n,sizeof(ts_vertex *));
a37ba2 241  for(i=0;i<vesicle->vlist->n;i++){
e9eab4 242         vesicle->vlist->vtx[i]=(ts_vertex *)malloc(sizeof(ts_vertex));
a37ba2 243  }
SP 244   /*restore vertices*/
245     for(i=0;i<vesicle->vlist->n;i++){
e9eab4 246         retval=fread(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
SP 247         /*restore neigh, bond, tristar. Ignoring cell */
248         vesicle->vlist->vtx[i]->neigh=(ts_vertex **)calloc(vesicle->vlist->vtx[i]->neigh_no, sizeof(ts_vertex *));
249         for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
3c1ac1 250             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 251             vesicle->vlist->vtx[i]->neigh[j]=vesicle->vlist->vtx[idx];
e9eab4 252         }
SP 253         vesicle->vlist->vtx[i]->bond=(ts_bond **)calloc(vesicle->vlist->vtx[i]->bond_no, sizeof(ts_bond *));
254         for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
3c1ac1 255             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 256 /* 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 */
257             vesicle->vlist->vtx[i]->bond[j]=vesicle->blist->bond[idx];    
e9eab4 258         }
SP 259
260         vesicle->vlist->vtx[i]->tristar=(ts_triangle **)calloc(vesicle->vlist->vtx[i]->tristar_no, sizeof(ts_triangle *));
261         for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
3c1ac1 262             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 263 /* same comment as above */
264             vesicle->vlist->vtx[i]->tristar[j]=vesicle->tlist->tria[idx];
e9eab4 265         }
SP 266
267     }
268
269     /*restore bonds*/
3c1ac1 270    // vesicle->blist->bond=(ts_bond **)calloc(vesicle->blist->n,sizeof(ts_bond *)); // done before.
e9eab4 271     for(i=0;i<vesicle->blist->n;i++){
3c1ac1 272      //   vesicle->blist->bond[i]=(ts_bond *)malloc(sizeof(ts_bond)); //done before.
e9eab4 273         retval=fread(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
SP 274         /* restore vtx1 and vtx2 */
3c1ac1 275         retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 276         vesicle->blist->bond[i]->vtx1=vesicle->vlist->vtx[idx];
277         retval=fread(&idx,sizeof(ts_uint),1,fh);
278         vesicle->blist->bond[i]->vtx2=vesicle->vlist->vtx[idx];
e9eab4 279     }
SP 280
281     /*restore triangles*/
3c1ac1 282 //    vesicle->tlist->tria=(ts_triangle **)calloc(vesicle->tlist->n,sizeof(ts_triangle *)); // done before
e9eab4 283     for(i=0;i<vesicle->tlist->n;i++){
3c1ac1 284  //       vesicle->tlist->tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle)); // done before
e9eab4 285         retval=fread(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
SP 286         /* restore pointers for vertices */
3c1ac1 287         retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 288         vesicle->tlist->tria[i]->vertex[0]=vesicle->vlist->vtx[idx];
289         retval=fread(&idx,sizeof(ts_uint),1,fh);
290         vesicle->tlist->tria[i]->vertex[1]=vesicle->vlist->vtx[idx];
291         retval=fread(&idx,sizeof(ts_uint),1,fh);
292         vesicle->tlist->tria[i]->vertex[2]=vesicle->vlist->vtx[idx];
e9eab4 293         /* restore pointers for neigh */
3c1ac1 294      vesicle->tlist->tria[i]->neigh=(ts_triangle **)malloc(vesicle->tlist->tria[i]->neigh_no*sizeof(ts_triangle *));
e9eab4 295         for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
3c1ac1 296             retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 297             vesicle->tlist->tria[i]->neigh[j]=vesicle->tlist->tria[idx];
e9eab4 298         }
SP 299
300     }
301    
302     /*restore cells */
303 /*TODO: do we need to recalculate cells here? */
304 /*    vesicle->clist->cell=(ts_cell **)malloc(vesicle->clist->cellno*sizeof(ts_cell *));
305     for(i=0;i<vesicle->clist->cellno;i++){
306         vesicle->clist->cell[i]=(ts_cell *)malloc(sizeof(ts_cell));
307         retval=fread(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
308     }
309 */
310     /*restore polymeres */
311     vesicle->poly_list->poly = (ts_poly **)calloc(vesicle->poly_list->n,sizeof(ts_poly *));
312     for(i=0;i<vesicle->poly_list->n;i++){
313         vesicle->poly_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly));
314         retval=fread(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
3c1ac1 315         vesicle->poly_list->poly[i]->vlist=(ts_vertex_list *)calloc(1,sizeof(ts_vertex_list));
SP 316         retval=fread(vesicle->poly_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
317         vesicle->poly_list->poly[i]->blist=(ts_bond_list *)calloc(1,sizeof(ts_bond_list));
318         retval=fread(vesicle->poly_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
319     /* initialize adress space for pointers that will hold specific vertices (monomers) and bonds */
320         vesicle->poly_list->poly[i]->vlist->vtx=(ts_vertex **)calloc(vesicle->poly_list->poly[i]->vlist->n,sizeof(ts_vertex *));
321         vesicle->poly_list->poly[i]->blist->bond=(ts_bond **)calloc(vesicle->poly_list->poly[i]->blist->n,sizeof(ts_bond *));
322      for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
323             vesicle->poly_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex));
324     }
325     for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
326             vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
327     }
328
e9eab4 329     } 
3c1ac1 330
e9eab4 331      
SP 332     /* restore poly vertex(monomer) list*/
333     for(i=0;i<vesicle->poly_list->n;i++){
334         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
335             retval=fread(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
3c1ac1 336                 
e9eab4 337             /* restore neigh and bonds */
SP 338             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 *));
339             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
3c1ac1 340                 retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 341                 vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->poly_list->poly[i]->vlist->vtx[idx];
e9eab4 342             }
SP 343             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 *));
344             for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
3c1ac1 345                 retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 346                 vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->poly_list->poly[i]->blist->bond[idx];
e9eab4 347             }
SP 348
349         }
3c1ac1 350     /* restore grafted vtx on vesicle and grafted_poly */
SP 351                 retval=fread(&idx,sizeof(ts_uint),1,fh);
352         vesicle->vlist->vtx[idx]->grafted_poly=vesicle->poly_list->poly[i];
353         vesicle->poly_list->poly[i]->grafted_vtx=vesicle->vlist->vtx[idx];    
e9eab4 354     }
3c1ac1 355
e9eab4 356     /* restore poly bonds between monomers list*/
SP 357     for(i=0;i<vesicle->poly_list->n;i++){
358         for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
3c1ac1 359        //     vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
e9eab4 360             retval=fread(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
SP 361             /* restore vtx1 and vtx2 */
3c1ac1 362                 retval=fread(&idx,sizeof(ts_uint),1,fh);
SP 363                 vesicle->poly_list->poly[i]->blist->bond[j]->vtx1=vesicle->poly_list->poly[i]->vlist->vtx[idx];
364                 retval=fread(&idx,sizeof(ts_uint),1,fh);
365                 vesicle->poly_list->poly[i]->blist->bond[j]->vtx2=vesicle->poly_list->poly[i]->vlist->vtx[idx];
e9eab4 366         }
SP 367     }
368
bcf455 369     /*restore filaments */
M 370     vesicle->filament_list->poly = (ts_poly **)calloc(vesicle->filament_list->n,sizeof(ts_poly *));
371     for(i=0;i<vesicle->filament_list->n;i++){
372         vesicle->filament_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly));
373         retval=fread(vesicle->filament_list->poly[i],sizeof(ts_poly),1,fh);
374         vesicle->filament_list->poly[i]->vlist=(ts_vertex_list *)calloc(1,sizeof(ts_vertex_list));
375         retval=fread(vesicle->filament_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
376         vesicle->filament_list->poly[i]->blist=(ts_bond_list *)calloc(1,sizeof(ts_bond_list));
377         retval=fread(vesicle->filament_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
378     /* initialize adress space for pointers that will hold specific vertices (monomers) and bonds */
379         vesicle->filament_list->poly[i]->vlist->vtx=(ts_vertex **)calloc(vesicle->filament_list->poly[i]->vlist->n,sizeof(ts_vertex *));
380         vesicle->filament_list->poly[i]->blist->bond=(ts_bond **)calloc(vesicle->filament_list->poly[i]->blist->n,sizeof(ts_bond *));
381      for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
382             vesicle->filament_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex));
383     }
384     for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
385             vesicle->filament_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
386     }
387
388     } 
389
390      
391     /* restore poly vertex(monomer) list*/
392     for(i=0;i<vesicle->filament_list->n;i++){
393         for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
394             retval=fread(vesicle->filament_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
395                 
396             /* restore neigh and bonds */
397             vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh=(ts_vertex **)calloc(vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no, sizeof(ts_vertex *));
398             for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
399                 retval=fread(&idx,sizeof(ts_uint),1,fh);
400                 vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->filament_list->poly[i]->vlist->vtx[idx];
401             }
402             vesicle->filament_list->poly[i]->vlist->vtx[j]->bond=(ts_bond **)calloc(vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no, sizeof(ts_bond *));
403             for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no;k++){
404                 retval=fread(&idx,sizeof(ts_uint),1,fh);
405                 vesicle->filament_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->filament_list->poly[i]->blist->bond[idx];
406             }
407
408         }
409     }
410
411     /* restore poly bonds between monomers list*/
412     for(i=0;i<vesicle->filament_list->n;i++){
413         for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
414        //     vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
415             retval=fread(vesicle->filament_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
416             /* restore vtx1 and vtx2 */
417                 retval=fread(&idx,sizeof(ts_uint),1,fh);
418                 vesicle->filament_list->poly[i]->blist->bond[j]->vtx1=vesicle->filament_list->poly[i]->vlist->vtx[idx];
419                 retval=fread(&idx,sizeof(ts_uint),1,fh);
420                 vesicle->filament_list->poly[i]->blist->bond[j]->vtx2=vesicle->filament_list->poly[i]->vlist->vtx[idx];
421         }
422     }
267db5 423     vesicle->tape=parsetape(command_line_args.tape_fullfilename);
063289 424 // recreating space for cells // 
055dd7 425     vesicle->clist=init_cell_list(vesicle->tape->ncxmax, vesicle->tape->ncymax, vesicle->tape->nczmax, vesicle->tape->stepsize);
364c30 426     vesicle->clist->max_occupancy=8;
47a7ac 427 //    vesicle->tape=(ts_tape *)malloc(sizeof(ts_tape));
SP 428 //    retval=fread(vesicle->tape, sizeof(ts_tape),1,fh);
1ab449 429     retval=fread(iteration,sizeof(ts_uint),1,fh);
e9eab4 430     if(retval); 
SP 431     fclose(fh);
432     return vesicle;
433 }
434
435
436
083e03 437 ts_bool parse_args(int argc, char **argv){
3f5c83 438     int c, retval;
SP 439     struct stat sb;
440     sprintf(command_line_args.path, "./"); //clear string;
7b0c07 441     sprintf(command_line_args.output_fullfilename,"output.pvd");
SP 442     sprintf(command_line_args.dump_fullfilename,"dump.bin");
443     sprintf(command_line_args.tape_fullfilename,"tape");
1bf3c3 444     sprintf(command_line_args.tape_templatefull,"./tape");
3f5c83 445             FILE *file;
SP 446     
083e03 447 while (1)
SP 448      {
449        static struct option long_options[] =
450          {
8a6614 451            {"force-from-tape", no_argument,       &(command_line_args.force_from_tape), 1},
SP 452        {"reset-iteration-count", no_argument, &(command_line_args.reset_iteration_count), 1},
083e03 453            {"tape",     no_argument,       0, 't'},
fd8126 454        {"version", no_argument, 0, 'v'},
083e03 455            {"output-file",  required_argument, 0, 'o'},
SP 456            {"directory",  required_argument, 0, 'd'},
3f5c83 457            {"dump-filename", required_argument,0, 'f'},
7b0c07 458            {"tape-options",required_argument,0,'c'},
1bf3c3 459            {"tape-template", required_argument,0,0},
ab798b 460             {"restore-from-vtk",required_argument,0,0},
083e03 461            {0, 0, 0, 0}
SP 462          };
463        /* getopt_long stores the option index here. */
464        int option_index = 0;
465
fd8126 466        c = getopt_long (argc, argv, "d:f:o:t:c:v",
083e03 467                         long_options, &option_index);
SP 468
469        /* Detect the end of the options. */
470        if (c == -1)
471          break;
472
473        switch (c)
474          {
475          case 0:
476            /* If this option set a flag, do nothing else now. */
477            if (long_options[option_index].flag != 0)
478              break;
1bf3c3 479 /*           printf ("option %s", long_options[option_index].name);
083e03 480            if (optarg)
1bf3c3 481              printf (" with arg %s", optarg); 
SP 482            printf ("\n"); */
483             //TODO: find a better way.
484             if(strcmp(long_options[option_index].name,"tape-template")==0){
485                 strcpy(command_line_args.tape_templatefull,optarg);
486             }
ab798b 487             if(strcmp(long_options[option_index].name,"restore-from-vtk")==0){
ee84bd 488                 strcpy(command_line_args.dump_from_vtk,optarg);
SP 489             }
083e03 490            break;
fd8126 491      case 'v':
SP 492         fprintf(stdout,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);
493             fprintf(stdout,"Programming done by: Samo Penic and Miha Fosnaric\n");
494             fprintf(stdout,"Released under terms of GPLv3\n");
495         exit(0);
083e03 496
7b0c07 497          case 'c':
SP 498               strcpy(command_line_args.tape_opts,optarg);
499             break;
500          case 't': //tape
3f5c83 501                 strcpy(command_line_args.tape_fullfilename,optarg);
083e03 502            break;
SP 503
7b0c07 504          case 'o':  //set filename of master pvd output file
3f5c83 505             strcpy(command_line_args.output_fullfilename, optarg);
SP 506             break;
083e03 507
SP 508          case 'd':
509             //check if directory exists. If not create one. If creation is
510             //successful, set directory for output files.
511             //printf ("option -d with value `%s'\n", optarg);
3f5c83 512             if (stat(optarg, &sb) == -1) {
SP 513                 //directory does not exist
083e03 514                 retval=mkdir(optarg, 0700);
SP 515                 if(retval){
3f5c83 516                     fatal("Could not create requested directory. Check if you have permissions",1);
083e03 517                 }
SP 518             }
3f5c83 519             //check if is a proper directory
SP 520             else if((sb.st_mode & S_IFMT) != S_IFDIR) {
521                 //it is not a directory. fatal error.
522                 ts_fprintf(stderr,"%s is not a directory!\n",optarg);
523                 fatal("Cannot continue",1);
083e03 524             }
3f5c83 525             strcpy(command_line_args.path, optarg);
083e03 526            break;
SP 527
528         case 'f':
3f5c83 529             strcpy(command_line_args.dump_fullfilename, optarg);
083e03 530             break;
SP 531
532          case '?':
533            /* getopt_long already printed an error message. */
534
535             ts_fprintf(stderr,"\n\nhere comes the help.\n\n");
536             fatal("Ooops, read help first",1);
537            break;
538
539          default:
540            exit (1);
541          }
542      }
543
7b0c07 544 //Here we set correct values for full filenames!
SP 545     char *buffer=(char *)malloc(10000*sizeof(char));
546     //correct the path and add trailing /
547     if(command_line_args.path[strlen(command_line_args.path)-1]!='/') strcat(command_line_args.path,"/");
548    
549 /* master pvd output file */ 
550     strcpy(buffer,command_line_args.path);
551     strcat(buffer,command_line_args.output_fullfilename);
552     if ((file = fopen(buffer, "w")) == NULL) {
553                 fprintf(stderr,"Could not create output file %s!\n", buffer);
554                 fatal("Please specify correct output file or check permissions of the file",1);
1bf3c3 555                 //there is a tape template. make a copy into desired directory
7b0c07 556                 
SP 557             } else {
558                 fclose(file);
559                 strcpy(command_line_args.output_fullfilename,buffer);
560             }
561
562 /* tape file */
563     strcpy(buffer,command_line_args.path);
564     strcat(buffer,command_line_args.tape_fullfilename);
565     if (stat(buffer, &sb) == -1) {
1bf3c3 566
SP 567                 //tape does not exist. does tape template exist?
568                 if(stat(command_line_args.tape_templatefull, &sb)==-1){ 
569                     ts_fprintf(stderr,"Tape '%s' does not exist and no tape template was specified (or does not exist)!\n",buffer);
570                     fatal("Please select correct tape or check permissions of the file",1);
571                 } else {
572                     //tape template found
573                     fatal("Samo did not program template copy yet",1); 
574                 }
7b0c07 575             } else {
SP 576                 strcpy(command_line_args.tape_fullfilename,buffer);
577             }
578
579
580 /* dump file */
581             strcpy(buffer,command_line_args.path);
582             strcat(buffer,command_line_args.dump_fullfilename);
583             //check if dump file exist first.
584             if (stat(buffer, &sb) == -1) {
585                 //no dump file. check if we can create one.
586                 if ((file = fopen(buffer, "w")) == NULL) {
587                     fprintf(stderr,"Could not create dump file '%s'!\n",buffer);
588                     fatal("Please specify correct dump file or check permissions of the file",1);
589                 } else {
590                     fclose(file);
591                     //good, file is writeable. delete it for now.
592                     remove(buffer);
593                 }
594             }
595             strcpy(command_line_args.dump_fullfilename, buffer);
596
597
598     free(buffer);
083e03 599     return TS_SUCCESS;
SP 600
601 }
602
603
604
d7639a 605 ts_bool print_vertex_list(ts_vertex_list *vlist){
SP 606     ts_uint i;
607     printf("Number of vertices: %u\n",vlist->n);
608     for(i=0;i<vlist->n;i++){
a6b1b5 609         printf("%u: %f %f %f\n",
8f6a69 610 vlist->vtx[i]->idx,vlist->vtx[i]->x, vlist->vtx[i]->y, vlist->vtx[i]->z);
d7639a 611     }
SP 612     return TS_SUCCESS;
613 }
614
615 ts_bool print_vertex_neighbours(ts_vertex_list *vlist){
616     ts_uint i,j;
a6b1b5 617     ts_vertex **vtx=vlist->vtx;
d7639a 618     printf("Vertex id(neigh no): (neighvertex coord) (neighvertex coord) ...\n");
SP 619     for(i=0;i<vlist->n;i++){
8f6a69 620         printf("%u(%u): ",vtx[i]->idx,vtx[i]->neigh_no);
SP 621         for(j=0;j<vtx[i]->neigh_no;j++){
622             printf("(%f,%f,%f)",vtx[i]->neigh[j]->x,
623 vtx[i]->neigh[j]->y,vtx[i]->neigh[j]->z);
d7639a 624         }
SP 625         printf("\n");
626     }
627
628 return TS_SUCCESS;
629 }
630
631 ts_bool write_vertex_fcompat_file(ts_vertex_list *vlist,ts_char *filename){
a6b1b5 632     ts_vertex **vtx=vlist->vtx;
d7639a 633     ts_uint i;
SP 634     FILE *fh;
635     
636     fh=fopen(filename, "w");
637     if(fh==NULL){
638         err("Cannot open file %s for writing");
639         return TS_FAIL;
640     }
641     for(i=0;i<vlist->n;i++)
8f6a69 642         fprintf(fh," %E\t%E\t%E\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a 643
SP 644     fclose(fh);
645 return TS_SUCCESS;
646 }
647
648
649 ts_bool fprint_vertex_list(FILE *fh,ts_vertex_list *vlist){
650     ts_uint i,j;
651     for(i=0;i<vlist->n;i++){
8f6a69 652         fprintf(fh," %.17E\t%.17E\t%.17E\t%u\n",vlist->vtx[i]->x,
SP 653             vlist->vtx[i]->y, vlist->vtx[i]->z,
654             vlist->vtx[i]->neigh_no);
655         for(j=0;j<vlist->vtx[i]->neigh_no;j++){
656             fprintf(fh,"\t%u",(ts_uint)(vlist->vtx[i]->neigh[j]->idx));
a6b1b5 657         //-vlist->vtx+1));
d7639a 658         }
SP 659         fprintf(fh,"\n");
660     }
661     return TS_SUCCESS;
662 }
663
664 ts_bool fprint_tristar(FILE *fh, ts_vesicle *vesicle){
665     ts_uint i,j;
a6b1b5 666     for(i=0;i<vesicle->vlist->n;i++){
8f6a69 667         fprintf(fh,"\t%u",vesicle->vlist->vtx[i]->tristar_no);
SP 668         for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
669             fprintf(fh,"\t%u",(ts_uint)(vesicle->vlist->vtx[i]->tristar[j]->idx));//-vesicle->tlist->tria+1));
d7639a 670         }
SP 671         fprintf(fh,"\n");
672     }
673     return TS_SUCCESS;
674 }
675
676 ts_bool fprint_triangle_list(FILE *fh, ts_vesicle *vesicle){
a6b1b5 677         ts_triangle_list *tlist=vesicle->tlist;
d7639a 678       ts_uint i,j;
SP 679     for(i=0;i<tlist->n;i++){
41a035 680         fprintf(fh,"\t%u",tlist->tria[i]->neigh_no);
SP 681         for(j=0;j<tlist->tria[i]->neigh_no;j++){
682             fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->neigh[j]->idx));//-tlist->tria+1)); 
d7639a 683         }
SP 684         fprintf(fh,"\n");
685             for(j=0;j<3;j++){
41a035 686             fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->vertex[j]->idx));//-vesicle->vlist->vtx+1)); 
d7639a 687             }
SP 688         fprintf(fh,"\n");
41a035 689         fprintf(fh,"%.17E\t%.17E\t%.17E\n",tlist->tria[i]->xnorm,
SP 690 tlist->tria[i]->ynorm,tlist->tria[i]->znorm);
d7639a 691         fprintf(fh,"0.00000000000000000\n0.00000000000000000\n");
SP 692     }
693     return TS_SUCCESS;
694 }
695
696 ts_bool fprint_vertex_data(FILE *fh,ts_vertex_list *vlist){
697     ts_uint i,j;
698     for(i=0;i<vlist->n;i++){
699         fprintf(fh," %.17E\t%.17E\t%.17E\t%.17E\t%.17E\t%u\n",
8f6a69 700         vlist->vtx[i]->xk,vlist->vtx[i]->c,vlist->vtx[i]->energy,
SP 701         vlist->vtx[i]->energy_h, vlist->vtx[i]->curvature, 0);
cef6ca 702         for(j=0;j<vlist->vtx[i]->bond_no;j++){
8f6a69 703             fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length_dual);
d7639a 704         }
SP 705             fprintf(fh,"\n");
cef6ca 706         for(j=0;j<vlist->vtx[i]->bond_no;j++){
8f6a69 707             fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length);
d7639a 708         }
SP 709             fprintf(fh,"\n");
710     }
711     return TS_SUCCESS;
712 }
713
714 ts_bool fprint_bonds(FILE *fh,ts_vesicle *vesicle){
715     ts_uint i;
a6b1b5 716     for(i=0;i<vesicle->blist->n;i++){
e016c4 717         fprintf(fh,"\t%u\t%u\n",(ts_uint)(vesicle->blist->bond[i]->vtx1->idx),
a6b1b5 718 //-vesicle->vlist->vtx+1),
e016c4 719         (ts_uint)(vesicle->blist->bond[i]->vtx2->idx));
a6b1b5 720     //-vesicle->vlist.vtx+1));
d7639a 721     }
SP 722     return TS_SUCCESS;
723 }
724
725
726 ts_bool write_dout_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
727     FILE *fh;
728     fh=fopen(filename, "w");
729     if(fh==NULL){
730         err("Cannot open file %s for writing");
731         return TS_FAIL;
732     }
733     fprintf(fh,"%.17E\n%.17E\n",vesicle->stepsize,vesicle->dmax);
a6b1b5 734     fprint_vertex_list(fh,vesicle->vlist);
d7639a 735     fprint_tristar(fh,vesicle);
SP 736     fprint_triangle_list(fh,vesicle);
a6b1b5 737     fprint_vertex_data(fh,vesicle->vlist);
d7639a 738     fprint_bonds(fh,vesicle);
SP 739     fclose(fh);    
740     return TS_SUCCESS;
741 }
742
743 ts_bool read_tape_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
744     FILE *fh;
745     char line[255];
746     fh=fopen(filename, "r");
747         if(fh==NULL){
748                 err("Cannot open file for reading... Nonexistant file?");
749                 return TS_FAIL;
750         }
a6b1b5 751     ts_uint retval=1;
d7639a 752     while(retval!=EOF){
SP 753         retval=fscanf(fh,"%s",line);
754         
755         fprintf(stderr,"%s",line);
756     }    
757     fclose(fh);    
758     return TS_SUCCESS;
759 }
760
761 ts_bool write_master_xml_file(ts_char *filename){
762      FILE *fh;
763     ts_char *i,*j;
764     ts_uint tstep;
765     ts_char *number;
766         fh=fopen(filename, "w");
767         if(fh==NULL){
768                 err("Cannot open file %s for writing");
769                 return TS_FAIL;
770         }
771
772     fprintf(fh,"<?xml version=\"1.0\"?>\n<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n<Collection>");
267db5 773     DIR *dir = opendir(command_line_args.path);
d7639a 774     if(dir){
SP 775         struct dirent *ent;
776         tstep=0;
777         while((ent = readdir(dir)) != NULL)
778         {
779             i=rindex(ent->d_name,'.');
780             if(i==NULL) continue;
781             if(strcmp(i+1,"vtu")==0){
782                     j=rindex(ent->d_name,'_');
783                     if(j==NULL) continue;
784                     number=strndup(j+1,j-i); 
a6b1b5 785                 fprintf(fh,"<DataSet timestep=\"%u\" group=\"\" part=\"0\" file=\"%s\"/>\n",atoi(number),ent->d_name);
d7639a 786                 tstep++;
SP 787                     free(number);
788             }  
789         }
790     }
f74313 791     free(dir);
d7639a 792     fprintf(fh,"</Collection>\n</VTKFile>\n");
SP 793     fclose(fh);
794     return TS_SUCCESS;
795 }
796
797 ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno){
a6b1b5 798     ts_vertex_list *vlist=vesicle->vlist;
SP 799     ts_bond_list *blist=vesicle->blist;
800     ts_vertex **vtx=vlist->vtx;
40aa5b 801     ts_uint i,j;
267db5 802         char filename[10000];
SP 803         char just_name[255];
d7639a 804     FILE *fh;
267db5 805         strcpy(filename,command_line_args.path);
SP 806         sprintf(just_name,"timestep_%.6u.vtu",timestepno);
807         strcat(filename,just_name);
d7639a 808
SP 809     fh=fopen(filename, "w");
810     if(fh==NULL){
811         err("Cannot open file %s for writing");
812         return TS_FAIL;
813     }
814     /* Here comes header of the file */
40aa5b 815
SP 816     //find number of extra vtxs and bonds of polymeres
58230a 817     ts_uint monono=0, polyno=0, poly_idx=0, filno=0, fonono=0;
M 818     ts_bool poly=0, fil=0;
40aa5b 819     if(vesicle->poly_list!=NULL){
SP 820         if(vesicle->poly_list->poly[0]!=NULL){
821         polyno=vesicle->poly_list->n;
822         monono=vesicle->poly_list->poly[0]->vlist->n;
823         poly=1;
824         }
825     }
58230a 826
M 827     if(vesicle->filament_list!=NULL){
828         if(vesicle->filament_list->poly[0]!=NULL){
829         filno=vesicle->filament_list->n;
830         fonono=vesicle->filament_list->poly[0]->vlist->n;
831         fil=1;
832         }
833     }
834
854cb6 835     fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n");
SP 836     xml_trisurf_data(fh,vesicle);
837     fprintf(fh, " <UnstructuredGrid>\n");
58230a 838     fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1));
d7639a 839     fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
SP 840        for(i=0;i<vlist->n;i++){
a6b1b5 841         fprintf(fh,"%u ",vtx[i]->idx);
d7639a 842     }
40aa5b 843     //polymeres
SP 844     if(poly){
3c1ac1 845         poly_idx=vlist->n;
40aa5b 846         for(i=0;i<vesicle->poly_list->n;i++){
3c1ac1 847             for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
58230a 848                 fprintf(fh,"%u ", poly_idx);
M 849             }
850         }
851     }
852     //filaments
853     if(fil){
854         poly_idx=vlist->n+monono*polyno;
855         for(i=0;i<vesicle->filament_list->n;i++){
856             for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
857     //    fprintf(stderr,"was here\n");
3c1ac1 858                 fprintf(fh,"%u ", poly_idx);
40aa5b 859             }
SP 860         }
861     }
d7639a 862
SP 863     fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
864     for(i=0;i<vlist->n;i++){
aad500 865         fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
40aa5b 866     }
SP 867     //polymeres
868     if(poly){
869         for(i=0;i<vesicle->poly_list->n;i++){
870             for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
aad500 871                 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 );
58230a 872             }
M 873         }
874     }
875     //filaments
876     if(fil){
877         for(i=0;i<vesicle->filament_list->n;i++){
878             for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
aad500 879                 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 );
40aa5b 880             }
SP 881         }
d7639a 882     }
SP 883
884     fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
885     for(i=0;i<blist->n;i++){
e016c4 886             fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a 887     }
40aa5b 888     //polymeres
SP 889     if(poly){
3c1ac1 890         poly_idx=vlist->n;
40aa5b 891         for(i=0;i<vesicle->poly_list->n;i++){
SP 892             for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
3c1ac1 893 //                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 894                 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 895             }
SP 896     //grafted bonds
3c1ac1 897         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 898         }
SP 899
900     }
901     
58230a 902     //filaments
M 903     if(fil){
904         poly_idx=vlist->n+monono*polyno;
905         for(i=0;i<vesicle->filament_list->n;i++){
906             for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
907                 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);
908 //        fprintf(stderr,"was here\n");
909             
910             }
911         }
912
913     }
40aa5b 914
d7639a 915     fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
58230a 916     for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2){
d7639a 917     fprintf(fh,"%u ",i);
SP 918     }
919     fprintf(fh,"\n");
920     fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
58230a 921      for (i=0;i<blist->n+monono*polyno+fonono*filno;i++){
d7639a 922         fprintf(fh,"3 ");
SP 923     }
924
925     fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
926     fclose(fh);
927     return TS_SUCCESS;
928
929 }
930
931 ts_bool write_vertex_vtk_file(ts_vesicle *vesicle,ts_char *filename, ts_char *text){
a6b1b5 932     ts_vertex_list *vlist=vesicle->vlist;
SP 933     ts_bond_list *blist=vesicle->blist;
934     ts_vertex **vtx=vlist->vtx;
935     ts_uint i;
d7639a 936     FILE *fh;
SP 937     
938     fh=fopen(filename, "w");
939     if(fh==NULL){
940         err("Cannot open file %s for writing");
941         return TS_FAIL;
942     }
943     /* Here comes header of the file */
944 //    fprintf(stderr,"NSHELL=%u\n",nshell);
945     fprintf(fh, "# vtk DataFile Version 2.0\n");
946     /* TODO: Do a sanity check on text. Max 255 char, must not me \n terminated */ 
947     fprintf(fh, "%s\n", text);
948     fprintf(fh,"ASCII\n");
949     fprintf(fh,"DATASET UNSTRUCTURED_GRID\n");
950     fprintf(fh,"POINTS %u double\n", vlist->n);
951     for(i=0;i<vlist->n;i++){
8f6a69 952         fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a 953     }
SP 954     
955     fprintf(fh,"CELLS %u %u\n",blist->n,3*blist->n);
956     for(i=0;i<blist->n;i++){
e016c4 957             fprintf(fh,"2 %u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a 958     }
SP 959     fprintf(fh,"CELL_TYPES %u\n",blist->n);
960     for(i=0;i<blist->n;i++)
961         fprintf(fh,"3\n");
962
963     fprintf(fh,"POINT_DATA %u\n", vlist->n);
964     fprintf(fh,"SCALARS scalars long 1\n");
965     fprintf(fh,"LOOKUP_TABLE default\n");
966
967     for(i=0;i<vlist->n;i++)
8f6a69 968         fprintf(fh,"%u\n",vtx[i]->idx);
d7639a 969
SP 970     fclose(fh);
971     return TS_SUCCESS;
972 }
973
974
975
144784 976 ts_bool write_pov_file(ts_vesicle *vesicle, char *filename){
SP 977     FILE *fh;
978     ts_uint i;
979     
980     fh=fopen(filename, "w");
981     if(fh==NULL){
982         err("Cannot open file %s for writing");
983         return TS_FAIL;
984     }
985
986     for(i=0;i<vesicle->tlist->n;i++){
987     
988     fprintf(fh,"\ttriangle {");
989     fprintf(fh,"\t<%e,%e,%e> <%e,%e,%e> <%e,%e,%e> }\n", 
990     vesicle->tlist->tria[i]->vertex[0]->x,
991     vesicle->tlist->tria[i]->vertex[0]->y,
992     vesicle->tlist->tria[i]->vertex[0]->z,
993
994     vesicle->tlist->tria[i]->vertex[1]->x,
995     vesicle->tlist->tria[i]->vertex[1]->y,
996     vesicle->tlist->tria[i]->vertex[1]->z,
997
998     vesicle->tlist->tria[i]->vertex[2]->x,
999     vesicle->tlist->tria[i]->vertex[2]->y,
1000     vesicle->tlist->tria[i]->vertex[2]->z
1001     );
1002     }
1003         
1004     fclose(fh);
1005     return TS_SUCCESS;
1006 }
1007
1008
1ab449 1009 ts_tape *parsetape(char *filename){
698ae1 1010     FILE *fd = fopen (filename, "r");
SP 1011     long length;
1012     size_t size;
1013     fseek (fd, 0, SEEK_END);
1014       length = ftell (fd);
1015     fseek (fd, 0, SEEK_SET);
1016     size=fread (tapetxt, 1, length, fd);
1017     fclose(fd);
1018     if(size);
1019     ts_tape *tape=parsetapebuffer(tapetxt);
1020     return tape;
1021 }
1022
1023 ts_tape *parsetapebuffer(char *buffer){
1ab449 1024     ts_tape *tape=(ts_tape *)calloc(1,sizeof(ts_tape));
SP 1025     tape->multiprocessing=calloc(255,sizeof(char));
40aa5b 1026
d7639a 1027     cfg_opt_t opts[] = {
1ab449 1028         CFG_SIMPLE_INT("nshell", &tape->nshell),
SP 1029         CFG_SIMPLE_INT("npoly", &tape->npoly),
1030         CFG_SIMPLE_INT("nmono", &tape->nmono),
58230a 1031     CFG_SIMPLE_INT("nfil",&tape->nfil),
M 1032     CFG_SIMPLE_INT("nfono",&tape->nfono),
1033     CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
37791b 1034     CFG_SIMPLE_FLOAT("R_nucleusX",&tape->R_nucleusX),
SP 1035     CFG_SIMPLE_FLOAT("R_nucleusY",&tape->R_nucleusY),
1036     CFG_SIMPLE_FLOAT("R_nucleusZ",&tape->R_nucleusZ),
58230a 1037     CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
ea1cce 1038     CFG_SIMPLE_FLOAT("dmin_interspecies", &tape->dmin_interspecies),
1ab449 1039         CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
SP 1040     CFG_SIMPLE_INT("pswitch",&tape->pswitch),
9166cb 1041     CFG_SIMPLE_INT("constvolswitch",&tape->constvolswitch),
c0ae90 1042     CFG_SIMPLE_INT("constareaswitch",&tape->constareaswitch),
43c042 1043     CFG_SIMPLE_FLOAT("constvolprecision",&tape->constvolprecision),
1ab449 1044     CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
SP 1045     CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
b30f45 1046     CFG_SIMPLE_FLOAT("xi",&tape->xi),
1ab449 1047         CFG_SIMPLE_FLOAT("stepsize",&tape->stepsize),
SP 1048         CFG_SIMPLE_INT("nxmax", &tape->ncxmax),
1049         CFG_SIMPLE_INT("nymax", &tape->ncymax),
1050         CFG_SIMPLE_INT("nzmax", &tape->nczmax),
1051         CFG_SIMPLE_INT("iterations",&tape->iterations),
1052     CFG_SIMPLE_INT("mcsweeps",&tape->mcsweeps),
1053     CFG_SIMPLE_INT("inititer", &tape->inititer),
1054         CFG_SIMPLE_BOOL("quiet",&tape->quiet),
1055         CFG_SIMPLE_STR("multiprocessing",tape->multiprocessing),
1056         CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
1057         CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
1058         CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
dc77e8 1059     CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
d7639a 1060         CFG_END()
SP 1061     };
1062     cfg_t *cfg;    
1063     ts_uint retval;
1064     cfg = cfg_init(opts, 0);
698ae1 1065     retval=cfg_parse_buf(cfg, buffer);
d7639a 1066     if(retval==CFG_FILE_ERROR){
a6b1b5 1067     fatal("No tape file.",100);
d7639a 1068     }
SP 1069     else if(retval==CFG_PARSE_ERROR){
1070     fatal("Invalid tape!",100);
1071     }
a2db52 1072
7b0c07 1073     /* here we override all values read from tape with values from commandline*/
SP 1074     getcmdline_tape(cfg,command_line_args.tape_opts);
d7639a 1075     cfg_free(cfg);
40aa5b 1076
SP 1077
1ab449 1078     /* global variables are set automatically */
SP 1079     quiet=tape->quiet;
1080     return tape;
1081 }
d7639a 1082
1ab449 1083 ts_bool tape_free(ts_tape *tape){
SP 1084     free(tape->multiprocessing);
1085     free(tape);
1086     return TS_SUCCESS;
d7639a 1087 }
f74313 1088
SP 1089
7b0c07 1090
SP 1091 ts_bool getcmdline_tape(cfg_t *cfg, char *opts){
1092
1093     char *commands, *backup, *saveptr, *saveopptr, *command, *operator[2];
1094     ts_uint i,j;
1095     commands=(char *)malloc(10000*sizeof(char));
1096     backup=commands; //since the pointer to commands will be lost, we acquire a pointer that will serve as backup.
1097     strcpy(commands,opts);
1098     for(i=0; ;i++, commands=NULL){
1099         //breaks comma separated list of commands into specific commands.
1100         command=strtok_r(commands,",",&saveptr);    
1101         if(command==NULL) break;
1102 //        fprintf(stdout,"Command %d: %s\n",i,command);    
1103         //extracts name of command and value of command into operator[2] array.
1104         for(j=0; j<2;j++,command=NULL){
1105             operator[j]=strtok_r(command,"=",&saveopptr);
1106             if(operator[j]==NULL) break;
1107 //            fprintf(stdout," ---> Operator %d: %s\n",j,operator[j]);        
1108         }
1109         //1. check: must have 2 operators.
1110         if(j!=2) fatal("Error. Command line tape options are not formatted properly",1);
1111
1112     //    cfg_setstr(cfg,operator[0],operator[1]);
1113         cmdline_to_tape(cfg,operator[0],operator[1]);
1114         //2. check: must be named properly.
1115         //3. check: must be of right format (integer, double, string, ...)
1116         
1117     }
1118     free(backup);
1119     return TS_SUCCESS;
1120 }
1121
1122
1123 ts_bool cmdline_to_tape(cfg_t *cfg, char *key, char *val){
1124
1125     cfg_opt_t *cfg_opt=cfg_getopt(cfg,key);
1126     if(cfg_opt==NULL) fatal("Commandline tape option not recognised",1); //return TS_FAIL; 
1127     switch (cfg_opt->type){
1128         case CFGT_INT:
1129             cfg_setint(cfg,key,atol(val));
1130             break;
1131         case CFGT_FLOAT:
1132             cfg_setfloat(cfg,key,atof(val));
1133             break;
1134 /*        case CFGT_BOOL:
1135             cfg_setbool(cfg,operator[0],operator[1]);
1136             break; */
1137         case CFGT_STR:
1138             cfg_setstr(cfg,key,val);
1139             break;
1140         default:
1141             break;
1142
1143     }
1144     return TS_SUCCESS;
1145 }
1146
1147
1148
f74313 1149 ts_bool read_geometry_file(char *fname, ts_vesicle *vesicle){
SP 1150     FILE *fh;
1151     ts_uint i, nvtx,nedges,ntria;
1152     ts_uint vtxi1,vtxi2;
1153     float x,y,z;
1154     ts_vertex_list *vlist;
1155     fh=fopen(fname, "r");
1156         if(fh==NULL){
1157                 err("Cannot open file for reading... Nonexistant file?");
1158                 return TS_FAIL;
1159         }
1160     ts_uint retval;
1161     retval=fscanf(fh,"%u %u %u",&nvtx, &nedges, &ntria);
1162     vesicle->vlist=init_vertex_list(nvtx);
1163     vlist=vesicle->vlist;
1164     for(i=0;i<nvtx;i++){
8f6a69 1165    //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->x,&vlist->vtx[i]->y,&vlist->vtx[i]->z);
f74313 1166        retval=fscanf(fh,"%F %F %F",&x,&y,&z);
8f6a69 1167         vlist->vtx[i]->x=x;
SP 1168         vlist->vtx[i]->y=y;
1169         vlist->vtx[i]->z=z;
f74313 1170     }
SP 1171     for(i=0;i<nedges;i++){
1172         retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
1173         bond_add(vesicle->blist,vesicle->vlist->vtx[vtxi1-1],vesicle->vlist->vtx[vtxi2-1]);
1174     }
1175     //TODO: neighbours from bonds,
1176     //TODO: triangles from neigbours
1177
1178 //    Don't need to read triangles. Already have enough data
1179     /*
1180     for(i=0;i<ntria;i++){
1181         retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
8f6a69 1182         vtxi1=vesicle->blist->vertex1->idx;
SP 1183         vtxi2=vesicle->blist->vertex1->idx;
f74313 1184         
SP 1185     }
1186     */
41a035 1187     if(retval);
f74313 1188     fclose(fh);    
SP 1189
1190
1191
1192     return TS_SUCCESS;
1193 }