Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2013-11-24 c85b513fcee3891a2c421925eb46179f44b55577
commit | author | age
d7639a 1 #include<stdlib.h>
SP 2 #include<math.h>
3 #include<string.h>
4 #include "general.h"
5 #include "vertex.h"
a10dd5 6 #include "bond.h"
d7639a 7 #include<stdio.h>
SP 8
73f967 9 ts_vertex_list *init_vertex_list(ts_uint N){    
d7639a 10     ts_int i;
a10dd5 11     ts_vertex_list *vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
73f967 12     
d7639a 13     if(N==0){
SP 14         err("Initialized vertex list with zero elements. Pointer set to NULL");
73f967 15         vlist->n=0;
SP 16         vlist->vtx=NULL;
17         return vlist;
d7639a 18     }
73f967 19     
8f6a69 20     vlist->vtx=(ts_vertex **)calloc(N,sizeof(ts_vertex *));
SP 21     if(vlist->vtx==NULL)
73f967 22         fatal("Fatal error reserving memory space for vertex list! Could number of requsted vertices be too large?", 100);
a10dd5 23     for(i=0;i<N;i++) {
8f6a69 24         vlist->vtx[i]=(ts_vertex *)calloc(1,sizeof(ts_vertex));
a10dd5 25         vlist->vtx[i]->idx=i;
c85b51 26     vlist->vtx[i]->locked=(ts_uint *)calloc(1,sizeof(ts_vertex));
074a17 27     /* initialize Ylm for spherical hamonics DONE in sh.c */
SP 28 /*    for(i=0;i<l;i++){
29         vlist->vtx[i]->Ylm[i]=(ts_double **)calloc(2*i+1,sizeof(ts_double *));
30         for(j=0;j<(2*i+1);j++){
31             clist->vtx[i]->Ylm[i][j]=(ts_double *)calloc(sizeof(ts_double));
32         }
33     }
34 */
35
36
a10dd5 37     }
73f967 38     vlist->n=N;
d7639a 39     return vlist;
SP 40 }
41
737714 42 ts_bool vtx_add_neighbour(ts_vertex *vtx, ts_vertex *nvtx){
d7639a 43     ts_uint i;
SP 44     /* no neighbour can be null! */
45     if(vtx==NULL || nvtx==NULL) return TS_FAIL;
46     
47     /*if it is already a neighbour don't add it to the list */
8f6a69 48     for(i=0; i<vtx->neigh_no;i++){
SP 49         if(vtx->neigh[i]==nvtx) return TS_FAIL;
d7639a 50     }
8f6a69 51     ts_uint nn=++vtx->neigh_no;
SP 52     vtx->neigh=(ts_vertex **)realloc(vtx->neigh, nn*sizeof(ts_vertex *));
53     vtx->neigh[nn-1]=nvtx;
3131dc 54 /* This was a bug in creating DIPYRAMID (the neighbours were not in right
SP 55  * order).
56  */
d7639a 57     /* pa se sosedu dodamo vertex */
SP 58     /*if it is already a neighbour don't add it to the list */
3131dc 59 /*
737714 60     for(i=0; i<nvtx->data->neigh_no;i++){
SP 61         if(nvtx->data->neigh[i]==vtx) return TS_FAIL;
d7639a 62     } 
a10dd5 63     nn=++nvtx->data->neigh_no;
737714 64     nvtx->data->neigh=(ts_vertex **)realloc(nvtx->data->neigh, nn*sizeof(ts_vertex *));
a10dd5 65     nvtx->data->neigh[nn-1]=vtx;
3131dc 66 */
d7639a 67
SP 68     return TS_SUCCESS;
69 }
a10dd5 70
9802f1 71 /* TODO: optimize this. test this. */
SP 72 ts_bool vtx_remove_neighbour(ts_vertex *vtx, ts_vertex *nvtx){
73 /* find a neighbour */
74 /* remove it from the list while shifting remaining neighbours up */
e19e79 75     ts_uint i,j=0;
8f6a69 76     for(i=0;i<vtx->neigh_no;i++){
e19e79 77 //        fprintf(stderr,"neigh_addr=%ld\n", (long)vtx->neigh[i]);
8f6a69 78         if(vtx->neigh[i]!=nvtx){
SP 79             vtx->neigh[j]=vtx->neigh[i];
9802f1 80             j++;
SP 81         }
82     }
e19e79 83 //    fprintf(stderr,"remove_neighbour: vtx1_addr=%ld, vtx2_addr=%ld\n",(long)vtx,(long)nvtx);
9802f1 84 /* resize memory. potentionally time consuming */
8f6a69 85     vtx->neigh_no--;
SP 86     vtx->neigh=(ts_vertex **)realloc(vtx->neigh,vtx->neigh_no*sizeof(ts_vertex *));
87     if(vtx->neigh == NULL && vtx->neigh_no!=0)
e19e79 88         fatal("(1) Reallocation of memory failed during removal of vertex neighbour in vtx_remove_neighbour",100);
SP 89 //fprintf(stderr,"first alloc");
9802f1 90 /* repeat for the neighbour */
SP 91 /* find a neighbour */
92 /* remove it from the list while shifting remaining neighbours up */
e19e79 93     j=0;
8f6a69 94     for(i=0;i<nvtx->neigh_no;i++){
SP 95         if(nvtx->neigh[i]!=vtx){
96             nvtx->neigh[j]=nvtx->neigh[i];
9802f1 97             j++;
SP 98         }
99     }
100 /* resize memory. potentionally time consuming. */
e19e79 101 //    fprintf(stderr,"Neigbours=%d\n",nvtx->neigh_no);
8f6a69 102     nvtx->neigh_no--;
e19e79 103             nvtx->neigh=(ts_vertex **)realloc(nvtx->neigh,nvtx->neigh_no*sizeof(ts_vertex *));
SP 104 //    fprintf(stderr,"Neigbours=%d\n",nvtx->neigh_no);
8f6a69 105     if(nvtx->neigh == NULL && nvtx->neigh_no!=0)
e19e79 106         fatal("(2) Reallocation of memory failed during removal of vertex neighbour in vtx_remove_neighbour",100);
9802f1 107
SP 108     return TS_SUCCESS;
109 }
7958e9 110
9802f1 111
SP 112
a10dd5 113 ts_bool vtx_add_bond(ts_bond_list *blist,ts_vertex *vtx1,ts_vertex *vtx2){
SP 114     ts_bond *bond;
115     bond=bond_add(blist,vtx1,vtx2);
116     if(bond==NULL) return TS_FAIL;
8f6a69 117     vtx1->bond_no++;
30ee9c 118     vtx2->bond_no++;
3131dc 119    // vtx2->data->bond_no++;
a10dd5 120
8f6a69 121     vtx1->bond=(ts_bond **)realloc(vtx1->bond, vtx1->bond_no*sizeof(ts_bond *)); 
30ee9c 122     vtx2->bond=(ts_bond **)realloc(vtx2->bond, vtx2->bond_no*sizeof(ts_bond *)); 
3131dc 123    // vtx2->data->bond=(ts_bond **)realloc(vtx2->data->bond, vtx2->data->bond_no*sizeof(ts_bond *)); 
8f6a69 124     vtx1->bond[vtx1->bond_no-1]=bond;
30ee9c 125     vtx2->bond[vtx2->bond_no-1]=bond;
8f6a69 126    // vtx2->ata->bond[vtx2->data->bond_no-1]=bond;
a10dd5 127     return TS_SUCCESS;
SP 128 }
129
130 ts_bool vtx_add_cneighbour(ts_bond_list *blist, ts_vertex *vtx1, ts_vertex *vtx2){
131     ts_bool retval;
132     retval=vtx_add_neighbour(vtx1,vtx2);
e19e79 133   //  retval=vtx_add_neighbour(vtx2,vtx1);
a10dd5 134     if(retval==TS_SUCCESS)
SP 135     retval=vtx_add_bond(blist,vtx1,vtx2); 
136     return retval;
137 }
138
9802f1 139 /*TODO: write and optimize this urgently before use! */
SP 140 ts_bool vtx_remove_cneighbour(ts_bond_list *blist, ts_vertex *vtx1, ts_vertex
141 *vtx2){
142 //    ts_bool retval;
143 /* remove the bond */
144 //retval=vtx_remove_bond(blist,vtx1,vtx2);
145 /* remove the vertices */
146     return TS_SUCCESS;
147 }
a10dd5 148
d7639a 149
737714 150 ts_bool vtx_free(ts_vertex  *vtx){
8f6a69 151     if(vtx->neigh!=NULL)   free(vtx->neigh);
SP 152     if(vtx->tristar!=NULL) free(vtx->tristar);
153     if(vtx->bond!=NULL)    free(vtx->bond);
737714 154     free(vtx);
d7639a 155     return TS_SUCCESS;
SP 156 }
157
73f967 158 ts_bool vtx_list_free(ts_vertex_list *vlist){
SP 159     int i;
160     for(i=0;i<vlist->n;i++){
8f6a69 161         if(vlist->vtx[i]!=NULL) vtx_free(vlist->vtx[i]);
73f967 162     }
8f6a69 163     //free(*(vlist->vtx));
737714 164     free(vlist->vtx);
73f967 165     free(vlist);
SP 166     return TS_SUCCESS;
167 }
d7639a 168
bb77ca 169 inline ts_double vtx_distance_sq(ts_vertex *vtx1, ts_vertex *vtx2){
d7639a 170     ts_double dist;
SP 171 #ifdef TS_DOUBLE_DOUBLE
8f6a69 172     dist=pow(vtx1->x-vtx2->x,2) + pow(vtx1->y-vtx2->y,2) + pow(vtx1->z-vtx2->z,2);
d7639a 173 #endif
SP 174 #ifdef TS_DOUBLE_LONGDOUBLE
8f6a69 175     dist=powl(vtx1->x-vtx2->x,2) + powl(vtx1->y-vtx2->y,2) + powl(vtx1->z-vtx2->z,2);
d7639a 176 #endif
SP 177 #ifdef TS_DOUBLE_FLOAT
8f6a69 178     dist=powf(vtx1->x-vtx2->x,2) + powf(vtx1->y-vtx2->y,2) + powf(vtx1->z-vtx2->z,2);
d7639a 179 #endif
SP 180     return(dist);
181 }
bb77ca 182
7958e9 183
SP 184
185 ts_bool vtx_set_global_values(ts_vesicle *vesicle){ 
186     ts_double xk=vesicle->bending_rigidity;
187     ts_uint i; 
188     for(i=0;i<vesicle->vlist->n;i++){
8f6a69 189         vesicle->vlist->vtx[i]->xk=xk;
7958e9 190     }
SP 191     return TS_SUCCESS;
192 }
193
194 inline ts_double vtx_direct(ts_vertex *vtx1, ts_vertex *vtx2, ts_vertex *vtx3){
8f6a69 195     ts_double dX2=vtx2->x-vtx1->x;
SP 196     ts_double dY2=vtx2->y-vtx1->y;
197     ts_double dZ2=vtx2->z-vtx1->z;
198     ts_double dX3=vtx3->x-vtx1->x;
199     ts_double dY3=vtx3->y-vtx1->y;
200     ts_double dZ3=vtx3->z-vtx1->z;
201     ts_double direct=vtx1->x*(dY2*dZ3 -dZ2*dY3)+ 
202         vtx1->y*(dZ2*dX3-dX2*dZ3)+
203         vtx1->z*(dX2*dY3-dY2*dX3);
7958e9 204     return(direct);    
SP 205 }
206
207
208 inline ts_bool vertex_add_tristar(ts_vertex *vtx, ts_triangle *tristarmem){
8f6a69 209     vtx->tristar_no++;
SP 210     vtx->tristar=(ts_triangle **)realloc(vtx->tristar,vtx->tristar_no*sizeof(ts_triangle *));
211     if(vtx->tristar==NULL){
7958e9 212             fatal("Reallocation of memory while adding tristar failed.",3);
SP 213     }
8f6a69 214     vtx->tristar[vtx->tristar_no-1]=tristarmem;
7958e9 215     return TS_SUCCESS;
SP 216 }
217
f74313 218
30ee9c 219 /* Insert neighbour is a function that is required in bondflip. It inserts a
SP 220  * neighbour exactly in the right place. */
221 inline ts_bool vtx_insert_neighbour(ts_vertex *vtx, ts_vertex *nvtx, ts_vertex *vtxm){
222 //nvtx is a vertex that is to be inserted after vtxm!
223         ts_uint i,j,midx;
224         vtx->neigh_no++;
225         if(vtxm==NULL ||  nvtx==NULL || vtx==NULL)
226                 fatal("vertex_insert_neighbour: one of pointers has been zero.. Cannot proceed.",3);
227         //We need to reallocate space! The pointer *neight must be zero if not having neighbours jey (if neigh_no was 0 at thime of calling
228         vtx->neigh=realloc(vtx->neigh,vtx->neigh_no*sizeof(ts_vertex *));
229         if(vtx->neigh == NULL){
230             fatal("Reallocation of memory failed during insertion of vertex neighbour in vertex_insert_neighbour",3);
231         }
232         midx=0;
233         for(i=0;i<vtx->neigh_no-1;i++) if(vtx->neigh[i]==vtxm) {midx=i; break;}
234      //   fprintf(stderr,"midx=%d, vseh=%d\n",midx,vtx->neigh_no-2);
235         if(midx==vtx->neigh_no-2) {
236             vtx->neigh[vtx->neigh_no-1]=nvtx;
237         } else {
238             for(j=vtx->neigh_no-2;j>midx;j--) {
239                 vtx->neigh[j+1]=vtx->neigh[j];
240 //                vtx->bond_length[j+1]=vtx->bond_length[j];
241 //                vtx->bond_length_dual[j+1]=vtx->bond_length_dual[j];
242             }
243             vtx->neigh[midx+1]=nvtx;
244         }
245     return TS_SUCCESS;
246 }
247
248
249 /* vtx remove tristar is required in  bondflip. */
250 /* TODO: Check whether it is important to keep the numbering of tristar
251  * elements in some order or not! */
252 inline ts_bool vtx_remove_tristar(ts_vertex *vtx, ts_triangle *tristar){
253     ts_uint i,j=0;
254     for(i=0;i<vtx->tristar_no;i++){
255         if(vtx->tristar[i]!=tristar){
256             vtx->tristar[j]=vtx->tristar[i];
257             j++;
258         }
259     }
260     vtx->tristar_no--;
261     vtx->tristar=realloc(vtx->tristar,vtx->tristar_no*sizeof(ts_triangle *));
262     if(vtx->neigh == NULL){
263             fatal("Reallocation of memory failed during insertion of vertex neighbour in vertex_add_neighbour",3);
264         }
265     return TS_SUCCESS;
266 }
267
268
269
f74313 270 /* ****************************************************************** */
SP 271 /* ***** New vertex copy operations. Inherently they are slow.  ***** */
272 /* ****************************************************************** */
273
b01cc1 274 ts_bool vtx_copy(ts_vertex *cvtx, ts_vertex *ovtx){
f74313 275     memcpy((void *)cvtx,(void *)ovtx,sizeof(ts_vertex));
8f6a69 276     cvtx->neigh=NULL;
SP 277     cvtx->neigh_no=0;
278     cvtx->tristar_no=0;
279     cvtx->bond_no=0;
280     cvtx->tristar=NULL;
281     cvtx->bond=NULL;
282     cvtx->cell=NULL;
b01cc1 283     return TS_SUCCESS;
f74313 284 }
dac2e5 285
SP 286 ts_bool vtx_duplicate(ts_vertex *cvtx, ts_vertex *ovtx){
287     memcpy((void *)cvtx,(void *)ovtx,sizeof(ts_vertex));
288     return TS_SUCCESS;
289 }
290
f74313 291 //TODO: needs to be done
SP 292 ts_vertex **vtx_neigh_copy(ts_vertex_list *vlist,ts_vertex *ovtx){
293         return NULL;
294 }
295
dac2e5 296
SP 297
f74313 298 ts_vertex_list *vertex_list_copy(ts_vertex_list *ovlist){
SP 299     ts_uint i;
300     ts_vertex_list *vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
b01cc1 301     vlist=memcpy((void *)vlist, (void *)ovlist, sizeof(ts_vertex_list));
f74313 302     ts_vertex **vtx=(ts_vertex **)malloc(vlist->n*sizeof(ts_vertex *));
b01cc1 303     vlist->vtx=vtx;
8f6a69 304     if(vlist->vtx==NULL)
b01cc1 305         fatal("Fatal error reserving memory space for vertex list! Could number of requsted vertices be too large?", 100);
SP 306     for(i=0;i<vlist->n;i++) {
8f6a69 307         vlist->vtx[i]=(ts_vertex *)calloc(1,sizeof(ts_vertex));
b01cc1 308         vlist->vtx[i]->idx=i;
SP 309         vtx_copy(vlist->vtx[i],ovlist->vtx[i]);
f74313 310     }
b01cc1 311
f74313 312     return vlist;
SP 313 }
200815 314
SP 315
316
317 ts_bool vertex_taint(ts_vertex *vtx, ts_uint level){
4e56fe 318     ts_uint i;
c85b51 319     (*(vtx->locked))++; //lock current vertex
4e56fe 320     if(level==0){ //if we reach last in recursion exit
200815 321         return TS_SUCCESS;    
SP 322     }    
4e56fe 323     for(i=0; i<vtx->neigh_no; i++){ //else recursive call self with decreased level
200815 324         vertex_taint(vtx->neigh[i], level-1);
SP 325     }
326     return TS_SUCCESS;
327 }
328
329 ts_bool vertex_untaint(ts_vertex *vtx, ts_uint level){
c85b51 330         (*(vtx->locked))--;
4e56fe 331     if(level==0){
200815 332         return TS_SUCCESS;    
SP 333     }    
334     ts_uint i;
335     for(i=0; i<vtx->neigh_no; i++){
336         vertex_untaint(vtx->neigh[i], level-1);
337     }
338     return TS_SUCCESS;
339 }
340
341 inline ts_bool vertex_tainted(ts_vertex *vtx, ts_uint level, ts_uint amount){
c85b51 342     if(*(vtx->locked)>amount) return 1;
c78004 343     if(level==0){
SP 344         return 0;
345     }
346     ts_uint i;
347     for(i=0; i<vtx->neigh_no;i++){
348         if(vertex_tainted(vtx->neigh[i], level-1, amount)) return 1;
349     }
4e56fe 350     
SP 351     return 0;
200815 352 }