| 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]); |
| d7639a | 16 | } |
| SP | 17 | |
| 18 | return TS_SUCCESS; | |
| 19 | } | |
| 20 | ||
| 21 | ||
| 22 | inline ts_bool energy_vertex(ts_vertex *vtx){ | |
| 23 | // ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value! | |
| 24 | // ts_triangle *tristar=vtx->tristar-1; | |
| f74313 | 25 | ts_vertex_data *data=vtx->data; |
| d7639a | 26 | ts_uint jj; |
| SP | 27 | ts_uint jjp,jjm; |
| 28 | ts_vertex *j,*jp, *jm; | |
| 29 | ts_triangle *jt; | |
| 30 | ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0; | |
| 31 | ts_double x1,x2,x3,ctp,ctm,tot,xlen; | |
| 32 | ts_double h,ht; | |
| f74313 | 33 | for(jj=1; jj<=data->neigh_no;jj++){ |
| d7639a | 34 | jjp=jj+1; |
| f74313 | 35 | if(jjp>data->neigh_no) jjp=1; |
| d7639a | 36 | jjm=jj-1; |
| f74313 | 37 | if(jjm<1) jjm=data->neigh_no; |
| SP | 38 | j=data->neigh[jj-1]; |
| 39 | jp=data->neigh[jjp-1]; | |
| 40 | jm=data->neigh[jjm-1]; | |
| 41 | jt=data->tristar[jj-1]; | |
| 42 | x1=vtx_distance_sq(vtx,jp); //shouldn't be zero! | |
| 43 | x2=vtx_distance_sq(j,jp); // shouldn't be zero! | |
| 44 | x3=(j->data->x-jp->data->x)*(data->x-jp->data->x)+ | |
| 45 | (j->data->y-jp->data->y)*(data->y-jp->data->y)+ | |
| 46 | (j->data->z-jp->data->z)*(data->z-jp->data->z); | |
| d7639a | 47 | |
| SP | 48 | #ifdef TS_DOUBLE_DOUBLE |
| 49 | ctp=x3/sqrt(x1*x2-x3*x3); | |
| 50 | #endif | |
| 51 | #ifdef TS_DOUBLE_FLOAT | |
| 52 | ctp=x3/sqrtf(x1*x2-x3*x3); | |
| 53 | #endif | |
| 54 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 55 | ctp=x3/sqrtl(x1*x2-x3*x3); | |
| 56 | #endif | |
| f74313 | 57 | x1=vtx_distance_sq(vtx,jm); |
| SP | 58 | x2=vtx_distance_sq(j,jm); |
| 59 | x3=(j->data->x-jm->data->x)*(data->x-jm->data->x)+ | |
| 60 | (j->data->y-jm->data->y)*(data->y-jm->data->y)+ | |
| 61 | (j->data->z-jm->data->z)*(data->z-jm->data->z); | |
| d7639a | 62 | #ifdef TS_DOUBLE_DOUBLE |
| SP | 63 | ctm=x3/sqrt(x1*x2-x3*x3); |
| 64 | #endif | |
| 65 | #ifdef TS_DOUBLE_FLOAT | |
| 66 | ctm=x3/sqrtf(x1*x2-x3*x3); | |
| 67 | #endif | |
| 68 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 69 | ctm=x3/sqrtl(x1*x2-x3*x3); | |
| 70 | #endif | |
| 71 | tot=ctp+ctm; | |
| 72 | tot=0.5*tot; | |
| f74313 | 73 | xlen=vtx_distance_sq(j,vtx); |
| d7639a | 74 | #ifdef TS_DOUBLE_DOUBLE |
| f74313 | 75 | data->bond[jj-1]->data->bond_length=sqrt(xlen); |
| d7639a | 76 | #endif |
| SP | 77 | #ifdef TS_DOUBLE_FLOAT |
| f74313 | 78 | data->bond[jj-1]->data->bond_length=sqrtf(xlen); |
| d7639a | 79 | #endif |
| SP | 80 | #ifdef TS_DOUBLE_LONGDOUBLE |
| f74313 | 81 | data->bond[jj-1]->data->bond_length=sqrtl(xlen); |
| d7639a | 82 | #endif |
| SP | 83 | |
| f74313 | 84 | data->bond[jj-1]->data->bond_length_dual=tot*data->bond[jj-1]->data->bond_length; |
| d7639a | 85 | |
| SP | 86 | s+=tot*xlen; |
| f74313 | 87 | xh+=tot*(j->data->x - data->x); |
| SP | 88 | yh+=tot*(j->data->y - data->y); |
| 89 | zh+=tot*(j->data->z - data->z); | |
| 90 | txn+=jt->data->xnorm; | |
| 91 | tyn+=jt->data->ynorm; | |
| 92 | tzn+=jt->data->znorm; | |
| d7639a | 93 | } |
| SP | 94 | |
| 95 | h=xh*xh+yh*yh+zh*zh; | |
| 96 | ht=txn*xh+tyn*yh + tzn*zh; | |
| 97 | s=s/4.0; | |
| 98 | #ifdef TS_DOUBLE_DOUBLE | |
| 99 | if(ht>=0.0) { | |
| f74313 | 100 | data->curvature=sqrt(h); |
| d7639a | 101 | } else { |
| f74313 | 102 | data->curvature=-sqrt(h); |
| d7639a | 103 | } |
| SP | 104 | #endif |
| 105 | #ifdef TS_DOUBLE_FLOAT | |
| 106 | if(ht>=0.0) { | |
| f74313 | 107 | data->curvature=sqrtf(h); |
| d7639a | 108 | } else { |
| f74313 | 109 | data->curvature=-sqrtf(h); |
| d7639a | 110 | } |
| SP | 111 | #endif |
| 112 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 113 | if(ht>=0.0) { | |
| f74313 | 114 | data->curvature=sqrtl(h); |
| d7639a | 115 | } else { |
| f74313 | 116 | data->curvature=-sqrtl(h); |
| d7639a | 117 | } |
| SP | 118 | #endif |
| f74313 | 119 | //TODO: MAJOR!!!! What is vtx->data->c?????????????? Here it is 0! |
| SP | 120 | data->energy=0.5*s*(data->curvature/s-data->c)*(data->curvature/s-data->c); |
| d7639a | 121 | |
| SP | 122 | return TS_SUCCESS; |
| 123 | } | |
