| commit | author | age | ||
| d7639a | 1 | #include<stdlib.h> |
| SP | 2 | #include "general.h" |
| 3 | #include "energy.h" | |
| 4 | #include "vertex.h" | |
| 5 | #include<math.h> | |
| 6 | #include<stdio.h> | |
| 7 | ts_bool mean_curvature_and_energy(ts_vesicle *vesicle){ | |
| 8 | ||
| f74313 | 9 | ts_uint i; |
| d7639a | 10 | |
| f74313 | 11 | ts_vertex_list *vlist=vesicle->vlist; |
| SP | 12 | ts_vertex **vtx=vlist->vtx; |
| d7639a | 13 | |
| SP | 14 | for(i=0;i<vlist->n;i++){ |
| f74313 | 15 | energy_vertex(vtx[i]); |
| b01cc1 | 16 | |
| d7639a | 17 | } |
| SP | 18 | |
| 19 | return TS_SUCCESS; | |
| 20 | } | |
| 21 | ||
| 22 | ||
| 23 | inline ts_bool energy_vertex(ts_vertex *vtx){ | |
| 24 | // ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value! | |
| 25 | // ts_triangle *tristar=vtx->tristar-1; | |
| 8f6a69 | 26 | //ts_vertex_data *data=vtx->data; |
| d7639a | 27 | ts_uint jj; |
| SP | 28 | ts_uint jjp,jjm; |
| 29 | ts_vertex *j,*jp, *jm; | |
| 30 | ts_triangle *jt; | |
| 31 | ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0; | |
| 32 | ts_double x1,x2,x3,ctp,ctm,tot,xlen; | |
| 33 | ts_double h,ht; | |
| 8f6a69 | 34 | for(jj=1; jj<=vtx->neigh_no;jj++){ |
| d7639a | 35 | jjp=jj+1; |
| 8f6a69 | 36 | if(jjp>vtx->neigh_no) jjp=1; |
| d7639a | 37 | jjm=jj-1; |
| 8f6a69 | 38 | if(jjm<1) jjm=vtx->neigh_no; |
| SP | 39 | j=vtx->neigh[jj-1]; |
| 40 | jp=vtx->neigh[jjp-1]; | |
| 41 | jm=vtx->neigh[jjm-1]; | |
| b01cc1 | 42 | // printf("tristar_no=%u, neigh_no=%u, jj=%u\n",data->tristar_no,data->neigh_no,jj); |
| 8f6a69 | 43 | jt=vtx->tristar[jj-1]; |
| f74313 | 44 | x1=vtx_distance_sq(vtx,jp); //shouldn't be zero! |
| SP | 45 | x2=vtx_distance_sq(j,jp); // shouldn't be zero! |
| 8f6a69 | 46 | x3=(j->x-jp->x)*(vtx->x-jp->x)+ |
| SP | 47 | (j->y-jp->y)*(vtx->y-jp->y)+ |
| 48 | (j->z-jp->z)*(vtx->z-jp->z); | |
| d7639a | 49 | |
| SP | 50 | #ifdef TS_DOUBLE_DOUBLE |
| 51 | ctp=x3/sqrt(x1*x2-x3*x3); | |
| 52 | #endif | |
| 53 | #ifdef TS_DOUBLE_FLOAT | |
| 54 | ctp=x3/sqrtf(x1*x2-x3*x3); | |
| 55 | #endif | |
| 56 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 57 | ctp=x3/sqrtl(x1*x2-x3*x3); | |
| 58 | #endif | |
| f74313 | 59 | x1=vtx_distance_sq(vtx,jm); |
| SP | 60 | x2=vtx_distance_sq(j,jm); |
| 8f6a69 | 61 | x3=(j->x-jm->x)*(vtx->x-jm->x)+ |
| SP | 62 | (j->y-jm->y)*(vtx->y-jm->y)+ |
| 63 | (j->z-jm->z)*(vtx->z-jm->z); | |
| d7639a | 64 | #ifdef TS_DOUBLE_DOUBLE |
| SP | 65 | ctm=x3/sqrt(x1*x2-x3*x3); |
| 66 | #endif | |
| 67 | #ifdef TS_DOUBLE_FLOAT | |
| 68 | ctm=x3/sqrtf(x1*x2-x3*x3); | |
| 69 | #endif | |
| 70 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 71 | ctm=x3/sqrtl(x1*x2-x3*x3); | |
| 72 | #endif | |
| 73 | tot=ctp+ctm; | |
| 74 | tot=0.5*tot; | |
| f74313 | 75 | xlen=vtx_distance_sq(j,vtx); |
| d7639a | 76 | #ifdef TS_DOUBLE_DOUBLE |
| 8f6a69 | 77 | vtx->bond[jj-1]->bond_length=sqrt(xlen); |
| d7639a | 78 | #endif |
| SP | 79 | #ifdef TS_DOUBLE_FLOAT |
| 8f6a69 | 80 | vtx->bond[jj-1]->bond_length=sqrtf(xlen); |
| d7639a | 81 | #endif |
| SP | 82 | #ifdef TS_DOUBLE_LONGDOUBLE |
| 8f6a69 | 83 | vtx->bond[jj-1]->bond_length=sqrtl(xlen); |
| d7639a | 84 | #endif |
| SP | 85 | |
| 8f6a69 | 86 | vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length; |
| d7639a | 87 | |
| SP | 88 | s+=tot*xlen; |
| 8f6a69 | 89 | xh+=tot*(j->x - vtx->x); |
| SP | 90 | yh+=tot*(j->y - vtx->y); |
| 91 | zh+=tot*(j->z - vtx->z); | |
| 41a035 | 92 | txn+=jt->xnorm; |
| SP | 93 | tyn+=jt->ynorm; |
| 94 | tzn+=jt->znorm; | |
| d7639a | 95 | } |
| SP | 96 | |
| 97 | h=xh*xh+yh*yh+zh*zh; | |
| 98 | ht=txn*xh+tyn*yh + tzn*zh; | |
| 99 | s=s/4.0; | |
| 100 | #ifdef TS_DOUBLE_DOUBLE | |
| 101 | if(ht>=0.0) { | |
| 8f6a69 | 102 | vtx->curvature=sqrt(h); |
| d7639a | 103 | } else { |
| 8f6a69 | 104 | vtx->curvature=-sqrt(h); |
| d7639a | 105 | } |
| SP | 106 | #endif |
| 107 | #ifdef TS_DOUBLE_FLOAT | |
| 108 | if(ht>=0.0) { | |
| 8f6a69 | 109 | vtx->curvature=sqrtf(h); |
| d7639a | 110 | } else { |
| 8f6a69 | 111 | vtx->curvature=-sqrtf(h); |
| d7639a | 112 | } |
| SP | 113 | #endif |
| 114 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 115 | if(ht>=0.0) { | |
| 8f6a69 | 116 | vtx->curvature=sqrtl(h); |
| d7639a | 117 | } else { |
| 8f6a69 | 118 | vtx->curvature=-sqrtl(h); |
| d7639a | 119 | } |
| SP | 120 | #endif |
| 8f6a69 | 121 | // What is vtx->c?????????????? Here it is 0! |
| dac2e5 | 122 | // c is forced curvature energy for each vertex. Should be set to zero for |
| 8f6a69 | 123 | // normal circumstances. |
| SP | 124 | vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c); |
| d7639a | 125 | |
| SP | 126 | return TS_SUCCESS; |
| 127 | } | |
