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 |
} |