| commit | author | age | ||
| 7f6076 | 1 | /* vim: set ts=4 sts=4 sw=4 noet : */ |
| 460c2a | 2 | #include<stdlib.h> |
| SP | 3 | #include<stdio.h> |
| 4 | #include<string.h> | |
| 5 | #include<math.h> | |
| 6 | #include "general.h" | |
| 7 | #include "constvol.h" | |
| 8 | #include "triangle.h" | |
| 9 | #include "energy.h" | |
| 10 | #include "vertex.h" | |
| 11 | #include "cell.h" | |
| 12 | ||
| fbcbdf | 13 | ts_bool constvolume(ts_vesicle *vesicle, ts_vertex *vtx_avoid, ts_double Vol, ts_double *retEnergy, ts_vertex **vtx_moved_retval, ts_vertex **vtx_backup){ |
| SP | 14 | ts_vertex *vtx_moved; |
| 460c2a | 15 | ts_uint vtxind,i,j; |
| 41637a | 16 | ts_uint Ntries=100; |
| 460c2a | 17 | ts_vertex *backupvtx; |
| fe5069 | 18 | ts_double Rv, dh, dvol, volFirst, voldiff, oenergy,delta_energy; |
| 460c2a | 19 | backupvtx=(ts_vertex *)calloc(sizeof(ts_vertex),10); |
| fbcbdf | 20 | ts_double l0 = (1.0 + sqrt(vesicle->dmax))/2.0; //make this a global constant if necessary |
| 460c2a | 21 | for(i=0;i<Ntries;i++){ |
| SP | 22 | vtxind=rand() % vesicle->vlist->n; |
| 23 | vtx_moved=vesicle->vlist->vtx[vtxind]; | |
| 41637a | 24 | /* chosen vertex must not be a nearest neighbour. WASTODO: probably must. |
| SP | 25 | * DONE: solved below, when using different algorithm. |
| fe5069 | 26 | * extend search in case of bondflip */ |
| 41637a | 27 | /* if(vtx_moved==vtx_avoid) continue; |
| 460c2a | 28 | for(j=0;j<vtx_moved->neigh_no;j++){ |
| SP | 29 | if(vtx_moved->neigh[j]==vtx_avoid) continue; |
| 30 | } | |
| 41637a | 31 | */ |
| SP | 32 | /* different check of nearest neighbour (and second nearest neighbour) check. |
| 33 | * Checking the distance between vertex and vertex to be moved to assure | |
| 34 | * constant volume. Solves upper todo problem. */ | |
| 35 | ||
| 36 | if(vtx_distance_sq(vtx_moved,vtx_avoid)<16.0*vesicle->dmax){ | |
| 37 | continue; | |
| 38 | } | |
| 39 | ||
| 460c2a | 40 | memcpy((void *)&backupvtx[0],(void *)vtx_moved,sizeof(ts_vertex)); |
| SP | 41 | |
| fe5069 | 42 | //move vertex in specified direction. first try, test move! |
| 460c2a | 43 | Rv=sqrt(pow(vtx_moved->x,2)+pow(vtx_moved->y,2)+pow(vtx_moved->z,2)); |
| fbcbdf | 44 | dh=2.0*Vol/(sqrt(3.0)*l0*l0); |
| SP | 45 | vtx_moved->x=vtx_moved->x*(1.0-dh/Rv); |
| 46 | vtx_moved->y=vtx_moved->y*(1.0-dh/Rv); | |
| 47 | vtx_moved->z=vtx_moved->z*(1.0-dh/Rv); | |
| 460c2a | 48 | |
| 41637a | 49 | |
| SP | 50 | //SKIPPING FIRST CHECK of CONSTRAINTS. This is not a final move. |
| 460c2a | 51 | //check for constraints |
| 41637a | 52 | /* if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){ |
| 460c2a | 53 | vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| 41637a | 54 | // continue; |
| SP | 55 | break; |
| 56 | } */ | |
| 460c2a | 57 | // All checks OK! |
| SP | 58 | |
| 59 | dvol=0.0; | |
| 60 | for(j=0;j<vtx_moved->tristar_no;j++){ | |
| 61 | dvol-=vtx_moved->tristar[j]->volume; | |
| fe5069 | 62 | } |
| SP | 63 | volFirst=dvol; |
| 64 | for(j=0;j<vtx_moved->tristar_no;j++){ | |
| 460c2a | 65 | triangle_normal_vector(vtx_moved->tristar[j]); |
| SP | 66 | dvol+=vtx_moved->tristar[j]->volume; |
| 67 | } | |
| fe5069 | 68 | //TODO: here there is a bug. Don't know where, but preliminary success can |
| SP | 69 | //happen sometimes. And when I do this checks, constant value is not achieved |
| 70 | //anymore | |
| 71 | /* voldiff=dvol-Vol; | |
| 460c2a | 72 | if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){ |
| SP | 73 | //calculate energy, return change in energy... |
| 74 | oenergy=vtx_moved->energy; | |
| 75 | energy_vertex(vtx_moved); | |
| 76 | delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy); | |
| 77 | //the same is done for neighbouring vertices | |
| fe5069 | 78 | for(j=0;j<vtx_moved->neigh_no;j++){ |
| SP | 79 | oenergy=vtx_moved->neigh[j]->energy; |
| 80 | energy_vertex(vtx_moved->neigh[j]); | |
| 81 | delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy); | |
| 460c2a | 82 | } |
| SP | 83 | *retEnergy=delta_energy; |
| fbcbdf | 84 | *vtx_backup=backupvtx; |
| SP | 85 | *vtx_moved_retval=vtx_moved; |
| b2fa8c | 86 | fprintf(stderr, "Preliminary success\n"); |
| 460c2a | 87 | return TS_SUCCESS; |
| fe5069 | 88 | } */ |
| SP | 89 | |
| 90 | // fprintf(stderr, "Step 2 success\n"); | |
| 460c2a | 91 | //do it again ;) |
| SP | 92 | dh=Vol*dh/dvol; |
| 93 | vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex)); | |
| 94 | vtx_moved->x=vtx_moved->x*(1-dh/Rv); | |
| 95 | vtx_moved->y=vtx_moved->y*(1-dh/Rv); | |
| 96 | vtx_moved->z=vtx_moved->z*(1-dh/Rv); | |
| 97 | //check for constraints | |
| 98 | if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){ | |
| 99 | vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex)); | |
| fe5069 | 100 | //also, restore normals |
| SP | 101 | for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]); |
| 41637a | 102 | // continue; |
| SP | 103 | break; |
| 460c2a | 104 | } |
| SP | 105 | |
| fe5069 | 106 | dvol=volFirst; |
| b2fa8c | 107 | for(j=0;j<vtx_moved->tristar_no;j++){ |
| SP | 108 | triangle_normal_vector(vtx_moved->tristar[j]); |
| 109 | dvol+=vtx_moved->tristar[j]->volume; | |
| 110 | } | |
| 460c2a | 111 | voldiff=dvol-Vol; |
| SP | 112 | if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){ |
| 113 | //calculate energy, return change in energy... | |
| fe5069 | 114 | // fprintf(stderr, "Constvol success! %e\n",voldiff); |
| 3de289 | 115 | for(j=0;j<vtx_moved->neigh_no;j++){ |
| SP | 116 | memcpy((void *)&backupvtx[j+1],(void *)vtx_moved->neigh[j],sizeof(ts_vertex)); |
| 117 | } | |
| 118 | ||
| 460c2a | 119 | oenergy=vtx_moved->energy; |
| SP | 120 | energy_vertex(vtx_moved); |
| 121 | delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy); | |
| 122 | //the same is done for neighbouring vertices | |
| fe5069 | 123 | for(j=0;j<vtx_moved->neigh_no;j++){ |
| SP | 124 | oenergy=vtx_moved->neigh[j]->energy; |
| 125 | energy_vertex(vtx_moved->neigh[j]); | |
| 126 | delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy); | |
| 460c2a | 127 | } |
| SP | 128 | *retEnergy=delta_energy; |
| fbcbdf | 129 | *vtx_backup=backupvtx; |
| SP | 130 | *vtx_moved_retval=vtx_moved; |
| 460c2a | 131 | return TS_SUCCESS; |
| SP | 132 | } |
| 133 | } | |
| 134 | free(backupvtx); | |
| 135 | return TS_FAIL; | |
| 136 | } | |
| 137 | ||
| 138 | ||
| 139 | ts_bool constvolConstraintCheck(ts_vesicle *vesicle, ts_vertex *vtx){ | |
| 140 | ts_uint i; | |
| 141 | ts_double dist; | |
| 142 | ts_uint cellidx; | |
| 143 | //distance with neighbours check | |
| 144 | for(i=0;i<vtx->neigh_no;i++){ | |
| 145 | dist=vtx_distance_sq(vtx,vtx->neigh[i]); | |
| 146 | if(dist<1.0 || dist>vesicle->dmax) { | |
| 147 | return TS_FAIL; | |
| 148 | } | |
| 149 | } | |
| 150 | // Distance with grafted poly-vertex check: | |
| 151 | if(vtx->grafted_poly!=NULL){ | |
| 152 | dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]); | |
| 153 | if(dist<1.0 || dist>vesicle->dmax) { | |
| 154 | return TS_FAIL; | |
| 155 | } | |
| 156 | } | |
| 157 | ||
| 158 | // Nucleus penetration check: | |
| 159 | if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){ | |
| 160 | return TS_FAIL; | |
| 161 | } | |
| 162 | ||
| 163 | //self avoidance check with distant vertices | |
| 164 | cellidx=vertex_self_avoidance(vesicle, vtx); | |
| 165 | //check occupation number | |
| 166 | return cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx); | |
| 167 | } | |
| 168 | ||
| 169 | ||
| 170 | ||
| 171 | ts_bool constvolumerestore(ts_vertex *vtx_moved,ts_vertex *vtx_backup){ | |
| 172 | ts_uint j; | |
| fbcbdf | 173 | memcpy((void *)vtx_moved,(void *)&vtx_backup[0],sizeof(ts_vertex)); |
| fe5069 | 174 | for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]); |
| 3de289 | 175 | for(j=0;j<vtx_moved->neigh_no;j++){ |
| SP | 176 | // memcpy((void *)vtx_moved->neigh[j],(void *)&vtx_backup[j+1],sizeof(ts_vertex)); |
| 177 | energy_vertex(vtx_moved->neigh[j]); | |
| 178 | } | |
| fe5069 | 179 | |
| fbcbdf | 180 | free(vtx_backup); |
| 460c2a | 181 | return TS_SUCCESS; |
| SP | 182 | } |
| 183 | ||
| fbcbdf | 184 | ts_bool constvolumeaccept(ts_vesicle *vesicle,ts_vertex *vtx_moved, ts_vertex *vtx_backup){ |
| SP | 185 | ts_bool retval; |
| 186 | ts_uint cellidx=vertex_self_avoidance(vesicle, vtx_moved); | |
| 187 | if(vtx_moved->cell!=vesicle->clist->cell[cellidx]){ | |
| 188 | retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx_moved); | |
| 189 | if(retval==TS_SUCCESS) cell_remove_vertex(vtx_backup[0].cell,vtx_moved); | |
| 190 | ||
| 191 | } | |
| 192 | free(vtx_backup); | |
| 460c2a | 193 | |
| SP | 194 | return TS_SUCCESS; |
| 195 | } | |
