| 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; | |
| f74313 | 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; | |
| f74313 | 34 | for(jj=1; jj<=data->neigh_no;jj++){ |
| d7639a | 35 | jjp=jj+1; |
| f74313 | 36 | if(jjp>data->neigh_no) jjp=1; |
| d7639a | 37 | jjm=jj-1; |
| f74313 | 38 | if(jjm<1) jjm=data->neigh_no; |
| SP | 39 | j=data->neigh[jj-1]; |
| 40 | jp=data->neigh[jjp-1]; | |
| 41 | jm=data->neigh[jjm-1]; | |
| b01cc1 | 42 | // printf("tristar_no=%u, neigh_no=%u, jj=%u\n",data->tristar_no,data->neigh_no,jj); |
| f74313 | 43 | jt=data->tristar[jj-1]; |
| SP | 44 | x1=vtx_distance_sq(vtx,jp); //shouldn't be zero! |
| 45 | x2=vtx_distance_sq(j,jp); // shouldn't be zero! | |
| 46 | x3=(j->data->x-jp->data->x)*(data->x-jp->data->x)+ | |
| 47 | (j->data->y-jp->data->y)*(data->y-jp->data->y)+ | |
| 48 | (j->data->z-jp->data->z)*(data->z-jp->data->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); |
| 61 | x3=(j->data->x-jm->data->x)*(data->x-jm->data->x)+ | |
| 62 | (j->data->y-jm->data->y)*(data->y-jm->data->y)+ | |
| 63 | (j->data->z-jm->data->z)*(data->z-jm->data->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 |
| f74313 | 77 | data->bond[jj-1]->data->bond_length=sqrt(xlen); |
| d7639a | 78 | #endif |
| SP | 79 | #ifdef TS_DOUBLE_FLOAT |
| f74313 | 80 | data->bond[jj-1]->data->bond_length=sqrtf(xlen); |
| d7639a | 81 | #endif |
| SP | 82 | #ifdef TS_DOUBLE_LONGDOUBLE |
| f74313 | 83 | data->bond[jj-1]->data->bond_length=sqrtl(xlen); |
| d7639a | 84 | #endif |
| SP | 85 | |
| f74313 | 86 | data->bond[jj-1]->data->bond_length_dual=tot*data->bond[jj-1]->data->bond_length; |
| d7639a | 87 | |
| SP | 88 | s+=tot*xlen; |
| f74313 | 89 | xh+=tot*(j->data->x - data->x); |
| SP | 90 | yh+=tot*(j->data->y - data->y); |
| 91 | zh+=tot*(j->data->z - data->z); | |
| 92 | txn+=jt->data->xnorm; | |
| 93 | tyn+=jt->data->ynorm; | |
| 94 | tzn+=jt->data->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) { | |
| f74313 | 102 | data->curvature=sqrt(h); |
| d7639a | 103 | } else { |
| f74313 | 104 | data->curvature=-sqrt(h); |
| d7639a | 105 | } |
| SP | 106 | #endif |
| 107 | #ifdef TS_DOUBLE_FLOAT | |
| 108 | if(ht>=0.0) { | |
| f74313 | 109 | data->curvature=sqrtf(h); |
| d7639a | 110 | } else { |
| f74313 | 111 | data->curvature=-sqrtf(h); |
| d7639a | 112 | } |
| SP | 113 | #endif |
| 114 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 115 | if(ht>=0.0) { | |
| f74313 | 116 | data->curvature=sqrtl(h); |
| d7639a | 117 | } else { |
| f74313 | 118 | data->curvature=-sqrtl(h); |
| d7639a | 119 | } |
| SP | 120 | #endif |
| dac2e5 | 121 | // What is vtx->data->c?????????????? Here it is 0! |
| SP | 122 | // c is forced curvature energy for each vertex. Should be set to zero for |
| 123 | // norman circumstances. | |
| f74313 | 124 | data->energy=0.5*s*(data->curvature/s-data->c)*(data->curvature/s-data->c); |
| d7639a | 125 | |
| SP | 126 | return TS_SUCCESS; |
| 127 | } | |
