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