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