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