Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2013-11-27 0652ee2715994b48a9bbfd1132bf0b65ad206289
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"
d7639a 12
SP 13 ts_bool print_vertex_list(ts_vertex_list *vlist){
14     ts_uint i;
15     printf("Number of vertices: %u\n",vlist->n);
16     for(i=0;i<vlist->n;i++){
a6b1b5 17         printf("%u: %f %f %f\n",
8f6a69 18 vlist->vtx[i]->idx,vlist->vtx[i]->x, vlist->vtx[i]->y, vlist->vtx[i]->z);
d7639a 19     }
SP 20     return TS_SUCCESS;
21 }
22
23 ts_bool print_vertex_neighbours(ts_vertex_list *vlist){
24     ts_uint i,j;
a6b1b5 25     ts_vertex **vtx=vlist->vtx;
d7639a 26     printf("Vertex id(neigh no): (neighvertex coord) (neighvertex coord) ...\n");
SP 27     for(i=0;i<vlist->n;i++){
8f6a69 28         printf("%u(%u): ",vtx[i]->idx,vtx[i]->neigh_no);
SP 29         for(j=0;j<vtx[i]->neigh_no;j++){
30             printf("(%f,%f,%f)",vtx[i]->neigh[j]->x,
31 vtx[i]->neigh[j]->y,vtx[i]->neigh[j]->z);
d7639a 32         }
SP 33         printf("\n");
34     }
35
36 return TS_SUCCESS;
37 }
38
39 ts_bool write_vertex_fcompat_file(ts_vertex_list *vlist,ts_char *filename){
a6b1b5 40     ts_vertex **vtx=vlist->vtx;
d7639a 41     ts_uint i;
SP 42     FILE *fh;
43     
44     fh=fopen(filename, "w");
45     if(fh==NULL){
46         err("Cannot open file %s for writing");
47         return TS_FAIL;
48     }
49     for(i=0;i<vlist->n;i++)
8f6a69 50         fprintf(fh," %E\t%E\t%E\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a 51
SP 52     fclose(fh);
53 return TS_SUCCESS;
54 }
55
56
57 ts_bool fprint_vertex_list(FILE *fh,ts_vertex_list *vlist){
58     ts_uint i,j;
59     for(i=0;i<vlist->n;i++){
8f6a69 60         fprintf(fh," %.17E\t%.17E\t%.17E\t%u\n",vlist->vtx[i]->x,
SP 61             vlist->vtx[i]->y, vlist->vtx[i]->z,
62             vlist->vtx[i]->neigh_no);
63         for(j=0;j<vlist->vtx[i]->neigh_no;j++){
64             fprintf(fh,"\t%u",(ts_uint)(vlist->vtx[i]->neigh[j]->idx));
a6b1b5 65         //-vlist->vtx+1));
d7639a 66         }
SP 67         fprintf(fh,"\n");
68     }
69     return TS_SUCCESS;
70 }
71
72 ts_bool fprint_tristar(FILE *fh, ts_vesicle *vesicle){
73     ts_uint i,j;
a6b1b5 74     for(i=0;i<vesicle->vlist->n;i++){
8f6a69 75         fprintf(fh,"\t%u",vesicle->vlist->vtx[i]->tristar_no);
SP 76         for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
77             fprintf(fh,"\t%u",(ts_uint)(vesicle->vlist->vtx[i]->tristar[j]->idx));//-vesicle->tlist->tria+1));
d7639a 78         }
SP 79         fprintf(fh,"\n");
80     }
81     return TS_SUCCESS;
82 }
83
84 ts_bool fprint_triangle_list(FILE *fh, ts_vesicle *vesicle){
a6b1b5 85         ts_triangle_list *tlist=vesicle->tlist;
d7639a 86       ts_uint i,j;
SP 87     for(i=0;i<tlist->n;i++){
41a035 88         fprintf(fh,"\t%u",tlist->tria[i]->neigh_no);
SP 89         for(j=0;j<tlist->tria[i]->neigh_no;j++){
90             fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->neigh[j]->idx));//-tlist->tria+1)); 
d7639a 91         }
SP 92         fprintf(fh,"\n");
93             for(j=0;j<3;j++){
41a035 94             fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->vertex[j]->idx));//-vesicle->vlist->vtx+1)); 
d7639a 95             }
SP 96         fprintf(fh,"\n");
41a035 97         fprintf(fh,"%.17E\t%.17E\t%.17E\n",tlist->tria[i]->xnorm,
SP 98 tlist->tria[i]->ynorm,tlist->tria[i]->znorm);
d7639a 99         fprintf(fh,"0.00000000000000000\n0.00000000000000000\n");
SP 100     }
101     return TS_SUCCESS;
102 }
103
104 ts_bool fprint_vertex_data(FILE *fh,ts_vertex_list *vlist){
105     ts_uint i,j;
106     for(i=0;i<vlist->n;i++){
107         fprintf(fh," %.17E\t%.17E\t%.17E\t%.17E\t%.17E\t%u\n",
8f6a69 108         vlist->vtx[i]->xk,vlist->vtx[i]->c,vlist->vtx[i]->energy,
SP 109         vlist->vtx[i]->energy_h, vlist->vtx[i]->curvature, 0);
110         for(j=0;j<vlist->vtx[i]->neigh_no;j++){
111             fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length_dual);
d7639a 112         }
SP 113             fprintf(fh,"\n");
8f6a69 114         for(j=0;j<vlist->vtx[i]->neigh_no;j++){
SP 115             fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length);
d7639a 116         }
SP 117             fprintf(fh,"\n");
118     }
119     return TS_SUCCESS;
120 }
121
122 ts_bool fprint_bonds(FILE *fh,ts_vesicle *vesicle){
123     ts_uint i;
a6b1b5 124     for(i=0;i<vesicle->blist->n;i++){
e016c4 125         fprintf(fh,"\t%u\t%u\n",(ts_uint)(vesicle->blist->bond[i]->vtx1->idx),
a6b1b5 126 //-vesicle->vlist->vtx+1),
e016c4 127         (ts_uint)(vesicle->blist->bond[i]->vtx2->idx));
a6b1b5 128     //-vesicle->vlist.vtx+1));
d7639a 129     }
SP 130     return TS_SUCCESS;
131 }
132
133
134 ts_bool write_dout_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
135     FILE *fh;
136     fh=fopen(filename, "w");
137     if(fh==NULL){
138         err("Cannot open file %s for writing");
139         return TS_FAIL;
140     }
141     fprintf(fh,"%.17E\n%.17E\n",vesicle->stepsize,vesicle->dmax);
a6b1b5 142     fprint_vertex_list(fh,vesicle->vlist);
d7639a 143     fprint_tristar(fh,vesicle);
SP 144     fprint_triangle_list(fh,vesicle);
a6b1b5 145     fprint_vertex_data(fh,vesicle->vlist);
d7639a 146     fprint_bonds(fh,vesicle);
SP 147     fclose(fh);    
148     return TS_SUCCESS;
149 }
150
151 ts_bool read_tape_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
152     FILE *fh;
153     char line[255];
154     fh=fopen(filename, "r");
155         if(fh==NULL){
156                 err("Cannot open file for reading... Nonexistant file?");
157                 return TS_FAIL;
158         }
a6b1b5 159     ts_uint retval=1;
d7639a 160     while(retval!=EOF){
SP 161         retval=fscanf(fh,"%s",line);
162         
163         fprintf(stderr,"%s",line);
164     }    
165     fclose(fh);    
166     return TS_SUCCESS;
167 }
168
169 ts_bool write_master_xml_file(ts_char *filename){
170      FILE *fh;
171     ts_char *i,*j;
172     ts_uint tstep;
173     ts_char *number;
174         fh=fopen(filename, "w");
175         if(fh==NULL){
176                 err("Cannot open file %s for writing");
177                 return TS_FAIL;
178         }
179
180     fprintf(fh,"<?xml version=\"1.0\"?>\n<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n<Collection>");
181     DIR *dir = opendir(".");
182     if(dir){
183         struct dirent *ent;
184         tstep=0;
185         while((ent = readdir(dir)) != NULL)
186         {
187             i=rindex(ent->d_name,'.');
188             if(i==NULL) continue;
189             if(strcmp(i+1,"vtu")==0){
190                     j=rindex(ent->d_name,'_');
191                     if(j==NULL) continue;
192                     number=strndup(j+1,j-i); 
a6b1b5 193                 fprintf(fh,"<DataSet timestep=\"%u\" group=\"\" part=\"0\" file=\"%s\"/>\n",atoi(number),ent->d_name);
d7639a 194                 tstep++;
SP 195                     free(number);
196             }  
197         }
198     }
f74313 199     free(dir);
d7639a 200     fprintf(fh,"</Collection>\n</VTKFile>\n");
SP 201     fclose(fh);
202     return TS_SUCCESS;
203 }
204
205 ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno){
a6b1b5 206     ts_vertex_list *vlist=vesicle->vlist;
SP 207     ts_bond_list *blist=vesicle->blist;
208     ts_vertex **vtx=vlist->vtx;
209     ts_uint i;
d7639a 210         char filename[255];
SP 211     FILE *fh;
212
213         sprintf(filename,"timestep_%.6u.vtu",timestepno);
214     fh=fopen(filename, "w");
215     if(fh==NULL){
216         err("Cannot open file %s for writing");
217         return TS_FAIL;
218     }
219     /* Here comes header of the file */
220     fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n <UnstructuredGrid>\n");
221     fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n, blist->n);
222     fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
223        for(i=0;i<vlist->n;i++){
a6b1b5 224         fprintf(fh,"%u ",vtx[i]->idx);
d7639a 225     }
SP 226
227     fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
228     for(i=0;i<vlist->n;i++){
8f6a69 229         fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a 230     }
SP 231
232     fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
233     for(i=0;i<blist->n;i++){
e016c4 234             fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a 235     }
SP 236     fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
237     for (i=2;i<blist->n*2+1;i+=2){
238     fprintf(fh,"%u ",i);
239     }
240     fprintf(fh,"\n");
241     fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
242      for (i=0;i<blist->n;i++){
243         fprintf(fh,"3 ");
244     }
245
246     fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
247     fclose(fh);
248     return TS_SUCCESS;
249
250 }
251
252 ts_bool write_vertex_vtk_file(ts_vesicle *vesicle,ts_char *filename, ts_char *text){
a6b1b5 253     ts_vertex_list *vlist=vesicle->vlist;
SP 254     ts_bond_list *blist=vesicle->blist;
255     ts_vertex **vtx=vlist->vtx;
256     ts_uint i;
d7639a 257     FILE *fh;
SP 258     
259     fh=fopen(filename, "w");
260     if(fh==NULL){
261         err("Cannot open file %s for writing");
262         return TS_FAIL;
263     }
264     /* Here comes header of the file */
265 //    fprintf(stderr,"NSHELL=%u\n",nshell);
266     fprintf(fh, "# vtk DataFile Version 2.0\n");
267     /* TODO: Do a sanity check on text. Max 255 char, must not me \n terminated */ 
268     fprintf(fh, "%s\n", text);
269     fprintf(fh,"ASCII\n");
270     fprintf(fh,"DATASET UNSTRUCTURED_GRID\n");
271     fprintf(fh,"POINTS %u double\n", vlist->n);
272     for(i=0;i<vlist->n;i++){
8f6a69 273         fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a 274     }
SP 275     
276     fprintf(fh,"CELLS %u %u\n",blist->n,3*blist->n);
277     for(i=0;i<blist->n;i++){
e016c4 278             fprintf(fh,"2 %u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a 279     }
SP 280     fprintf(fh,"CELL_TYPES %u\n",blist->n);
281     for(i=0;i<blist->n;i++)
282         fprintf(fh,"3\n");
283
284     fprintf(fh,"POINT_DATA %u\n", vlist->n);
285     fprintf(fh,"SCALARS scalars long 1\n");
286     fprintf(fh,"LOOKUP_TABLE default\n");
287
288     for(i=0;i<vlist->n;i++)
8f6a69 289         fprintf(fh,"%u\n",vtx[i]->idx);
d7639a 290
SP 291     fclose(fh);
292     return TS_SUCCESS;
293 }
294
295
296
d7a113 297 ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations){
d7639a 298     long int nshell=17,ncxmax=60, ncymax=60, nczmax=60;  // THIS IS DUE TO CONFUSE BUG!
b14a8d 299     char *buf=malloc(255*sizeof(char));
50f10f 300     buf[0]=0;
SP 301     long int smps=1;
b14a8d 302     long int brezveze1=1;
SP 303     long int brezveze2=1;
d7639a 304     ts_double xk0=25.0, dmax=1.67,stepsize=0.15;
d7a113 305     long int iter=1000, init=1000, mcsw=1000;
d7639a 306     cfg_opt_t opts[] = {
SP 307         CFG_SIMPLE_INT("nshell", &nshell),
308         CFG_SIMPLE_FLOAT("dmax", &dmax),
309         CFG_SIMPLE_FLOAT("xk0",&xk0),
310         CFG_SIMPLE_FLOAT("stepsize",&stepsize),
311         CFG_SIMPLE_INT("nxmax", &ncxmax),
312         CFG_SIMPLE_INT("nymax", &ncymax),
313         CFG_SIMPLE_INT("nzmax", &nczmax),
d7a113 314         CFG_SIMPLE_INT("iterations",&iter),
SP 315     CFG_SIMPLE_INT("mcsweeps",&mcsw),
316     CFG_SIMPLE_INT("inititer", &init),
d7639a 317         CFG_SIMPLE_BOOL("quiet",&quiet),
50f10f 318         CFG_SIMPLE_STR("multiprocessing",&buf),
SP 319         CFG_SIMPLE_INT("smp_cores",&smps),
b14a8d 320         CFG_SIMPLE_INT("cluster_nodes",&brezveze1),
SP 321         CFG_SIMPLE_INT("distributed_processes",&brezveze2),
d7639a 322         CFG_END()
SP 323     };
324     cfg_t *cfg;    
325     ts_uint retval;
326     cfg = cfg_init(opts, 0);
a6b1b5 327     retval=cfg_parse(cfg, "tape");
d7639a 328     if(retval==CFG_FILE_ERROR){
a6b1b5 329     fatal("No tape file.",100);
d7639a 330     }
SP 331     else if(retval==CFG_PARSE_ERROR){
332     fatal("Invalid tape!",100);
333     }
b14a8d 334     ts_vesicle *vesicle;
d7a113 335         *iterations=iter;
SP 336     *inititer=init;
337     *mcsweeps=mcsw;
b14a8d 338     vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize);
d7639a 339     vesicle->nshell=nshell;
SP 340     vesicle->dmax=dmax*dmax;
341     vesicle->bending_rigidity=xk0;
342     vesicle->stepsize=stepsize;
a6b1b5 343     vesicle->clist->ncmax[0]=ncxmax;
SP 344     vesicle->clist->ncmax[1]=ncymax;
345     vesicle->clist->ncmax[2]=nczmax;
346     vesicle->clist->max_occupancy=8;
50f10f 347     if(strcmp(buf,"smp")==0){    
SP 348     vesicle->threads=smps;
349     } else {
350     vesicle->threads=1;
351     }
d7639a 352     cfg_free(cfg);
b14a8d 353     free(buf);
SP 354   //  fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
355     return vesicle;
d7639a 356
SP 357 }
f74313 358
SP 359
360 ts_bool read_geometry_file(char *fname, ts_vesicle *vesicle){
361     FILE *fh;
362     ts_uint i, nvtx,nedges,ntria;
363     ts_uint vtxi1,vtxi2;
364     float x,y,z;
365     ts_vertex_list *vlist;
366     fh=fopen(fname, "r");
367         if(fh==NULL){
368                 err("Cannot open file for reading... Nonexistant file?");
369                 return TS_FAIL;
370         }
371     ts_uint retval;
372     retval=fscanf(fh,"%u %u %u",&nvtx, &nedges, &ntria);
373     vesicle->vlist=init_vertex_list(nvtx);
374     vlist=vesicle->vlist;
375     for(i=0;i<nvtx;i++){
8f6a69 376    //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->x,&vlist->vtx[i]->y,&vlist->vtx[i]->z);
f74313 377        retval=fscanf(fh,"%F %F %F",&x,&y,&z);
8f6a69 378         vlist->vtx[i]->x=x;
SP 379         vlist->vtx[i]->y=y;
380         vlist->vtx[i]->z=z;
f74313 381     }
SP 382     for(i=0;i<nedges;i++){
383         retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
384         bond_add(vesicle->blist,vesicle->vlist->vtx[vtxi1-1],vesicle->vlist->vtx[vtxi2-1]);
385     }
386     //TODO: neighbours from bonds,
387     //TODO: triangles from neigbours
388
389 //    Don't need to read triangles. Already have enough data
390     /*
391     for(i=0;i<ntria;i++){
392         retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
8f6a69 393         vtxi1=vesicle->blist->vertex1->idx;
SP 394         vtxi2=vesicle->blist->vertex1->idx;
f74313 395         
SP 396     }
397     */
41a035 398     if(retval);
f74313 399     fclose(fh);    
SP 400
401
402
403     return TS_SUCCESS;
404 }