From 83b03a38e0ab88224d3f8ce9ddfbbdd82038455c Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@gmail.com> Date: Thu, 23 Feb 2012 21:55:40 +0000 Subject: [PATCH] Seems to work. Lookout for debug code in main.c and initial_distribution.c --- src/energy.c | 60 ++++++++++++++++++++++++++++++------------------------------ 1 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/energy.c b/src/energy.c index 7cd96f5..4d93753 100644 --- a/src/energy.c +++ b/src/energy.c @@ -23,7 +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_vertex_data *data=vtx->data; ts_uint jj; ts_uint jjp,jjm; ts_vertex *j,*jp, *jm; @@ -31,21 +31,21 @@ ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=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 @@ -74,21 +74,21 @@ tot=0.5*tot; xlen=vtx_distance_sq(j,vtx); #ifdef TS_DOUBLE_DOUBLE - data->bond[jj-1]->bond_length=sqrt(xlen); + vtx->bond[jj-1]->bond_length=sqrt(xlen); #endif #ifdef TS_DOUBLE_FLOAT - data->bond[jj-1]->bond_length=sqrtf(xlen); + vtx->bond[jj-1]->bond_length=sqrtf(xlen); #endif #ifdef TS_DOUBLE_LONGDOUBLE - data->bond[jj-1]->bond_length=sqrtl(xlen); + vtx->bond[jj-1]->bond_length=sqrtl(xlen); #endif - data->bond[jj-1]->bond_length_dual=tot*data->bond[jj-1]->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); + 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; @@ -99,29 +99,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 -// What is vtx->data->c?????????????? Here it is 0! +// What is vtx->c?????????????? Here it is 0! // c is forced curvature energy for each vertex. Should be set to zero for -// norman circumstances. - data->energy=0.5*s*(data->curvature/s-data->c)*(data->curvature/s-data->c); +// normal circumstances. + vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c); return TS_SUCCESS; } -- Gitblit v1.9.3