Trisurf Monte Carlo simulator
Samo Penic
2012-07-10 e19e790e95f14ca69a7ce9c5e45d815fe21df36e
commit | author | age
7958e9 1 #include<stdlib.h>
SP 2 #include<math.h>
3 #include<stdio.h>
4 #include "general.h"
5 #include "vertex.h"
6 #include "bond.h"
7 #include "vesicle.h"
8 #include "vertex.h"
9 #include "triangle.h"
10 #include "initial_distribution.h"
f74313 11 #include "energy.h"
7958e9 12
SP 13 ts_vesicle *initial_distribution_dipyramid(ts_uint nshell, ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
14     ts_fprintf(stderr,"Starting initial_distribution on vesicle with %u shells!...\n",nshell);
15     ts_bool retval;
16     ts_uint no_vertices=5*nshell*nshell+2;
314f2d 17
SP 18
7958e9 19     
SP 20     ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize);
314f2d 21
SP 22 //TODO: debugging only. Please remove ASAP!
83b03a 23     vesicle->bending_rigidity=25.0*25.0;
314f2d 24
7958e9 25     vesicle->nshell=nshell;
SP 26     retval = vtx_set_global_values(vesicle);
27     retval = pentagonal_dipyramid_vertex_distribution(vesicle->vlist);
28     retval = init_vertex_neighbours(vesicle->vlist);
b01cc1 29     vesicle->vlist = init_sort_neighbours(vesicle->blist,vesicle->vlist);
SP 30    // retval = init_vesicle_bonds(vesicle); // bonds are created in sort_neigh
7958e9 31     retval = init_triangles(vesicle);
SP 32     retval = init_triangle_neighbours(vesicle);
33     retval = init_common_vertex_triangle_neighbours(vesicle);
dac2e5 34     retval = init_normal_vectors(vesicle->tlist);
f74313 35     retval = mean_curvature_and_energy(vesicle);
7958e9 36  ts_fprintf(stderr,"initial_distribution finished!\n");
41a035 37     if(retval);
7958e9 38     return vesicle;
SP 39
40
41
42 ts_bool pentagonal_dipyramid_vertex_distribution(ts_vertex_list *vlist){
43     /* Some often used relations */
44     const ts_double s1= sin(2.0*M_PI/5.0);
45     const ts_double s2= sin(4.0*M_PI/5.0);
46     const ts_double c1= cos(2.0*M_PI/5.0);
47     const ts_double c2= cos(4.0*M_PI/5.0);
48
49     /* Calculates projection lenght of an edge bond to pentagram plane */
50     const ts_double xl0=A0/(2.0*sin(M_PI/5.0));
51 #ifdef TS_DOUBLE_DOUBLE
52     const ts_double z0=sqrt(pow(A0,2)-pow(xl0,2));
53 #endif
54 #ifdef TS_DOUBLE_FLOAT
55     const ts_double z0=sqrtf(powf(A0,2)-powf(xl0,2));
56 #endif
57 #ifdef TS_DOUBLE_LONGDOUBLE
58     const ts_double z0=sqrtl(powl(A0,2)-powl(xl0,2));
59 #endif
60 //    const z0=sqrt(A0*A0 -xl0*xl0); /* I could use pow function but if pow is used make a check on the float type. If float then powf, if long double use powl */
61
62 /*placeholder for the pointer to vertex datastructure list... DIRTY: actual pointer points towards invalid address, one position before actual beginning of the list... This is to solve the difference between 1 based indexing in original program in fortran and 0 based indexing in C. All algorithms remain unchanged because of this!*/
63     ts_vertex **vtx=vlist->vtx -1 ; 
64
65
66     ts_uint nshell=(ts_uint)( sqrt((ts_double)(vlist->n-2)/5));
67 //    printf("nshell=%u\n",nshell);
68     ts_uint i,n0; // some for loop prereq
69     ts_int j,k;
70     ts_double dx,dy; // end loop prereq
71
72     /* topmost vertex */
8f6a69 73     vtx[1]->x=0.0;
SP 74     vtx[1]->y=0.0;
75     vtx[1]->z=z0*(ts_double)nshell;
7958e9 76     
SP 77     /* starting from to in circular order on pentagrams */    
78     for(i=1;i<=nshell;i++){
79         n0=2+5*i*(i-1)/2; //-1 would be for the reason that C index starts from 0 
8f6a69 80         vtx[n0]->x=0.0;
SP 81         vtx[n0]->y=(ts_double)i*xl0;
82         vtx[n0+i]->x=vtx[n0]->y*s1;
83         vtx[n0+i]->y=vtx[n0]->y*c1;
84         vtx[n0+2*i]->x=vtx[n0]->y*s2;
85         vtx[n0+2*i]->y=vtx[n0]->y*c2;
86         vtx[n0+3*i]->x=-vtx[n0+2*i]->x;
87         vtx[n0+3*i]->y=vtx[n0+2*i]->y;
88         vtx[n0+4*i]->x=-vtx[n0+i]->x;
89         vtx[n0+4*i]->y=vtx[n0+i]->y;
7958e9 90     }
SP 91
92     /* vertexes on the faces of the dipyramid */
93     for(i=1;i<=nshell;i++){
94         n0=2+5*i*(i-1)/2; // -1 would be because of C!
95         for(j=1;j<=i-1;j++){
8f6a69 96             dx=(vtx[n0]->x-vtx[n0+4*i]->x)/(ts_double)i;
SP 97             dy=(vtx[n0]->y-vtx[n0+4*i]->y)/(ts_double)i;
98             vtx[n0+4*i+j]->x=(ts_double)j*dx+vtx[n0+4*i]->x;
99             vtx[n0+4*i+j]->y=(ts_double)j*dy+vtx[n0+4*i]->y;
7958e9 100         }
SP 101         for(k=0;k<=3;k++){ // I would be worried about zero starting of for
8f6a69 102             dx=(vtx[n0+(k+1)*i]->x - vtx[n0+k*i]->x)/(ts_double) i;
SP 103             dy=(vtx[n0+(k+1)*i]->y - vtx[n0+k*i]->y)/(ts_double) i;
7958e9 104             for(j=1; j<=i-1;j++){
8f6a69 105                 vtx[n0+k*i+j]->x= (ts_double)j*dx+vtx[n0+k*i]->x;
SP 106                 vtx[n0+k*i+j]->y= (ts_double)j*dy+vtx[n0+k*i]->y;
7958e9 107             } 
SP 108         } 
109     }
110
111     for(i=1;i<=nshell;i++){
112         n0= 2+ 5*i*(i-1)/2;
113         for(j=0;j<=5*i-1;j++){
8f6a69 114         vtx[n0+j]->z= z0*(ts_double)(nshell-i);   // I would be worried about zero starting of for
7958e9 115         }
SP 116     }
117
118 /* for botom part of dipyramide we calculate the positions of vertices */
119     for(i=2+5*nshell*(nshell+1)/2;i<=vlist->n;i++){
8f6a69 120         vtx[i]->x=vtx[vlist->n - i +1]->x;
SP 121         vtx[i]->y=vtx[vlist->n - i +1]->y;
122         vtx[i]->z=-vtx[vlist->n - i +1]->z;
7958e9 123     }
SP 124
125     for(i=1;i<=vlist->n;i++){
126         for(j=1;j<=vlist->n;j++){
127             if(i!=j && vtx_distance_sq(vtx[i],vtx[j])<0.001){
128                 printf("Vertices %u and %u are the same!\n",i,j);
129             }
130         }
131     }
132     return TS_SUCCESS;
133 }
134
135
136
137 ts_bool init_vertex_neighbours(ts_vertex_list *vlist){
138     ts_vertex **vtx=vlist->vtx -1; // take a look at dipyramid function for comment.
139     const ts_double eps=0.001; //TODO: find out if you can use EPS from math.h
140     ts_uint i,j;
141     ts_double dist2; // Square of distance of neighbours
142     /*this is not required if we zero all data in vertex structure at initialization */
143     /*if we force zeroing at initialization this for loop can safely be deleted */
144     //for(i=1;i<=vlist->n;i++){
145     //    vtx[i].neigh_no=0;
146     //}
147     for(i=1;i<=vlist->n;i++){
148         for(j=1;j<=vlist->n;j++){
149             dist2=vtx_distance_sq(vtx[i],vtx[j]);
150             if( (dist2>eps) && (dist2<(A0*A0+eps))){ 
151     //if it is close enough, but not too much close (solves problem of comparing when i==j)
152                 vtx_add_neighbour(vtx[i],vtx[j]);
153             }
154         }
155     //        printf ("vertex %u ima %u sosedov!\n",i,vtx[i]->data->neigh_no);
156     }
157
158     return TS_SUCCESS;
159 }
160
b01cc1 161 // TODO: with new datastructure can be rewritten. Partially it is done, but it is complicated.
SP 162 ts_vertex_list *init_sort_neighbours(ts_bond_list *blist,ts_vertex_list *vlist){
7958e9 163     ts_vertex **vtx=vlist->vtx -1; // take a look at dipyramid function for comment.
SP 164     ts_uint i,l,j,jj,jjj,k=0;   
165     ts_double eps=0.001; // Take a look if EPS from math.h can be used
166
167 /*lets initialize memory for temporary vertex_list. Should we write a function instead */
b01cc1 168     ts_vertex_list *tvlist=vertex_list_copy(vlist);
7958e9 169     ts_vertex **tvtx=tvlist->vtx -1;  /* again to compensate for 0-indexing */
SP 170
171     ts_double dist2; // Square of distance of neighbours
172     ts_double direct; // Something, dont know what, but could be normal of some kind
173     for(i=1;i<=vlist->n;i++){
174         k++; // WHY i IS NOT GOOD??
8f6a69 175            vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->neigh[0]->idx+1]); //always add 1st
7958e9 176            jjj=1;
SP 177            jj=1;
8f6a69 178            for(l=2;l<=vtx[i]->neigh_no;l++){
SP 179                for(j=2;j<=vtx[i]->neigh_no;j++){
180                    dist2=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
181                    direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
182 // TODO: check if fabs can be used with all floating point types!!
7958e9 183                    if( (fabs(dist2-A0*A0)<=eps) && (direct>0.0) && (j!=jjj) ){
8f6a69 184                        vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->neigh[j-1]->idx+1]);
7958e9 185                        jjj=jj;
SP 186                        jj=j;
187                        break;
188                    }
189                }
190            }    
191     }
b01cc1 192 /* We use the temporary vertex for our main vertices and we abandon main
SP 193  * vertices, because their neighbours are not correctly ordered */
194    // tvtx=vlist->vtx;
195    // vlist->vtx=tvtx;
196    // tvlist->vtx=vtx;
197     vtx_list_free(vlist);
198 /* Let's make a check if the number of bonds is correct */
199     if((blist->n)!=3*(tvlist->n-2)){
200         ts_fprintf(stderr,"Number of bonds is %u should be %u!\n", blist->n, 3*(tvlist->n-2));
201         fatal("Number of bonds is not 3*(no_vertex-2).",4);
7958e9 202     }
SP 203
b01cc1 204     return tvlist;
7958e9 205 }
SP 206
207
208 ts_bool init_vesicle_bonds(ts_vesicle *vesicle){
209     ts_vertex_list *vlist=vesicle->vlist;
210     ts_bond_list *blist=vesicle->blist;
211     ts_vertex **vtx=vesicle->vlist->vtx - 1; // Because of 0 indexing
212 /* lets make correct clockwise ordering of in nearest neighbour list */
213     ts_uint i,j,k;
214     for(i=1;i<=vlist->n;i++){
215         for(j=i+1;j<=vlist->n;j++){
8f6a69 216             for(k=0;k<vtx[i]->neigh_no;k++){ // has changed 0 to < instead of 1 and <=
SP 217                 if(vtx[i]->neigh[k]==vtx[j]){  //if addresses matches it is the same
7958e9 218                     bond_add(blist,vtx[i],vtx[j]);
SP 219                     break;
220                 }
221             }
222         }
223     } 
224 /* Let's make a check if the number of bonds is correct */
225     if((blist->n)!=3*(vlist->n-2)){
226         ts_fprintf(stderr,"Number of bonds is %u should be %u!\n", blist->n, 3*(vlist->n-2));
227         fatal("Number of bonds is not 3*(no_vertex-2).",4);
228     }
229     return TS_SUCCESS;
230 }
231
232
233
234 ts_bool init_triangles(ts_vesicle *vesicle){
235     ts_uint i,j,jj,k;
236     ts_vertex **vtx=vesicle->vlist->vtx -1; // difference between 0 indexing and 1 indexing
237     ts_triangle_list *tlist=vesicle->tlist;
238     ts_double dist, direct;
239     ts_double eps=0.001; // can we use EPS from math.h?
240     k=0;
241     for(i=1;i<=vesicle->vlist->n;i++){
8f6a69 242         for(j=1;j<=vtx[i]->neigh_no;j++){
SP 243             for(jj=1;jj<=vtx[i]->neigh_no;jj++){
7958e9 244         //        ts_fprintf(stderr,"%u: (%u,%u) neigh_no=%u ",i,j,jj,vtx[i].neigh_no);
SP 245         //      ts_fprintf(stderr,"%e, %e",vtx[i].neigh[j-1]->x,vtx[i].neigh[jj-1]->x);
8f6a69 246                 dist=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
SP 247                 direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);                
248 // TODO: same as above                
249                 if(fabs(dist-A0*A0)<=eps && direct < 0.0 && vtx[i]->neigh[j-1]->idx+1 > i && vtx[i]->neigh[jj-1]->idx+1 >i){
250                     triangle_add(tlist,vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
7958e9 251                 }    
SP 252             }    
253         }
254     }
255 /* We check if all triangles have 3 vertices and if the number of triangles
256  * matches the theoretical value.
257  */
258     for(i=0;i<tlist->n;i++){
259         k=0;
260         for(j=0;j<3;j++){
41a035 261             if(tlist->tria[i]->vertex[j]!=NULL)
7958e9 262             k++;
SP 263         }
264             if(k!=3){
8f6a69 265                 fatal("Some triangles have less than 3 vertices..",4);
7958e9 266             }   
SP 267     } 
268     if(tlist->n!=2*(vesicle->vlist->n -2)){
269         ts_fprintf(stderr,"The number of triangles is %u but should be %u!\n",tlist->n,2*(vesicle->vlist->n -2));
270         fatal("The number of triangles doesn't match 2*(no_vertex -2).",4);
271     }
272     return TS_SUCCESS;
273 }
274
275
276
277 ts_bool init_triangle_neighbours(ts_vesicle *vesicle){
278     ts_uint i,j,nobo;
279     ts_vertex *i1,*i2,*i3,*j1,*j2,*j3;
280 //    ts_vertex **vtx=vesicle->vlist->vtx -1; // difference between 0 indexing and 1 indexing
281     ts_triangle_list *tlist=vesicle->tlist;
282     ts_triangle **tria=tlist->tria -1;
283     nobo=0;
284     for(i=1;i<=tlist->n;i++){
41a035 285         i1=tria[i]->vertex[0]; 
SP 286         i2=tria[i]->vertex[1]; 
287         i3=tria[i]->vertex[2]; 
7958e9 288         for(j=1;j<=tlist->n;j++){
SP 289             if(j==i) continue;
41a035 290             j1=tria[j]->vertex[0]; 
SP 291             j2=tria[j]->vertex[1]; 
292             j3=tria[j]->vertex[2]; 
7958e9 293             if((i1==j1 && i3==j2) || (i1==j2 && i3==j3) || (i1==j3 && i3==j1)){
SP 294                     triangle_add_neighbour(tria[i],tria[j]);
295                     nobo++;
296             }
297         }
298     }
299     for(i=1;i<=tlist->n;i++){
41a035 300         i1=tria[i]->vertex[0]; 
SP 301         i2=tria[i]->vertex[1]; 
302         i3=tria[i]->vertex[2]; 
7958e9 303         for(j=1;j<=tlist->n;j++){
SP 304             if(j==i) continue;
41a035 305             j1=tria[j]->vertex[0]; 
SP 306             j2=tria[j]->vertex[1]; 
307             j3=tria[j]->vertex[2]; 
7958e9 308             if((i1==j1 && i2==j3) || (i1==j3 && i2==j2) || (i1==j2 && i2==j1)){
SP 309                 triangle_add_neighbour(tria[i],tria[j]);
310                 nobo++;
311             }
312         }
313     }
314     for(i=1;i<=tlist->n;i++){
41a035 315         i1=tria[i]->vertex[0]; 
SP 316         i2=tria[i]->vertex[1]; 
317         i3=tria[i]->vertex[2]; 
7958e9 318         for(j=1;j<=tlist->n;j++){
SP 319             if(j==i) continue;
41a035 320             j1=tria[j]->vertex[0]; 
SP 321             j2=tria[j]->vertex[1]; 
322             j3=tria[j]->vertex[2]; 
7958e9 323             if((i2==j1 && i3==j3) || (i2==j3 && i3==j2) || (i2==j2 && i3==j1)){
SP 324                 triangle_add_neighbour(tria[i],tria[j]);
325                 nobo++;
326             }
327         }
328     }
329     if(nobo != vesicle->blist->n*2) {
330             ts_fprintf(stderr,"Number of triangles= %u, number of bonds= %u\n",nobo/2, vesicle->blist->n);
331             fatal("Number of triangle neighbour pairs differs from double the number of bonds!",4);
332     }
333     return TS_SUCCESS;
334 }
335
336
337 ts_bool init_common_vertex_triangle_neighbours(ts_vesicle *vesicle){
338     ts_uint i,j,jp,k;
339     ts_vertex *k1,*k2,*k3,*k4,*k5;
340     ts_vertex **vtx=vesicle->vlist->vtx -1; // difference between 0 indexing and 1 indexing
341     ts_triangle_list *tlist=vesicle->tlist;
342     ts_triangle **tria=tlist->tria -1;
343
344     for(i=1;i<=vesicle->vlist->n;i++){
8f6a69 345         for(j=1;j<=vtx[i]->neigh_no;j++){
SP 346             k1=vtx[i]->neigh[j-1];
7958e9 347             jp=j+1;
8f6a69 348             if(j == vtx[i]->neigh_no) jp=1;
SP 349             k2=vtx[i]->neigh[jp-1];
7958e9 350             for(k=1;k<=tlist->n;k++){        // VERY NON-OPTIMAL!!! too many loops (vlist.n * vtx.neigh * tlist.n )!
41a035 351                 k3=tria[k]->vertex[0];
SP 352                 k4=tria[k]->vertex[1];
353                 k5=tria[k]->vertex[2];
7958e9 354 //                ts_fprintf(stderr,"%u %u: k=(%u %u %u)\n",k1,k2,k3,k4,k5);
SP 355                 if((vtx[i]==k3 && k1==k4 && k2==k5) ||
356                 (vtx[i]==k4 && k1==k5 && k2==k3) ||
357                 (vtx[i]==k5 && k1==k3 && k2==k4)){
b01cc1 358
SP 359 //TODO: probably something wrong with neighbour distribution.
360 //                if(vtx[i]==k3 || vtx[i]==k4 || vtx[i]==k5){
dac2e5 361     //                    if(i==6) ts_fprintf(stdout, "Vtx[%u] > Added to tristar!\n",i);
7958e9 362                     vertex_add_tristar(vtx[i],tria[k]);
SP 363                 }
364             }
365         }
366 /*        ts_fprintf(stderr,"TRISTAR for %u (%u):",i-1,vtx[i].tristar_no);
367         for(j=0;j<vtx[i].tristar_no;j++){
368             ts_fprintf(stderr," %u,",vtx[i].tristar[j]->idx);
369         }
370         ts_fprintf(stderr,"\n"); */
371     }
372     return TS_SUCCESS;
373 }
374
375
376 ts_bool init_normal_vectors(ts_triangle_list *tlist){
377     /* Normals point INSIDE vesicle */
378     ts_uint k;
379     ts_triangle **tria=tlist->tria -1; //for 0 indexing
380     for(k=1;k<=tlist->n;k++){
381         triangle_normal_vector(tria[k]);    
382     }
383     return TS_SUCCESS;
384 }