Trisurf Monte Carlo simulator
Samo Penic
2012-02-23 e016c401c98b6a33392b551b2c40fccded21efe3
commit | author | age
d7639a 1 #include<stdlib.h>
SP 2 #include "general.h"
3 #include "vertex.h"
bb77ca 4
SP 5 ts_cell_list  *init_cell_list(ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
d7639a 6     ts_uint i;
bb77ca 7     ts_uint nocells=ncmax1*ncmax2*ncmax3;
SP 8     ts_cell_list *clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
9     if(clist==NULL) fatal("Error while allocating memory for cell list!",100);
10
11     clist->ncmax[0]=ncmax1;
12     clist->ncmax[1]=ncmax2;
13     clist->ncmax[2]=ncmax3;
14     clist->cellno=nocells;
d7639a 15     clist->dcell=1.0/(1.0 + stepsize);
SP 16     clist->shift=(ts_double) clist->ncmax[0]/2;
bb77ca 17
SP 18     clist->cell=(ts_cell **)malloc(nocells*sizeof(ts_cell *));
19     if(clist->cell==NULL) fatal("Error while allocating memory for cell list! ncmax too large?",101);
20
d7639a 21     for(i=0;i<nocells;i++){
bb77ca 22         clist->cell[i]=(ts_cell *)malloc(sizeof(ts_cell));
SP 23         if(clist->cell[i]==NULL) fatal("Error while allocating memory for cell list! ncmax too large?",102);
24         clist->cell[i]->idx=i+1; // We enumerate cells! Probably never required!
25         clist->cell[i]->data=(ts_cell_data *)calloc(1,sizeof(ts_cell_data));
d7639a 26     }
bb77ca 27     return clist;
SP 28 }
29
30 ts_bool cell_free(ts_cell* cell){
31     if(cell->data!=NULL){
32         if(cell->data->vertex!=NULL) free(cell->data->vertex);
33         free(cell->data);
34     }
35     free(cell);
d7639a 36     return TS_SUCCESS;
SP 37 }
38
39 ts_bool cell_list_free(ts_cell_list *clist){
40     ts_uint i;
bb77ca 41     if(clist==NULL) return TS_FAIL;
SP 42     ts_uint nocells=clist->cellno;
d7639a 43     for(i=0;i<nocells;i++)
bb77ca 44          if(clist->cell[i] != NULL) cell_free(clist->cell[i]);
d7639a 45     free(clist->cell);
bb77ca 46     free(clist);
d7639a 47     return TS_SUCCESS;
SP 48 }
49
50
bb77ca 51 //TODO: not debugged at all!
SP 52 inline ts_uint vertex_self_avoidance(ts_vesicle *vesicle, ts_vertex *vtx){
53     ts_uint cellidx;
54     ts_uint ncx, ncy,ncz;
55     ts_cell_list *clist=vesicle->clist;
56     ncx=(ts_uint)((vtx->data->x-vesicle->cm[0])*clist->dcell+clist->shift);
57     ncy=(ts_uint)((vtx->data->y-vesicle->cm[1])*clist->dcell+clist->shift);
58     ncz=(ts_uint)((vtx->data->z-vesicle->cm[2])*clist->dcell+clist->shift);
59
60     if(ncx == clist->ncmax[0]-1 || ncx == 2){
d7639a 61         fatal("Vesicle is positioned outside the cell covered area. Coordinate x is the problem.",1500);
SP 62     }
bb77ca 63     if(ncy == clist->ncmax[1]-1 || ncy == 2){
d7639a 64         fatal("Vesicle is positioned outside the cell covered area. Coordinate y is the problem.",1500);
SP 65     }
bb77ca 66     if(ncz == clist->ncmax[2]-1 || ncz == 2){
d7639a 67         fatal("Vesicle is positioned outside the cell covered area. Coordinate z is the problem.",1500);
SP 68     }
bb77ca 69     cellidx=ncz+(ncy-1)*clist->ncmax[2] + (ncx-1)*clist->ncmax[2]* 
SP 70                                     clist->ncmax[1] - 1; // -1 is because of 0 based indexing
71         return cellidx;
d7639a 72 }
SP 73
74
bb77ca 75 //TODO: looks ok, but debug anyway in the future
d7639a 76 ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){
bb77ca 77     cell->data->nvertex++;
SP 78     cell->data->vertex=(ts_vertex **)realloc(cell->data->vertex,cell->data->nvertex*sizeof(ts_vertex *));
dac2e5 79         if(cell->data->vertex == NULL){
b01cc1 80             fatal("Reallocation of memory failed during insertion of vertex in cell_add_vertex",3);
d7639a 81         }
bb77ca 82     cell->data->vertex[cell->data->nvertex-1]=vtx;
d7639a 83     return TS_SUCCESS;
SP 84 }
85
bb77ca 86
SP 87 ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist){
d7639a 88     ts_uint i;
SP 89     for(i=0;i<clist->cellno;i++){
bb77ca 90         if(clist->cell[i]->data->vertex != NULL){
SP 91             free(clist->cell[i]->data->vertex);
92             clist->cell[i]->data->vertex=NULL;
d7639a 93         }
bb77ca 94         clist->cell[i]->data->nvertex=0;
d7639a 95     }
SP 96     return TS_SUCCESS;
97 }
98
bb77ca 99 // TODO: compiles ok, but it is completely untested and undebugged. It was debugged before rewrite, but this was long time ago.
d7639a 100 ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx, ts_vertex *tvtx){
SP 101     ts_uint ncx,ncy,ncz,remainder,cell_occupation;
bb77ca 102     ts_uint i,j,k,l,neigh_cidx;
d7639a 103     ts_double dist;
bb77ca 104     ncx=(cellidx+1)/(clist->ncmax[2]*clist->ncmax[1])+1; //+1 because of zero indexing.
d7639a 105     remainder=(cellidx+1)%(clist->ncmax[2]*clist->ncmax[1]);
SP 106     ncy=remainder/clist->ncmax[2]+1;
107     ncz=remainder%clist->ncmax[2];
108 //    fprintf(stderr,"here are ncx,ncy,ncz=%i,%i,%i\n",ncx,ncy,ncz);
109
110     for(i=ncx-1;i<=ncx+1;i++){
111         for(j=ncy-1;j<=ncy+1;j++){
112             for(k=ncz-1;k<=ncz+1;k++){
113                 neigh_cidx=k+(j-1)*clist->ncmax[2]+(i-1)*clist->ncmax[2]*clist->ncmax[1] -1;
114           //      fprintf(stderr,"neigh_cell_index=%i\n",neigh_cidx);
bb77ca 115                 cell_occupation=clist->cell[neigh_cidx]->data->nvertex;
d7639a 116           //      fprintf(stderr, "cell_occupation=%i\n",cell_occupation);
SP 117                 if(cell_occupation>clist->max_occupancy){
118                     fatal("Neighbouring cell occupation more than set max_occupancy value.",2000);
119                 }
120 // Now we check whether we didn't come close to some other vertices in the same
121 // cell!
122                 if(cell_occupation>1){
123                     for(l=0;l<cell_occupation;l++){
bb77ca 124                         if(clist->cell[neigh_cidx]->data->vertex[l]!=vtx){
d7639a 125                     //        fprintf(stderr,"calling dist on vertex %i\n",l);
bb77ca 126                            dist=vtx_distance_sq(clist->cell[neigh_cidx]->data->vertex[l],tvtx);
d7639a 127                     //        fprintf(stderr,"dist was %f\n",dist);
SP 128                             if(dist<1) return TS_FAIL;
129                         }
130                     }
131                 }
132             }
133         }
134     }
135     
136     return TS_SUCCESS;
137 }