| 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; |
| SP | 15 | ts_uint Ntries=20; |
| 16 | ts_vertex *backupvtx; | |
| 17 | ts_double Rv, dh, dvol, voldiff, oenergy,delta_energy; | |
| 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]; | |
| 23 | if(vtx_moved==vtx_avoid) continue; | |
| 24 | ||
| 25 | for(j=0;j<vtx_moved->neigh_no;j++){ | |
| 26 | if(vtx_moved->neigh[j]==vtx_avoid) continue; | |
| b2fa8c | 27 | /* for(k=0;k<vtx_moved->neigh[j]->neigh_no;k++){ |
| SP | 28 | if(vtx_moved->neigh[j]->neigh[k]==vtx_avoid) continue; |
| 29 | } | |
| 30 | */ | |
| 31 | ||
| 460c2a | 32 | } |
| fbcbdf | 33 | |
| 460c2a | 34 | memcpy((void *)&backupvtx[0],(void *)vtx_moved,sizeof(ts_vertex)); |
| SP | 35 | //move vertex in specified direction. first try, test move! |
| 36 | ||
| 37 | Rv=sqrt(pow(vtx_moved->x,2)+pow(vtx_moved->y,2)+pow(vtx_moved->z,2)); | |
| fbcbdf | 38 | dh=2.0*Vol/(sqrt(3.0)*l0*l0); |
| SP | 39 | // fprintf(stderr,"Prej (x,y,z)=(%e,%e,%e).\n",vtx_moved->x,vtx_moved->y,vtx_moved->z); |
| 40 | vtx_moved->x=vtx_moved->x*(1.0-dh/Rv); | |
| 41 | vtx_moved->y=vtx_moved->y*(1.0-dh/Rv); | |
| 42 | vtx_moved->z=vtx_moved->z*(1.0-dh/Rv); | |
| 43 | // fprintf(stderr,"Potem (x,y,z)=(%e,%e,%e). Vol=%e\n",vtx_moved->x,vtx_moved->y,vtx_moved->z,Vol); | |
| 460c2a | 44 | |
| SP | 45 | //check for constraints |
| 46 | if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){ | |
| 47 | vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex)); | |
| 48 | continue; | |
| 49 | } | |
| fbcbdf | 50 | // fprintf(stderr,"Sprejet.\n"); |
| 460c2a | 51 | |
| SP | 52 | // All checks OK! |
| b2fa8c | 53 | fprintf(stderr, "Step 1 success\n"); |
| 460c2a | 54 | |
| SP | 55 | for(j=0;j<vtx_moved->neigh_no;j++){ |
| 56 | memcpy((void *)&backupvtx[j+1],(void *)vtx_moved->neigh[j],sizeof(ts_vertex)); | |
| 57 | } | |
| 58 | dvol=0.0; | |
| 59 | for(j=0;j<vtx_moved->tristar_no;j++){ | |
| 60 | dvol-=vtx_moved->tristar[j]->volume; | |
| 61 | triangle_normal_vector(vtx_moved->tristar[j]); | |
| 62 | dvol+=vtx_moved->tristar[j]->volume; | |
| 63 | } | |
| 64 | ||
| 65 | voldiff=dvol-Vol; | |
| 66 | ||
| 67 | if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){ | |
| 68 | //calculate energy, return change in energy... | |
| 69 | oenergy=vtx_moved->energy; | |
| 70 | energy_vertex(vtx_moved); | |
| 71 | delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy); | |
| 72 | //the same is done for neighbouring vertices | |
| 73 | for(i=0;i<vtx_moved->neigh_no;i++){ | |
| 74 | oenergy=vtx_moved->neigh[i]->energy; | |
| 75 | energy_vertex(vtx_moved->neigh[i]); | |
| 76 | delta_energy+=vtx_moved->neigh[i]->xk*(vtx_moved->neigh[i]->energy-oenergy); | |
| 77 | } | |
| 78 | *retEnergy=delta_energy; | |
| fbcbdf | 79 | *vtx_backup=backupvtx; |
| SP | 80 | *vtx_moved_retval=vtx_moved; |
| b2fa8c | 81 | fprintf(stderr, "Preliminary success\n"); |
| 460c2a | 82 | return TS_SUCCESS; |
| SP | 83 | } |
| b2fa8c | 84 | fprintf(stderr, "Step 2 success\n"); |
| 460c2a | 85 | //do it again ;) |
| SP | 86 | dh=Vol*dh/dvol; |
| 87 | vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex)); | |
| 88 | vtx_moved->x=vtx_moved->x*(1-dh/Rv); | |
| 89 | vtx_moved->y=vtx_moved->y*(1-dh/Rv); | |
| 90 | vtx_moved->z=vtx_moved->z*(1-dh/Rv); | |
| 91 | //check for constraints | |
| 92 | if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){ | |
| 93 | for(j=0;j<vtx_moved->neigh_no;j++){ | |
| 94 | memcpy((void *)vtx_moved->neigh[j],(void *)&backupvtx[j+1],sizeof(ts_vertex)); | |
| 95 | } | |
| 96 | vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex)); | |
| 97 | continue; | |
| 98 | } | |
| 99 | ||
| b2fa8c | 100 | dvol=0.0; |
| SP | 101 | for(j=0;j<vtx_moved->tristar_no;j++){ |
| 102 | dvol-=vtx_moved->tristar[j]->volume; | |
| 103 | triangle_normal_vector(vtx_moved->tristar[j]); | |
| 104 | dvol+=vtx_moved->tristar[j]->volume; | |
| 105 | } | |
| 106 | ||
| 107 | fprintf(stderr, "Step 3a success voldiff=%e\n",voldiff); | |
| 460c2a | 108 | voldiff=dvol-Vol; |
| b2fa8c | 109 | fprintf(stderr, "Step 3b success voldiff=%e\n",voldiff); |
| 460c2a | 110 | if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){ |
| SP | 111 | //calculate energy, return change in energy... |
| 112 | oenergy=vtx_moved->energy; | |
| 113 | energy_vertex(vtx_moved); | |
| 114 | delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy); | |
| 115 | //the same is done for neighbouring vertices | |
| 116 | for(i=0;i<vtx_moved->neigh_no;i++){ | |
| 117 | oenergy=vtx_moved->neigh[i]->energy; | |
| 118 | energy_vertex(vtx_moved->neigh[i]); | |
| 119 | delta_energy+=vtx_moved->neigh[i]->xk*(vtx_moved->neigh[i]->energy-oenergy); | |
| 120 | } | |
| 121 | *retEnergy=delta_energy; | |
| fbcbdf | 122 | *vtx_backup=backupvtx; |
| SP | 123 | *vtx_moved_retval=vtx_moved; |
| b2fa8c | 124 | fprintf(stderr, "DVOL=%e\n",voldiff); |
| 460c2a | 125 | return TS_SUCCESS; |
| SP | 126 | } |
| 127 | ||
| 128 | ||
| 129 | } | |
| 130 | free(backupvtx); | |
| b2fa8c | 131 | fprintf(stderr, "fail\n"); |
| 460c2a | 132 | return TS_FAIL; |
| SP | 133 | } |
| 134 | ||
| 135 | ||
| 136 | ts_bool constvolConstraintCheck(ts_vesicle *vesicle, ts_vertex *vtx){ | |
| 137 | ts_uint i; | |
| 138 | ts_double dist; | |
| 139 | ts_uint cellidx; | |
| 140 | //distance with neighbours check | |
| 141 | for(i=0;i<vtx->neigh_no;i++){ | |
| 142 | dist=vtx_distance_sq(vtx,vtx->neigh[i]); | |
| 143 | if(dist<1.0 || dist>vesicle->dmax) { | |
| 144 | return TS_FAIL; | |
| 145 | } | |
| 146 | } | |
| 147 | // Distance with grafted poly-vertex check: | |
| 148 | if(vtx->grafted_poly!=NULL){ | |
| 149 | dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]); | |
| 150 | if(dist<1.0 || dist>vesicle->dmax) { | |
| 151 | return TS_FAIL; | |
| 152 | } | |
| 153 | } | |
| 154 | ||
| 155 | // Nucleus penetration check: | |
| 156 | if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){ | |
| 157 | return TS_FAIL; | |
| 158 | } | |
| 159 | ||
| 160 | //self avoidance check with distant vertices | |
| 161 | cellidx=vertex_self_avoidance(vesicle, vtx); | |
| 162 | //check occupation number | |
| 163 | return cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx); | |
| 164 | } | |
| 165 | ||
| 166 | ||
| 167 | ||
| 168 | ts_bool constvolumerestore(ts_vertex *vtx_moved,ts_vertex *vtx_backup){ | |
| 169 | ts_uint j; | |
| fbcbdf | 170 | memcpy((void *)vtx_moved,(void *)&vtx_backup[0],sizeof(ts_vertex)); |
| SP | 171 | for(j=0;j<vtx_moved->neigh_no;j++){ |
| 460c2a | 172 | memcpy((void *)vtx_moved->neigh[j],(void *)&vtx_backup[j+1],sizeof(ts_vertex)); |
| fbcbdf | 173 | } |
| SP | 174 | free(vtx_backup); |
| 460c2a | 175 | return TS_SUCCESS; |
| SP | 176 | } |
| 177 | ||
| fbcbdf | 178 | ts_bool constvolumeaccept(ts_vesicle *vesicle,ts_vertex *vtx_moved, ts_vertex *vtx_backup){ |
| SP | 179 | ts_bool retval; |
| 180 | ts_uint cellidx=vertex_self_avoidance(vesicle, vtx_moved); | |
| 181 | if(vtx_moved->cell!=vesicle->clist->cell[cellidx]){ | |
| 182 | retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx_moved); | |
| 183 | if(retval==TS_SUCCESS) cell_remove_vertex(vtx_backup[0].cell,vtx_moved); | |
| 184 | ||
| 185 | } | |
| 186 | free(vtx_backup); | |
| 460c2a | 187 | |
| SP | 188 | return TS_SUCCESS; |
| 189 | } | |
