| commit | author | age | ||
| aec47d | 1 | #include<stdlib.h> |
| SP | 2 | #include<math.h> |
| 3 | #include "general.h" | |
| 4 | #include "vertex.h" | |
| 5 | #include "bond.h" | |
| 6 | #include "triangle.h" | |
| 7 | #include "vesicle.h" | |
| 8 | #include "energy.h" | |
| 9 | #include "timestep.h" | |
| 10 | #include "cell.h" | |
| 30ee9c | 11 | #include "bondflip.h" |
| aec47d | 12 | //#include "io.h" |
| SP | 13 | #include<stdio.h> |
| 14 | ||
| 15 | ts_bool single_bondflip_timestep(ts_vesicle *vesicle, ts_bond *bond, ts_double *rn){ | |
| 16 | /*c Vertex and triangle (lm and lp) indexing for bond flip: | |
| 17 | c +----- k-------+ +----- k ------+ | |
| 18 | c |lm1 / | \ lp1 | |lm1 / \ lp1 | | |
| 19 | c | / | \ | | / \ | | |
| 20 | c |/ | \ | FLIP |/ lm \ | | |
| 21 | c km lm | lp kp ---> km ---------- kp | |
| 22 | c |\ | / | |\ lp / | | |
| 23 | c | \ | / | | \ / | | |
| 24 | c |lm2 \ | / lp2 | |lm2 \ / lp2 | | |
| 25 | c +------it------+ +----- it -----+ | |
| 26 | c | |
| 27 | */ | |
| 28 | ts_vertex *it=bond->vtx1; | |
| 29 | ts_vertex *k=bond->vtx2; | |
| 30 | ts_uint nei,neip,neim; | |
| b866cf | 31 | ts_uint i,j; |
| aec47d | 32 | ts_double oldenergy, delta_energy; |
| b866cf | 33 | ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL; |
| aec47d | 34 | |
| SP | 35 | ts_vertex *kp,*km; |
| 36 | ||
| 37 | if(it->neigh_no< 3) return TS_FAIL; | |
| 38 | if(k->neigh_no< 3) return TS_FAIL; | |
| 39 | if(k==NULL || it==NULL){ | |
| 40 | fatal("In bondflip, number of neighbours of k or it is less than 3!",999); | |
| 41 | } | |
| 42 | ||
| 30ee9c | 43 | nei=0; |
| aec47d | 44 | for(i=0;i<it->neigh_no;i++){ // Finds the nn of it, that is k |
| SP | 45 | if(it->neigh[i]==k){ |
| 46 | nei=i; | |
| 47 | break; | |
| 48 | } | |
| 49 | } | |
| 50 | neip=nei+1; // I don't like it.. Smells like I must have it in correct order | |
| 51 | neim=nei-1; | |
| 52 | if(neip>=it->neigh_no) neip=0; | |
| 53 | if((ts_int)neim<0) neim=it->neigh_no-1; /* casting is essential... If not | |
| 54 | there the neim is never <0 !!! */ | |
| 55 | // fprintf(stderr,"The numbers are: %u %u\n",neip, neim); | |
| 56 | km=it->neigh[neim]; // We located km and kp | |
| 57 | kp=it->neigh[neip]; | |
| 58 | ||
| 59 | if(km==NULL || kp==NULL){ | |
| 60 | fatal("In bondflip, cannot determine km and kp!",999); | |
| 61 | } | |
| 62 | ||
| 63 | // fprintf(stderr,"I WAS HERE! after the 4 vertices are known!\n"); | |
| 64 | ||
| 65 | /* test if the membrane is wrapped too much, so that kp is nearest neighbour of | |
| 66 | * km. If it is true, then don't flip! */ | |
| 67 | for(i=0;i<km->neigh_no;i++){ | |
| 68 | if(km->neigh[i] == kp) return TS_FAIL; | |
| 69 | } | |
| 70 | // fprintf(stderr,"Membrane didn't wrap too much.. Continue.\n"); | |
| 71 | /* if bond would be too long, return... */ | |
| 30ee9c | 72 | if(vtx_distance_sq(km,kp) > vesicle->dmax ) return TS_FAIL; |
| aec47d | 73 | // fprintf(stderr,"Bond will not be too long.. Continue.\n"); |
| SP | 74 | |
| 75 | /* we make a bond flip. this is different than in original fortran */ | |
| b866cf | 76 | // find lm, lp |
| aec47d | 77 | // 1. step. We find lm and lp from k->tristar ! |
| SP | 78 | for(i=0;i<it->tristar_no;i++){ |
| 79 | for(j=0;j<k->tristar_no;j++){ | |
| 80 | if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik | |
| 81 | if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1] | |
| 82 | == km || it->tristar[i]->vertex[2]== km )){ | |
| 83 | lm=it->tristar[i]; | |
| 84 | // lmidx=i; | |
| 85 | } | |
| 86 | else | |
| 87 | { | |
| 88 | lp=it->tristar[i]; | |
| 89 | // lpidx=i; | |
| 90 | } | |
| 91 | ||
| 92 | } | |
| 93 | } | |
| 94 | } | |
| 95 | if(lm==NULL || lp==NULL) fatal("ts_flip_bond: Cannot find triangles lm and lp!",999); | |
| 96 | ||
| 97 | //we look for important triangles lp1 and lm2. | |
| 98 | ||
| 99 | for(i=0;i<k->tristar_no;i++){ | |
| 100 | for(j=0;j<kp->tristar_no;j++){ | |
| 101 | if((k->tristar[i] == kp->tristar[j]) && k->tristar[i]!=lp){ //ce gre za skupen trikotnik | |
| 102 | lp1=k->tristar[i]; | |
| 103 | } | |
| 104 | } | |
| 105 | } | |
| 106 | ||
| 107 | for(i=0;i<it->tristar_no;i++){ | |
| 108 | for(j=0;j<km->tristar_no;j++){ | |
| 109 | if((it->tristar[i] == km->tristar[j]) && it->tristar[i]!=lm){ //ce gre za skupen trikotnik | |
| 110 | lm2=it->tristar[i]; | |
| 111 | } | |
| 112 | } | |
| 113 | } | |
| 114 | ||
| 115 | if(lm2==NULL || lp1==NULL) fatal("ts_flip_bond: Cannot find triangles lm2 and lp1!",999); | |
| b866cf | 116 | // fprintf(stderr,"I WAS HERE! Before bondflip!\n"); |
| M | 117 | // fprintf(stderr,"%e, %e, %e\n", lm->xnorm, lm->ynorm, lm->znorm); |
| 118 | ||
| 119 | ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1); | |
| 120 | // fprintf(stderr,"I WAS HERE! Bondflip successful!\n"); | |
| aec47d | 121 | |
| b866cf | 122 | oldenergy=0; |
| M | 123 | oldenergy+=k->xk* k->energy; |
| 124 | oldenergy+=kp->xk* kp->energy; | |
| 125 | oldenergy+=km->xk* km->energy; | |
| 126 | oldenergy+=it->xk* it->energy; | |
| 127 | //Neigbours don't change its energy. | |
| aec47d | 128 | |
| SP | 129 | |
| 130 | ||
| b866cf | 131 | /* Calculating the new energy */ |
| M | 132 | delta_energy=0; |
| 133 | delta_energy+=k->xk* k->energy; | |
| 134 | delta_energy+=kp->xk* kp->energy; | |
| 135 | delta_energy+=km->xk* km->energy; | |
| 136 | delta_energy+=it->xk* it->energy; | |
| 137 | //Neigbours don't change its energy. | |
| 138 | delta_energy-=oldenergy; | |
| 139 | // fprintf(stderr,"I WAS HERE! Got energy!\n"); | |
| 140 | /* MONTE CARLO */ | |
| 141 | if(delta_energy>=0){ | |
| 142 | #ifdef TS_DOUBLE_DOUBLE | |
| 143 | if(exp(-delta_energy)< drand48() ) | |
| 144 | #endif | |
| 145 | #ifdef TS_DOUBLE_FLOAT | |
| 146 | if(expf(-delta_energy)< (ts_float)drand48()) | |
| 147 | #endif | |
| 148 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 149 | if(expl(-delta_energy)< (ts_ldouble)drand48()) | |
| 150 | #endif | |
| 151 | { | |
| 152 | //not accepted, reverting changes | |
| 153 | // fprintf(stderr,"Failed to move, due to MC\n"); | |
| aec47d | 154 | |
| b866cf | 155 | ts_flip_bond(kp,km,k,it, bond, lm,lp,lm2,lp1); |
| M | 156 | // fprintf(stderr,"%e, %e, %e\n", lp->xnorm, lp->ynorm, lp->znorm); |
| 157 | return TS_FAIL; | |
| 158 | } | |
| 159 | } | |
| 160 | /* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */ | |
| 161 | // fprintf(stderr,"SUCCESS!!!\n"); | |
| 162 | return TS_SUCCESS; | |
| aec47d | 163 | } |
| SP | 164 | |
| 165 | ||
| b866cf | 166 | ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp, |
| M | 167 | ts_bond *bond, ts_triangle *lm, ts_triangle *lp, ts_triangle *lm2, ts_triangle *lp1){ |
| 168 | ||
| 169 | ts_uint i; //lmidx, lpidx; | |
| 170 | if(k==NULL || it==NULL || km==NULL || kp==NULL){ | |
| 171 | fatal("ts_flip_bond: You called me with invalid pointers to vertices",999); | |
| 172 | } | |
| aec47d | 173 | // 2. step. We change the triangle vertices... (actual bond flip) |
| SP | 174 | for(i=0;i<3;i++) if(lm->vertex[i]== it) lm->vertex[i]= kp; |
| 175 | for(i=0;i<3;i++) if(lp->vertex[i]== k) lp->vertex[i]= km; | |
| e19e79 | 176 | //fprintf(stderr,"2. step: actual bondflip made\n"); |
| aec47d | 177 | // 2a. step. If any changes in triangle calculations must be done, do it here! |
| SP | 178 | // * normals are recalculated here |
| 179 | triangle_normal_vector(lp); | |
| 180 | triangle_normal_vector(lm); | |
| e19e79 | 181 | //fprintf(stderr,"2a. step: triangle normals recalculated\n"); |
| aec47d | 182 | // 3. step. Correct neighbours in vertex_list |
| SP | 183 | |
| 184 | ||
| 30ee9c | 185 | vtx_remove_neighbour(k,it); |
| 306559 | 186 | // vtx_remove_neighbour(it,k); |
| e19e79 | 187 | //fprintf(stderr,"3. step (PROGRESS): removed k and it neighbours\n"); |
| 306559 | 188 | |
| aec47d | 189 | //Tukaj pa nastopi tezava... Kam dodati soseda? |
| 30ee9c | 190 | vtx_insert_neighbour(km,kp,k); |
| SP | 191 | vtx_insert_neighbour(kp,km,it); |
| aec47d | 192 | // vertex_add_neighbour(km,kp); //pazi na vrstni red. |
| SP | 193 | // vertex_add_neighbour(kp,km); |
| e19e79 | 194 | //fprintf(stderr,"3. step: vertex neighbours corrected\n"); |
| aec47d | 195 | |
| SP | 196 | // 3a. step. If any changes to ts_vertex, do it here! |
| 197 | // bond_length calculatons not required for it is done in energy.c | |
| 198 | ||
| 199 | // 4. step. Correct bond_list (don't know why I still have it!) | |
| 200 | bond->vtx1=km; | |
| 201 | bond->vtx2=kp; | |
| 202 | //fprintf(stderr,"4. step: bondlist corrected\n"); | |
| 203 | ||
| 204 | ||
| 205 | // 5. step. Correct neighbouring triangles | |
| 206 | ||
| 207 | triangle_remove_neighbour(lp,lp1); | |
| e19e79 | 208 | // fprintf(stderr,".\n"); |
| aec47d | 209 | triangle_remove_neighbour(lp1,lp); |
| SP | 210 | // fprintf(stderr,".\n"); |
| 211 | triangle_remove_neighbour(lm,lm2); | |
| 212 | // fprintf(stderr,".\n"); | |
| 213 | triangle_remove_neighbour(lm2,lm); | |
| 214 | ||
| 215 | triangle_add_neighbour(lm,lp1); | |
| 216 | triangle_add_neighbour(lp1,lm); | |
| 217 | triangle_add_neighbour(lp,lm2); //Vrstni red?! | |
| 218 | triangle_add_neighbour(lm2,lp); | |
| 219 | ||
| 220 | //fprintf(stderr,"5. step: triangle neigbours corrected\n"); | |
| 221 | ||
| 222 | ||
| 223 | // 6. step. Correct tristar for vertices km, kp, k and it | |
| 224 | vertex_add_tristar(km,lp); // Preveri vrstni red! | |
| 225 | vertex_add_tristar(kp,lm); | |
| 30ee9c | 226 | vtx_remove_tristar(it,lm); |
| SP | 227 | vtx_remove_tristar(k,lp); |
| aec47d | 228 | //fprintf(stderr,"6. step: tristar corrected\n"); |
| SP | 229 | energy_vertex(k); |
| 230 | energy_vertex(kp); | |
| 231 | energy_vertex(km); | |
| 232 | energy_vertex(it); | |
| 233 | // END modifications to data structure! | |
| 234 | return TS_SUCCESS; | |
| 235 | } | |
