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