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