From 523bf18206f550a315c6c17e5a0a253381b0f8bf Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Thu, 07 Jun 2012 11:16:16 +0000
Subject: [PATCH] Spherical harmonics. Almost everyhing is done. Missing triangle area calculation when vertex is moved or bond is flipped. Also missing volume calculation on vertex move or bondflip. Calculation of co coefficient is not done completely yet. Problems are in numbering the coefficients. Newly added data structure ts_spharm is referenced from ts_vesicle. Missing function for initialization and freeing the memory of that datastructure -- but that memory is already used by some functions.
---
src/energy.c | 33 +++++++++++++++++++--------------
1 files changed, 19 insertions(+), 14 deletions(-)
diff --git a/src/energy.c b/src/energy.c
index 61720ee..4d93753 100644
--- a/src/energy.c
+++ b/src/energy.c
@@ -6,14 +6,14 @@
#include<stdio.h>
ts_bool mean_curvature_and_energy(ts_vesicle *vesicle){
- ts_uint i, jj, jjp;
+ ts_uint i;
- ts_vertex_list *vlist=&vesicle->vlist;
- ts_vertex *vtx=vlist->vertex;
+ ts_vertex_list *vlist=vesicle->vlist;
+ ts_vertex **vtx=vlist->vtx;
for(i=0;i<vlist->n;i++){
- //should call with zero index!!!
- energy_vertex(&vtx[i]);
+ energy_vertex(vtx[i]);
+
}
return TS_SUCCESS;
@@ -23,6 +23,7 @@
inline ts_bool energy_vertex(ts_vertex *vtx){
// ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value!
// ts_triangle *tristar=vtx->tristar-1;
+ //ts_vertex_data *data=vtx->data;
ts_uint jj;
ts_uint jjp,jjm;
ts_vertex *j,*jp, *jm;
@@ -38,9 +39,10 @@
j=vtx->neigh[jj-1];
jp=vtx->neigh[jjp-1];
jm=vtx->neigh[jjm-1];
+// printf("tristar_no=%u, neigh_no=%u, jj=%u\n",data->tristar_no,data->neigh_no,jj);
jt=vtx->tristar[jj-1];
- x1=vertex_distance_sq(vtx,jp); //shouldn't be zero!
- x2=vertex_distance_sq(j,jp); // shouldn't be zero!
+ x1=vtx_distance_sq(vtx,jp); //shouldn't be zero!
+ x2=vtx_distance_sq(j,jp); // shouldn't be zero!
x3=(j->x-jp->x)*(vtx->x-jp->x)+
(j->y-jp->y)*(vtx->y-jp->y)+
(j->z-jp->z)*(vtx->z-jp->z);
@@ -54,8 +56,8 @@
#ifdef TS_DOUBLE_LONGDOUBLE
ctp=x3/sqrtl(x1*x2-x3*x3);
#endif
- x1=vertex_distance_sq(vtx,jm);
- x2=vertex_distance_sq(j,jm);
+ x1=vtx_distance_sq(vtx,jm);
+ x2=vtx_distance_sq(j,jm);
x3=(j->x-jm->x)*(vtx->x-jm->x)+
(j->y-jm->y)*(vtx->y-jm->y)+
(j->z-jm->z)*(vtx->z-jm->z);
@@ -70,18 +72,18 @@
#endif
tot=ctp+ctm;
tot=0.5*tot;
- xlen=vertex_distance_sq(j,vtx);
+ xlen=vtx_distance_sq(j,vtx);
#ifdef TS_DOUBLE_DOUBLE
- vtx->bond_length[jj-1]=sqrt(xlen);
+ vtx->bond[jj-1]->bond_length=sqrt(xlen);
#endif
#ifdef TS_DOUBLE_FLOAT
- vtx->bond_length[jj-1]=sqrtf(xlen);
+ vtx->bond[jj-1]->bond_length=sqrtf(xlen);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
- vtx->bond_length[jj-1]=sqrtl(xlen);
+ vtx->bond[jj-1]->bond_length=sqrtl(xlen);
#endif
- vtx->bond_length_dual[jj-1]=tot*vtx->bond_length[jj-1];
+ vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length;
s+=tot*xlen;
xh+=tot*(j->x - vtx->x);
@@ -116,6 +118,9 @@
vtx->curvature=-sqrtl(h);
}
#endif
+// What is vtx->c?????????????? Here it is 0!
+// c is forced curvature energy for each vertex. Should be set to zero for
+// normal circumstances.
vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c);
return TS_SUCCESS;
--
Gitblit v1.9.3