Trisurf Monte Carlo simulator
Samo Penic
2018-12-09 514ebcb0e1b01fa41d022e3a92fe08a390a421dd
commit | author | age
22a18c 1 /* vim: set ts=4 sts=4 sw=4 noet : */
SP 2 #include<stdio.h>
3 #include<math.h>
4 #include<stdlib.h>
5 #include "general.h"
6 //#include "vertex.h"
7 //#include "bond.h"
8 //#include "triangle.h"
9 //#include "cell.h"
10 #include "vesicle.h"
11 #include "io.h"
12 //#include "initial_distribution.h"
13 //#include "frame.h"
14 //#include "timestep.h"
15 //#include "poly.h"
16 #include "stats.h"
17 #include "sh.h"
18 #include "shcomplex.h"
19 #include "dumpstate.h"
20 #include "restore.h"
8c91d2 21 #include "cluster.h"
22a18c 22 #include <string.h>
SP 23 #include <getopt.h>
24 #include <sys/stat.h>
25 #include <sys/types.h>
26 #include <dirent.h>
27 #include <errno.h>
28 #include <snapshot.h>
29 #include<gsl/gsl_complex.h>
30 #include<gsl/gsl_complex_math.h>
033407 31 #include<stdio.h>
22a18c 32
SP 33 ts_vesicle *restoreVesicle(char *filename){
34     ts_vesicle *vesicle = parseDump(filename);
35     return vesicle;
36 }
37
38 void vesicle_calculate_ulm2(ts_vesicle *vesicle){
39     //complex_sph_free(vesicle->sphHarmonics);
40
41     //vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,21);
42     vesicle_volume(vesicle);
43     preparationSh(vesicle,getR0(vesicle));
44     calculateUlmComplex(vesicle);
45     ts_int i,j;
46     for(i=0;i<vesicle->sphHarmonics->l;i++){
47             for(j=i;j<2*i+1;j++){
48             printf("%e ", gsl_complex_abs2(vesicle->sphHarmonics->ulmComplex[i][j]));
49             }
50     }
51         printf("\n");
52
53 }
54
89434d 55
SP 56
57 int count_bonds_with_energy(ts_bond_list *blist){
58
59     unsigned int i, cnt;
60     cnt=0;
61     for(i=0;i<blist->n;i++){
62         if(fabs(blist->bond[i]->energy)>1e-16) cnt++;
63     }
64     return cnt;
65 }
033407 66
SP 67
68
69 ts_bool write_histogram_data(ts_uint timestep_no, ts_vesicle *vesicle){
70     ts_cluster_list *cstlist=init_cluster_list();
71     clusterize_vesicle(vesicle,cstlist);
72     //printf("No clusters=%d\n",cstlist->n);
0a2c81 73     int k,i,cnt, test=0;
033407 74     int max_nvtx=0;
SP 75     char filename[255];
76     sprintf(filename,"histogram_%.6u.csv",timestep_no);
77     FILE *fd=fopen(filename,"w");
78     fprintf(fd,"Number_of_vertices_in cluster Number_of_clusters\n");
79     for(k=0;k<cstlist->n;k++)
80         if(cstlist->cluster[k]->nvtx>max_nvtx) max_nvtx=cstlist->cluster[k]->nvtx;
81     //printf("Max. number of vertices in cluster: %d\n",max_nvtx);
82     for(i=1;i<=max_nvtx;i++){
83         cnt=0;
84         for(k=0;k<cstlist->n;k++)
85             if(cstlist->cluster[k]->nvtx==i) cnt++;
86         fprintf(fd,"%d %d\n",i,cnt);
0a2c81 87         test+=cnt*i;
033407 88     }
SP 89     //for(k=0;k<cstlist->n;k++){
90 //        printf("*Cluster %d has %d vertices\n",k,cstlist->cluster[k]->nvtx);
91 //    }
92
93     fclose(fd);
0a2c81 94 //    printf("*Sum of all vertices in clusters: %d\n", test);
SP 95 //    write_vertex_xml_file(vesicle,timestep_no,cstlist);
033407 96     cluster_list_free(cstlist);
SP 97     
98     return TS_SUCCESS;
99 }
100
89434d 101
22a18c 102 int main(){
SP 103     ts_vesicle *vesicle;
104     ts_char *i,*j;
105     ts_uint tstep,n;
106         ts_char *number;
107     struct dirent **list;
810894 108     ts_double l1,l2,l3,hbar;
22a18c 109     int count;
SP 110     ts_fprintf(stderr,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);
111
810894 112     fprintf(stdout, "OuterLoop Volume Area lamdba1 lambda2 lambda3 Nbw/Nb hbar\n");
22a18c 113
SP 114
115     count=scandir(".",&list,0,alphasort);
116     if(count<0){
117         fatal("Error, cannot open directory.",1);
118     }
119         tstep=0;
120     for(n=0;n<count;n++){
121         struct dirent *ent;
122         ent=list[n];    
123                 i=rindex(ent->d_name,'.');
124                 if(i==NULL) {
125                 continue;
126         }
127                 if(strcmp(i+1,"vtu")==0){
128                     j=rindex(ent->d_name,'_');
129                     if(j==NULL) continue;
130                     number=strndup(j+1,j-i); 
131             quiet=1;
132                     ts_fprintf(stdout,"timestep: %u filename: %s\n",atoi(number),ent->d_name);
133 //            printf("%u ",atoi(number));
134             vesicle=restoreVesicle(ent->d_name);
135 //            vesicle_calculate_ulm2(vesicle);
136             vesicle_volume(vesicle);
137             vesicle_area(vesicle);
138             gyration_eigen(vesicle,&l1,&l2,&l3);
810894 139             hbar=vesicle_meancurvature(vesicle)/vesicle->area;            
M 140             fprintf(stdout,"%d %.17e %.17e %.17e %.17e %.17e %.17e %.17e\n",atoi(number),vesicle->volume, vesicle->area,l1,l2,l3, (ts_double)count_bonds_with_energy(vesicle->blist)/(ts_double)vesicle->blist->n,hbar);
22a18c 141                         tstep++;
033407 142             write_histogram_data(atoi(number), vesicle);
22a18c 143                     free(number);
SP 144             tape_free(vesicle->tape);
145             vesicle_free(vesicle);
146                 }
147         }
148     for (n = 0; n < count; n++)
149       {
150           free(list[n]);
151       }
152     
153     free(list);
154     return 0;
155 }
156