Trisurf Monte Carlo simulator
Samo Penic
2019-10-19 dba66b4f8a25cb7eb18de5f7ffca97e179d7a5d9
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
d7639a 2 #include<stdlib.h>
SP 3 #include "general.h"
4 #include "vertex.h"
bb77ca 5
SP 6 ts_cell_list  *init_cell_list(ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
d7639a 7     ts_uint i;
bb77ca 8     ts_uint nocells=ncmax1*ncmax2*ncmax3;
SP 9     ts_cell_list *clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
10     if(clist==NULL) fatal("Error while allocating memory for cell list!",100);
11
12     clist->ncmax[0]=ncmax1;
13     clist->ncmax[1]=ncmax2;
14     clist->ncmax[2]=ncmax3;
15     clist->cellno=nocells;
d7639a 16     clist->dcell=1.0/(1.0 + stepsize);
SP 17     clist->shift=(ts_double) clist->ncmax[0]/2;
bb77ca 18
SP 19     clist->cell=(ts_cell **)malloc(nocells*sizeof(ts_cell *));
20     if(clist->cell==NULL) fatal("Error while allocating memory for cell list! ncmax too large?",101);
21
d7639a 22     for(i=0;i<nocells;i++){
34d3de 23         clist->cell[i]=(ts_cell *)calloc(1,sizeof(ts_cell));
bb77ca 24         if(clist->cell[i]==NULL) fatal("Error while allocating memory for cell list! ncmax too large?",102);
SP 25         clist->cell[i]->idx=i+1; // We enumerate cells! Probably never required!
d7639a 26     }
bb77ca 27     return clist;
SP 28 }
29
30 ts_bool cell_free(ts_cell* cell){
34d3de 31     if(cell->vertex!=NULL) free(cell->vertex);
bb77ca 32     free(cell);
d7639a 33     return TS_SUCCESS;
SP 34 }
35
36 ts_bool cell_list_free(ts_cell_list *clist){
37     ts_uint i;
bb77ca 38     if(clist==NULL) return TS_FAIL;
SP 39     ts_uint nocells=clist->cellno;
d7639a 40     for(i=0;i<nocells;i++)
bb77ca 41          if(clist->cell[i] != NULL) cell_free(clist->cell[i]);
d7639a 42     free(clist->cell);
bb77ca 43     free(clist);
d7639a 44     return TS_SUCCESS;
SP 45 }
46
47
bb77ca 48 inline ts_uint vertex_self_avoidance(ts_vesicle *vesicle, ts_vertex *vtx){
SP 49     ts_uint cellidx;
50     ts_uint ncx, ncy,ncz;
51     ts_cell_list *clist=vesicle->clist;
8f6a69 52     ncx=(ts_uint)((vtx->x-vesicle->cm[0])*clist->dcell+clist->shift);
SP 53     ncy=(ts_uint)((vtx->y-vesicle->cm[1])*clist->dcell+clist->shift);
54     ncz=(ts_uint)((vtx->z-vesicle->cm[2])*clist->dcell+clist->shift);
bb77ca 55
a63f17 56     if(ncx >= clist->ncmax[0]-1 || ncx <= 2){
d7639a 57         fatal("Vesicle is positioned outside the cell covered area. Coordinate x is the problem.",1500);
SP 58     }
a63f17 59     if(ncy >= clist->ncmax[1]-1 || ncy <= 2){
d7639a 60         fatal("Vesicle is positioned outside the cell covered area. Coordinate y is the problem.",1500);
SP 61     }
a63f17 62     if(ncz >= clist->ncmax[2]-1 || ncz <= 2){
d7639a 63         fatal("Vesicle is positioned outside the cell covered area. Coordinate z is the problem.",1500);
SP 64     }
bb77ca 65     cellidx=ncz+(ncy-1)*clist->ncmax[2] + (ncx-1)*clist->ncmax[2]* 
SP 66                                     clist->ncmax[1] - 1; // -1 is because of 0 based indexing
67         return cellidx;
d7639a 68 }
SP 69
70
a63f17 71 inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){
SP 72     ts_uint i;
73     for(i=0;i<cell->nvertex;i++){
74         if(cell->vertex[i]==vtx){
75     //vertex is already in the cell!
76             //fprintf(stderr,"VTX in the cell!\n");
77             return TS_FAIL;
78         }
79     }
80             //fprintf(stderr,"VTX added to the cell!\n");
34d3de 81     cell->nvertex++;
SP 82     cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *));
83         if(cell->vertex == NULL){
b01cc1 84             fatal("Reallocation of memory failed during insertion of vertex in cell_add_vertex",3);
d7639a 85         }
34d3de 86     cell->vertex[cell->nvertex-1]=vtx;
a63f17 87     vtx->cell=cell;
d7639a 88     return TS_SUCCESS;
SP 89 }
90
a63f17 91 inline ts_bool cell_remove_vertex(ts_cell *cell, ts_vertex *vtx){
SP 92    ts_uint i,j=0;
93     for(i=0;i<cell->nvertex;i++){
94         if(cell->vertex[i]!=vtx){
95             cell->vertex[j]=cell->vertex[i];
96             j++;
97         }
98     }
99     if(j==i){
100     fatal("Vertex was not in the cell!",3);
101     } 
102     //fprintf(stderr, "Vertex deleted from the cell!\n");
103
104 /* resize memory. potentionally time consuming */
105     cell->nvertex--;
106     cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *));
107     if(vtx->neigh == NULL && vtx->neigh_no!=0)
108         if(cell->vertex == NULL){
109             fatal("Reallocation of memory failed during removal of vertex in cell_remove_vertex",3);
110         }
111     return TS_SUCCESS;
112 }
bb77ca 113
SP 114 ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist){
d7639a 115     ts_uint i;
SP 116     for(i=0;i<clist->cellno;i++){
34d3de 117         if(clist->cell[i]->vertex != NULL){
SP 118             free(clist->cell[i]->vertex);
119             clist->cell[i]->vertex=NULL;
d7639a 120         }
34d3de 121         clist->cell[i]->nvertex=0;
d7639a 122     }
SP 123     return TS_SUCCESS;
124 }
125
ea1cce 126
a63f17 127 ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx){
d7639a 128     ts_uint ncx,ncy,ncz,remainder,cell_occupation;
bb77ca 129     ts_uint i,j,k,l,neigh_cidx;
d7639a 130     ts_double dist;
bb77ca 131     ncx=(cellidx+1)/(clist->ncmax[2]*clist->ncmax[1])+1; //+1 because of zero indexing.
d7639a 132     remainder=(cellidx+1)%(clist->ncmax[2]*clist->ncmax[1]);
SP 133     ncy=remainder/clist->ncmax[2]+1;
134     ncz=remainder%clist->ncmax[2];
135 //    fprintf(stderr,"here are ncx,ncy,ncz=%i,%i,%i\n",ncx,ncy,ncz);
136
137     for(i=ncx-1;i<=ncx+1;i++){
138         for(j=ncy-1;j<=ncy+1;j++){
139             for(k=ncz-1;k<=ncz+1;k++){
140                 neigh_cidx=k+(j-1)*clist->ncmax[2]+(i-1)*clist->ncmax[2]*clist->ncmax[1] -1;
141           //      fprintf(stderr,"neigh_cell_index=%i\n",neigh_cidx);
34d3de 142                 cell_occupation=clist->cell[neigh_cidx]->nvertex;
d7639a 143           //      fprintf(stderr, "cell_occupation=%i\n",cell_occupation);
SP 144                 if(cell_occupation>clist->max_occupancy){
42190f 145             ts_fprintf(stderr,"max occupancy= %d, cell occupation= %d", clist->max_occupancy, cell_occupation);
d7639a 146                     fatal("Neighbouring cell occupation more than set max_occupancy value.",2000);
SP 147                 }
148 // Now we check whether we didn't come close to some other vertices in the same
149 // cell!
a18ed8 150                 if(cell_occupation>0){
d7639a 151                     for(l=0;l<cell_occupation;l++){
a63f17 152
SP 153                 //carefull with this checks!
555fd8 154                         if(clist->cell[neigh_cidx]->vertex[l]!=vtx){
d7639a 155                     //        fprintf(stderr,"calling dist on vertex %i\n",l);
a63f17 156                            dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],vtx);
555fd8 157
M 158 //                if(vtx->idx==1)
159 //                fprintf(stderr,"VTX(0) ima bliznji vertex z indeksom, %d, tipa %d \n", clist->cell[neigh_cidx]->vertex[l]->idx, clist->cell[neigh_cidx]->vertex[l]->id);
160 //                if(vtx->idx==0 && clist->cell[neigh_cidx]->vertex[l]->idx==0)
161 //                            fprintf(stderr,"*** dist was %f\n",dist);
162                 
ea1cce 163                             if(dist<=1.0 || (dist<=clist->dmin_interspecies && (clist->cell[neigh_cidx]->vertex[l]->id != vtx->id))) return TS_FAIL;
d7639a 164                         }
SP 165                     }
166                 }
167             }
168         }
a63f17 169     } 
d7639a 170     return TS_SUCCESS;
SP 171 }