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