Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2014-02-11 6aec63272fd6fe35a6ded508c2f2e3e19e107215
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>
11
12
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
297 ts_bool parsetape(ts_vesicle *vesicle,ts_uint *iterations){
298     long int nshell=17,ncxmax=60, ncymax=60, nczmax=60;  // THIS IS DUE TO CONFUSE BUG!
314f2d 299     char buf[255];
SP 300     long int brezveze=1;
d7639a 301     ts_double xk0=25.0, dmax=1.67,stepsize=0.15;
SP 302     *iterations=1000;
303     cfg_opt_t opts[] = {
304         CFG_SIMPLE_INT("nshell", &nshell),
305         CFG_SIMPLE_FLOAT("dmax", &dmax),
306         CFG_SIMPLE_FLOAT("xk0",&xk0),
307         CFG_SIMPLE_FLOAT("stepsize",&stepsize),
308         CFG_SIMPLE_INT("nxmax", &ncxmax),
309         CFG_SIMPLE_INT("nymax", &ncymax),
310         CFG_SIMPLE_INT("nzmax", &nczmax),
311         CFG_SIMPLE_INT("iterations",iterations),
312         CFG_SIMPLE_BOOL("quiet",&quiet),
314f2d 313         CFG_SIMPLE_STR("multiprocessing",buf),
SP 314         CFG_SIMPLE_INT("smp_cores",&brezveze),
315         CFG_SIMPLE_INT("cluster_nodes",&brezveze),
316         CFG_SIMPLE_INT("distributed_processes",&brezveze),
d7639a 317         CFG_END()
SP 318     };
319     cfg_t *cfg;    
320     ts_uint retval;
321     cfg = cfg_init(opts, 0);
a6b1b5 322     retval=cfg_parse(cfg, "tape");
d7639a 323     if(retval==CFG_FILE_ERROR){
a6b1b5 324     fatal("No tape file.",100);
d7639a 325     }
SP 326     else if(retval==CFG_PARSE_ERROR){
327     fatal("Invalid tape!",100);
328     }
329     vesicle->nshell=nshell;
330     vesicle->dmax=dmax*dmax;
331     vesicle->bending_rigidity=xk0;
332     vesicle->stepsize=stepsize;
a6b1b5 333     vesicle->clist->ncmax[0]=ncxmax;
SP 334     vesicle->clist->ncmax[1]=ncymax;
335     vesicle->clist->ncmax[2]=nczmax;
336     vesicle->clist->max_occupancy=8;
d7639a 337     cfg_free(cfg);
SP 338 //    fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
339     return TS_SUCCESS;
340
341 }
f74313 342
SP 343
344 ts_bool read_geometry_file(char *fname, ts_vesicle *vesicle){
345     FILE *fh;
346     ts_uint i, nvtx,nedges,ntria;
347     ts_uint vtxi1,vtxi2;
348     float x,y,z;
349     ts_vertex_list *vlist;
350     fh=fopen(fname, "r");
351         if(fh==NULL){
352                 err("Cannot open file for reading... Nonexistant file?");
353                 return TS_FAIL;
354         }
355     ts_uint retval;
356     retval=fscanf(fh,"%u %u %u",&nvtx, &nedges, &ntria);
357     vesicle->vlist=init_vertex_list(nvtx);
358     vlist=vesicle->vlist;
359     for(i=0;i<nvtx;i++){
8f6a69 360    //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->x,&vlist->vtx[i]->y,&vlist->vtx[i]->z);
f74313 361        retval=fscanf(fh,"%F %F %F",&x,&y,&z);
8f6a69 362         vlist->vtx[i]->x=x;
SP 363         vlist->vtx[i]->y=y;
364         vlist->vtx[i]->z=z;
f74313 365     }
SP 366     for(i=0;i<nedges;i++){
367         retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
368         bond_add(vesicle->blist,vesicle->vlist->vtx[vtxi1-1],vesicle->vlist->vtx[vtxi2-1]);
369     }
370     //TODO: neighbours from bonds,
371     //TODO: triangles from neigbours
372
373 //    Don't need to read triangles. Already have enough data
374     /*
375     for(i=0;i<ntria;i++){
376         retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
8f6a69 377         vtxi1=vesicle->blist->vertex1->idx;
SP 378         vtxi2=vesicle->blist->vertex1->idx;
f74313 379         
SP 380     }
381     */
41a035 382     if(retval);
f74313 383     fclose(fh);    
SP 384
385
386
387     return TS_SUCCESS;
388 }