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