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