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