| 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 | ||
| 9 | ts_uint i, jj, jjp; | |
| 10 | ||
| 11 | ts_vertex_list *vlist=&vesicle->vlist; | |
| 12 | ts_vertex *vtx=vlist->vertex; | |
| 13 | ||
| 14 | for(i=0;i<vlist->n;i++){ | |
| 15 | //should call with zero index!!! | |
| 16 | energy_vertex(&vtx[i]); | |
| 17 | } | |
| 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; | |
| 26 | ts_uint jj; | |
| 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; | |
| 33 | for(jj=1; jj<=vtx->neigh_no;jj++){ | |
| 34 | jjp=jj+1; | |
| 35 | if(jjp>vtx->neigh_no) jjp=1; | |
| 36 | jjm=jj-1; | |
| 37 | if(jjm<1) jjm=vtx->neigh_no; | |
| 38 | j=vtx->neigh[jj-1]; | |
| 39 | jp=vtx->neigh[jjp-1]; | |
| 40 | jm=vtx->neigh[jjm-1]; | |
| 41 | jt=vtx->tristar[jj-1]; | |
| 42 | x1=vertex_distance_sq(vtx,jp); //shouldn't be zero! | |
| 43 | x2=vertex_distance_sq(j,jp); // shouldn't be zero! | |
| 44 | x3=(j->x-jp->x)*(vtx->x-jp->x)+ | |
| 45 | (j->y-jp->y)*(vtx->y-jp->y)+ | |
| 46 | (j->z-jp->z)*(vtx->z-jp->z); | |
| 47 | ||
| 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 | |
| 57 | x1=vertex_distance_sq(vtx,jm); | |
| 58 | x2=vertex_distance_sq(j,jm); | |
| 59 | x3=(j->x-jm->x)*(vtx->x-jm->x)+ | |
| 60 | (j->y-jm->y)*(vtx->y-jm->y)+ | |
| 61 | (j->z-jm->z)*(vtx->z-jm->z); | |
| 62 | #ifdef TS_DOUBLE_DOUBLE | |
| 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; | |
| 73 | xlen=vertex_distance_sq(j,vtx); | |
| 74 | #ifdef TS_DOUBLE_DOUBLE | |
| 75 | vtx->bond_length[jj-1]=sqrt(xlen); | |
| 76 | #endif | |
| 77 | #ifdef TS_DOUBLE_FLOAT | |
| 78 | vtx->bond_length[jj-1]=sqrtf(xlen); | |
| 79 | #endif | |
| 80 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 81 | vtx->bond_length[jj-1]=sqrtl(xlen); | |
| 82 | #endif | |
| 83 | ||
| 84 | vtx->bond_length_dual[jj-1]=tot*vtx->bond_length[jj-1]; | |
| 85 | ||
| 86 | s+=tot*xlen; | |
| 87 | xh+=tot*(j->x - vtx->x); | |
| 88 | yh+=tot*(j->y - vtx->y); | |
| 89 | zh+=tot*(j->z - vtx->z); | |
| 90 | txn+=jt->xnorm; | |
| 91 | tyn+=jt->ynorm; | |
| 92 | tzn+=jt->znorm; | |
| 93 | } | |
| 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) { | |
| 100 | vtx->curvature=sqrt(h); | |
| 101 | } else { | |
| 102 | vtx->curvature=-sqrt(h); | |
| 103 | } | |
| 104 | #endif | |
| 105 | #ifdef TS_DOUBLE_FLOAT | |
| 106 | if(ht>=0.0) { | |
| 107 | vtx->curvature=sqrtf(h); | |
| 108 | } else { | |
| 109 | vtx->curvature=-sqrtf(h); | |
| 110 | } | |
| 111 | #endif | |
| 112 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 113 | if(ht>=0.0) { | |
| 114 | vtx->curvature=sqrtl(h); | |
| 115 | } else { | |
| 116 | vtx->curvature=-sqrtl(h); | |
| 117 | } | |
| 118 | #endif | |
| 119 | vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c); | |
| 120 | ||
| 121 | return TS_SUCCESS; | |
| 122 | } | |
