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