Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2012-02-23 34d3de02dc183d38bf661789e98e43d9c1b4c0bd
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",
SP 18 vlist->vtx[i]->idx,vlist->vtx[i]->data->x, vlist->vtx[i]->data->y, vlist->vtx[i]->data->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++){
a6b1b5 28         printf("%u(%u): ",vtx[i]->idx,vtx[i]->data->neigh_no);
SP 29         for(j=0;j<vtx[i]->data->neigh_no;j++){
30             printf("(%f,%f,%f)",vtx[i]->data->neigh[j]->data->x,
31 vtx[i]->data->neigh[j]->data->y,vtx[i]->data->neigh[j]->data->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++)
a6b1b5 50         fprintf(fh," %E\t%E\t%E\n",vtx[i]->data->x,vtx[i]->data->y, vtx[i]->data->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++){
a6b1b5 60         fprintf(fh," %.17E\t%.17E\t%.17E\t%u\n",vlist->vtx[i]->data->x,
SP 61             vlist->vtx[i]->data->y, vlist->vtx[i]->data->z,
62             vlist->vtx[i]->data->neigh_no);
63         for(j=0;j<vlist->vtx[i]->data->neigh_no;j++){
64             fprintf(fh,"\t%u",(ts_uint)(vlist->vtx[i]->data->neigh[j]->idx));
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++){
SP 75         fprintf(fh,"\t%u",vesicle->vlist->vtx[i]->data->tristar_no);
76         for(j=0;j<vesicle->vlist->vtx[i]->data->tristar_no;j++){
77             fprintf(fh,"\t%u",(ts_uint)(vesicle->vlist->vtx[i]->data->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",
a6b1b5 108         vlist->vtx[i]->data->xk,vlist->vtx[i]->data->c,vlist->vtx[i]->data->energy,
SP 109         vlist->vtx[i]->data->energy_h, vlist->vtx[i]->data->curvature, 0);
110         for(j=0;j<vlist->vtx[i]->data->neigh_no;j++){
e016c4 111             fprintf(fh," %.17E", vlist->vtx[i]->data->bond[j]->bond_length_dual);
d7639a 112         }
SP 113             fprintf(fh,"\n");
a6b1b5 114         for(j=0;j<vlist->vtx[i]->data->neigh_no;j++){
e016c4 115             fprintf(fh," %.17E", vlist->vtx[i]->data->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++){
a6b1b5 229         fprintf(fh,"%e %e %e\n",vtx[i]->data->x,vtx[i]->data->y, vtx[i]->data->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++){
a6b1b5 273         fprintf(fh,"%e %e %e\n",vtx[i]->data->x,vtx[i]->data->y, vtx[i]->data->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++)
a6b1b5 289         fprintf(fh,"%u\n",vtx[i]->data->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!
299     ts_double xk0=25.0, dmax=1.67,stepsize=0.15;
300     *iterations=1000;
301     cfg_opt_t opts[] = {
302         CFG_SIMPLE_INT("nshell", &nshell),
303         CFG_SIMPLE_FLOAT("dmax", &dmax),
304         CFG_SIMPLE_FLOAT("xk0",&xk0),
305         CFG_SIMPLE_FLOAT("stepsize",&stepsize),
306         CFG_SIMPLE_INT("nxmax", &ncxmax),
307         CFG_SIMPLE_INT("nymax", &ncymax),
308         CFG_SIMPLE_INT("nzmax", &nczmax),
309         CFG_SIMPLE_INT("iterations",iterations),
310         CFG_SIMPLE_BOOL("quiet",&quiet),
311         CFG_END()
312     };
313     cfg_t *cfg;    
314     ts_uint retval;
315     cfg = cfg_init(opts, 0);
a6b1b5 316     retval=cfg_parse(cfg, "tape");
d7639a 317     if(retval==CFG_FILE_ERROR){
a6b1b5 318     fatal("No tape file.",100);
d7639a 319     }
SP 320     else if(retval==CFG_PARSE_ERROR){
321     fatal("Invalid tape!",100);
322     }
323     vesicle->nshell=nshell;
324     vesicle->dmax=dmax*dmax;
325     vesicle->bending_rigidity=xk0;
326     vesicle->stepsize=stepsize;
a6b1b5 327     vesicle->clist->ncmax[0]=ncxmax;
SP 328     vesicle->clist->ncmax[1]=ncymax;
329     vesicle->clist->ncmax[2]=nczmax;
330     vesicle->clist->max_occupancy=8;
d7639a 331     cfg_free(cfg);
SP 332 //    fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
333     return TS_SUCCESS;
334
335 }
f74313 336
SP 337
338 ts_bool read_geometry_file(char *fname, ts_vesicle *vesicle){
339     FILE *fh;
340     ts_uint i, nvtx,nedges,ntria;
341     ts_uint vtxi1,vtxi2;
342     float x,y,z;
343     ts_vertex_list *vlist;
344     fh=fopen(fname, "r");
345         if(fh==NULL){
346                 err("Cannot open file for reading... Nonexistant file?");
347                 return TS_FAIL;
348         }
349     ts_uint retval;
350     retval=fscanf(fh,"%u %u %u",&nvtx, &nedges, &ntria);
351     vesicle->vlist=init_vertex_list(nvtx);
352     vlist=vesicle->vlist;
353     for(i=0;i<nvtx;i++){
354    //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->data->x,&vlist->vtx[i]->data->y,&vlist->vtx[i]->data->z);
355        retval=fscanf(fh,"%F %F %F",&x,&y,&z);
356         vlist->vtx[i]->data->x=x;
357         vlist->vtx[i]->data->y=y;
358         vlist->vtx[i]->data->z=z;
359     }
360     for(i=0;i<nedges;i++){
361         retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
362         bond_add(vesicle->blist,vesicle->vlist->vtx[vtxi1-1],vesicle->vlist->vtx[vtxi2-1]);
363     }
364     //TODO: neighbours from bonds,
365     //TODO: triangles from neigbours
366
367 //    Don't need to read triangles. Already have enough data
368     /*
369     for(i=0;i<ntria;i++){
370         retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
371         vtxi1=vesicle->blist->data->vertex1->idx;
372         vtxi2=vesicle->blist->data->vertex1->idx;
373         
374     }
375     */
41a035 376     if(retval);
f74313 377     fclose(fh);    
SP 378
379
380
381     return TS_SUCCESS;
382 }