From 142a67fe82b830e5c7816914afa62445959c87ca Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Tue, 05 Nov 2013 14:04:21 +0000
Subject: [PATCH] changes in bondflip call. No need to bondflip all the bonds, but only as many bonds as there are vertices. Also, rnvec seems to be not needed for bondflip, so it is commented out
---
src/energy.c | 72 +++++++++++++++++++-----------------
1 files changed, 38 insertions(+), 34 deletions(-)
diff --git a/src/energy.c b/src/energy.c
index e30295b..abbb8e7 100644
--- a/src/energy.c
+++ b/src/energy.c
@@ -23,29 +23,29 @@
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_vertex_data *data=vtx->data;
ts_uint jj;
ts_uint jjp,jjm;
ts_vertex *j,*jp, *jm;
ts_triangle *jt;
- ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0;
+ ts_double s=0.0,xh=0.0,yh=0.0,zh=0.0,txn=0.0,tyn=0.0,tzn=0.0;
ts_double x1,x2,x3,ctp,ctm,tot,xlen;
ts_double h,ht;
- for(jj=1; jj<=data->neigh_no;jj++){
+ for(jj=1; jj<=vtx->neigh_no;jj++){
jjp=jj+1;
- if(jjp>data->neigh_no) jjp=1;
+ if(jjp>vtx->neigh_no) jjp=1;
jjm=jj-1;
- if(jjm<1) jjm=data->neigh_no;
- j=data->neigh[jj-1];
- jp=data->neigh[jjp-1];
- jm=data->neigh[jjm-1];
+ if(jjm<1) jjm=vtx->neigh_no;
+ 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=data->tristar[jj-1];
+ jt=vtx->tristar[jj-1];
x1=vtx_distance_sq(vtx,jp); //shouldn't be zero!
x2=vtx_distance_sq(j,jp); // shouldn't be zero!
- x3=(j->data->x-jp->data->x)*(data->x-jp->data->x)+
- (j->data->y-jp->data->y)*(data->y-jp->data->y)+
- (j->data->z-jp->data->z)*(data->z-jp->data->z);
+ 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);
#ifdef TS_DOUBLE_DOUBLE
ctp=x3/sqrt(x1*x2-x3*x3);
@@ -58,9 +58,9 @@
#endif
x1=vtx_distance_sq(vtx,jm);
x2=vtx_distance_sq(j,jm);
- x3=(j->data->x-jm->data->x)*(data->x-jm->data->x)+
- (j->data->y-jm->data->y)*(data->y-jm->data->y)+
- (j->data->z-jm->data->z)*(data->z-jm->data->z);
+ 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);
#ifdef TS_DOUBLE_DOUBLE
ctm=x3/sqrt(x1*x2-x3*x3);
#endif
@@ -72,26 +72,28 @@
#endif
tot=ctp+ctm;
tot=0.5*tot;
+
xlen=vtx_distance_sq(j,vtx);
+/*
#ifdef TS_DOUBLE_DOUBLE
- data->bond[jj-1]->data->bond_length=sqrt(xlen);
+ vtx->bond[jj-1]->bond_length=sqrt(xlen);
#endif
#ifdef TS_DOUBLE_FLOAT
- data->bond[jj-1]->data->bond_length=sqrtf(xlen);
+ vtx->bond[jj-1]->bond_length=sqrtf(xlen);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
- data->bond[jj-1]->data->bond_length=sqrtl(xlen);
+ vtx->bond[jj-1]->bond_length=sqrtl(xlen);
#endif
- data->bond[jj-1]->data->bond_length_dual=tot*data->bond[jj-1]->data->bond_length;
-
+ vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length;
+*/
s+=tot*xlen;
- xh+=tot*(j->data->x - data->x);
- yh+=tot*(j->data->y - data->y);
- zh+=tot*(j->data->z - data->z);
- txn+=jt->data->xnorm;
- tyn+=jt->data->ynorm;
- tzn+=jt->data->znorm;
+ xh+=tot*(j->x - vtx->x);
+ yh+=tot*(j->y - vtx->y);
+ zh+=tot*(j->z - vtx->z);
+ txn+=jt->xnorm;
+ tyn+=jt->ynorm;
+ tzn+=jt->znorm;
}
h=xh*xh+yh*yh+zh*zh;
@@ -99,27 +101,29 @@
s=s/4.0;
#ifdef TS_DOUBLE_DOUBLE
if(ht>=0.0) {
- data->curvature=sqrt(h);
+ vtx->curvature=sqrt(h);
} else {
- data->curvature=-sqrt(h);
+ vtx->curvature=-sqrt(h);
}
#endif
#ifdef TS_DOUBLE_FLOAT
if(ht>=0.0) {
- data->curvature=sqrtf(h);
+ vtx->curvature=sqrtf(h);
} else {
- data->curvature=-sqrtf(h);
+ vtx->curvature=-sqrtf(h);
}
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
if(ht>=0.0) {
- data->curvature=sqrtl(h);
+ vtx->curvature=sqrtl(h);
} else {
- data->curvature=-sqrtl(h);
+ vtx->curvature=-sqrtl(h);
}
#endif
-//TODO: MAJOR!!!! What is vtx->data->c?????????????? Here it is 0!
- data->energy=0.5*s*(data->curvature/s-data->c)*(data->curvature/s-data->c);
+// 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