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