vertex data removed and memory leaks fixed

Samo Penic

2012-02-23 8f6a69d8fddcae6b54d25cf02dd56f0afce6a075

vertex data removed and memory leaks fixed

 src/cell.c

@@ -49,9 +49,9 @@
    ts_uint cellidx;
    ts_uint ncx, ncy,ncz;
    ts_cell_list *clist=vesicle->clist;
    ncx=(ts_uint)((vtx->data->x-vesicle->cm[0])*clist->dcell+clist->shift);
    ncy=(ts_uint)((vtx->data->y-vesicle->cm[1])*clist->dcell+clist->shift);
    ncz=(ts_uint)((vtx->data->z-vesicle->cm[2])*clist->dcell+clist->shift);
    ncx=(ts_uint)((vtx->x-vesicle->cm[0])*clist->dcell+clist->shift);
    ncy=(ts_uint)((vtx->y-vesicle->cm[1])*clist->dcell+clist->shift);
    ncz=(ts_uint)((vtx->z-vesicle->cm[2])*clist->dcell+clist->shift);

    if(ncx == clist->ncmax[0]-1 || ncx == 2){
        fatal("Vesicle is positioned outside the cell covered area. Coordinate x is the problem.",1500);

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

 src/frame.c

@@ -9,9 +9,9 @@
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;
    for(i=0;i<n;i++){
        vesicle->cm[0]+=vtx[i]->data->x;
        vesicle->cm[1]+=vtx[i]->data->y;
        vesicle->cm[2]+=vtx[i]->data->z; 
        vesicle->cm[0]+=vtx[i]->x;
        vesicle->cm[1]+=vtx[i]->y;
        vesicle->cm[2]+=vtx[i]->z; 
    } 
    vesicle->cm[0]/=(ts_float)n;
    vesicle->cm[1]/=(ts_float)n;
@@ -30,7 +30,7 @@
    cell_list_cell_occupation_clear(vesicle->clist);
    for(i=0;i<n;i++){
    cellidx=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]);
    vesicle->vlist->vtx[i]->data->cell=vesicle->clist->cell[cellidx];
    vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx];

    cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->vlist->vtx[i]);


 src/general.h

@@ -112,10 +112,10 @@

/** @brief Data structure of all data connected to a vertex
 *
 *  ts_vertex_data holds the data for one single point (bead, vertex). To understand how to use it
 *  ts_vertex holds the data for one single point (bead, vertex). To understand how to use it
 *  here is a detailed description of the fields in the data structure. */
struct ts_vertex_data {
        ts_uint idx; /**< Represents index of the vertex point. Should become obsolete, since it is also present in ts_vertex structure. */        
struct ts_vertex {
        ts_uint idx;
        ts_double x; /**< The x coordinate of vertex. */
        ts_double y; /**< The y coordinate of vertex. */
        ts_double z; /**< The z coordinate of vertex. */
@@ -134,12 +134,6 @@
        ts_double xk;
        ts_double c;
        ts_uint id;
};
typedef struct ts_vertex_data ts_vertex_data;

struct ts_vertex {
        ts_uint idx;
        ts_vertex_data *data;
};
typedef struct ts_vertex ts_vertex;


 src/initial_distribution.c

@@ -64,40 +64,40 @@
    ts_double dx,dy; // end loop prereq

    /* topmost vertex */
    vtx[1]->data->x=0.0;
    vtx[1]->data->y=0.0;
    vtx[1]->data->z=z0*(ts_double)nshell;
    vtx[1]->x=0.0;
    vtx[1]->y=0.0;
    vtx[1]->z=z0*(ts_double)nshell;
    
    /* starting from to in circular order on pentagrams */    
    for(i=1;i<=nshell;i++){
        n0=2+5*i*(i-1)/2; //-1 would be for the reason that C index starts from 0 
        vtx[n0]->data->x=0.0;
        vtx[n0]->data->y=(ts_double)i*xl0;
        vtx[n0+i]->data->x=vtx[n0]->data->y*s1;
        vtx[n0+i]->data->y=vtx[n0]->data->y*c1;
        vtx[n0+2*i]->data->x=vtx[n0]->data->y*s2;
        vtx[n0+2*i]->data->y=vtx[n0]->data->y*c2;
        vtx[n0+3*i]->data->x=-vtx[n0+2*i]->data->x;
        vtx[n0+3*i]->data->y=vtx[n0+2*i]->data->y;
        vtx[n0+4*i]->data->x=-vtx[n0+i]->data->x;
        vtx[n0+4*i]->data->y=vtx[n0+i]->data->y;
        vtx[n0]->x=0.0;
        vtx[n0]->y=(ts_double)i*xl0;
        vtx[n0+i]->x=vtx[n0]->y*s1;
        vtx[n0+i]->y=vtx[n0]->y*c1;
        vtx[n0+2*i]->x=vtx[n0]->y*s2;
        vtx[n0+2*i]->y=vtx[n0]->y*c2;
        vtx[n0+3*i]->x=-vtx[n0+2*i]->x;
        vtx[n0+3*i]->y=vtx[n0+2*i]->y;
        vtx[n0+4*i]->x=-vtx[n0+i]->x;
        vtx[n0+4*i]->y=vtx[n0+i]->y;
    }

    /* vertexes on the faces of the dipyramid */
    for(i=1;i<=nshell;i++){
        n0=2+5*i*(i-1)/2; // -1 would be because of C!
        for(j=1;j<=i-1;j++){
            dx=(vtx[n0]->data->x-vtx[n0+4*i]->data->x)/(ts_double)i;
            dy=(vtx[n0]->data->y-vtx[n0+4*i]->data->y)/(ts_double)i;
            vtx[n0+4*i+j]->data->x=(ts_double)j*dx+vtx[n0+4*i]->data->x;
            vtx[n0+4*i+j]->data->y=(ts_double)j*dy+vtx[n0+4*i]->data->y;
            dx=(vtx[n0]->x-vtx[n0+4*i]->x)/(ts_double)i;
            dy=(vtx[n0]->y-vtx[n0+4*i]->y)/(ts_double)i;
            vtx[n0+4*i+j]->x=(ts_double)j*dx+vtx[n0+4*i]->x;
            vtx[n0+4*i+j]->y=(ts_double)j*dy+vtx[n0+4*i]->y;
        }
        for(k=0;k<=3;k++){ // I would be worried about zero starting of for
            dx=(vtx[n0+(k+1)*i]->data->x - vtx[n0+k*i]->data->x)/(ts_double) i;
            dy=(vtx[n0+(k+1)*i]->data->y - vtx[n0+k*i]->data->y)/(ts_double) i;
            dx=(vtx[n0+(k+1)*i]->x - vtx[n0+k*i]->x)/(ts_double) i;
            dy=(vtx[n0+(k+1)*i]->y - vtx[n0+k*i]->y)/(ts_double) i;
            for(j=1; j<=i-1;j++){
                vtx[n0+k*i+j]->data->x= (ts_double)j*dx+vtx[n0+k*i]->data->x;
                vtx[n0+k*i+j]->data->y= (ts_double)j*dy+vtx[n0+k*i]->data->y;
                vtx[n0+k*i+j]->x= (ts_double)j*dx+vtx[n0+k*i]->x;
                vtx[n0+k*i+j]->y= (ts_double)j*dy+vtx[n0+k*i]->y;
            } 
        } 
    }
@@ -105,15 +105,15 @@
    for(i=1;i<=nshell;i++){
        n0= 2+ 5*i*(i-1)/2;
        for(j=0;j<=5*i-1;j++){
        vtx[n0+j]->data->z= z0*(ts_double)(nshell-i);   // I would be worried about zero starting of for
        vtx[n0+j]->z= z0*(ts_double)(nshell-i);   // I would be worried about zero starting of for
        }
    }

/* for botom part of dipyramide we calculate the positions of vertices */
    for(i=2+5*nshell*(nshell+1)/2;i<=vlist->n;i++){
        vtx[i]->data->x=vtx[vlist->n - i +1]->data->x;
        vtx[i]->data->y=vtx[vlist->n - i +1]->data->y;
        vtx[i]->data->z=-vtx[vlist->n - i +1]->data->z;
        vtx[i]->x=vtx[vlist->n - i +1]->x;
        vtx[i]->y=vtx[vlist->n - i +1]->y;
        vtx[i]->z=-vtx[vlist->n - i +1]->z;
    }

    for(i=1;i<=vlist->n;i++){
@@ -166,15 +166,16 @@
    ts_double direct; // Something, dont know what, but could be normal of some kind
    for(i=1;i<=vlist->n;i++){
        k++; // WHY i IS NOT GOOD??
           vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->data->neigh[0]->idx+1]); //always add 1st
           vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->neigh[0]->idx+1]); //always add 1st
           jjj=1;
           jj=1;
           for(l=2;l<=vtx[i]->data->neigh_no;l++){
               for(j=2;j<=vtx[i]->data->neigh_no;j++){
                   dist2=vtx_distance_sq(vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
                   direct=vtx_direct(vtx[i],vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
           for(l=2;l<=vtx[i]->neigh_no;l++){
               for(j=2;j<=vtx[i]->neigh_no;j++){
                   dist2=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
                   direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
// TODO: check if fabs can be used with all floating point types!!
                   if( (fabs(dist2-A0*A0)<=eps) && (direct>0.0) && (j!=jjj) ){
                       vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->data->neigh[j-1]->idx+1]);
                       vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->neigh[j-1]->idx+1]);
                       jjj=jj;
                       jj=j;
                       break;
@@ -206,8 +207,8 @@
    ts_uint i,j,k;
    for(i=1;i<=vlist->n;i++){
        for(j=i+1;j<=vlist->n;j++){
            for(k=0;k<vtx[i]->data->neigh_no;k++){ // has changed 0 to < instead of 1 and <=
                if(vtx[i]->data->neigh[k]==vtx[j]){  //if addresses matches it is the same
            for(k=0;k<vtx[i]->neigh_no;k++){ // has changed 0 to < instead of 1 and <=
                if(vtx[i]->neigh[k]==vtx[j]){  //if addresses matches it is the same
                    bond_add(blist,vtx[i],vtx[j]);
                    break;
                }
@@ -232,14 +233,15 @@
    ts_double eps=0.001; // can we use EPS from math.h?
    k=0;
    for(i=1;i<=vesicle->vlist->n;i++){
        for(j=1;j<=vtx[i]->data->neigh_no;j++){
            for(jj=1;jj<=vtx[i]->data->neigh_no;jj++){
        for(j=1;j<=vtx[i]->neigh_no;j++){
            for(jj=1;jj<=vtx[i]->neigh_no;jj++){
        //        ts_fprintf(stderr,"%u: (%u,%u) neigh_no=%u ",i,j,jj,vtx[i].neigh_no);
        //      ts_fprintf(stderr,"%e, %e",vtx[i].neigh[j-1]->x,vtx[i].neigh[jj-1]->x);
                dist=vtx_distance_sq(vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
                direct=vtx_direct(vtx[i],vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);				
                if(fabs(dist-A0*A0)<=eps && direct < 0.0 && vtx[i]->data->neigh[j-1]->idx+1 > i && vtx[i]->data->neigh[jj-1]->idx+1 >i){
                    triangle_add(tlist,vtx[i],vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
                dist=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
                direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);				
// TODO: same as above				
                if(fabs(dist-A0*A0)<=eps && direct < 0.0 && vtx[i]->neigh[j-1]->idx+1 > i && vtx[i]->neigh[jj-1]->idx+1 >i){
                    triangle_add(tlist,vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
                }    
            }    
        }
@@ -254,7 +256,7 @@
            k++;
        }
            if(k!=3){
                fatal("Some triangles has less than 3 vertices..",4);
                fatal("Some triangles have less than 3 vertices..",4);
            }   
    } 
    if(tlist->n!=2*(vesicle->vlist->n -2)){
@@ -334,11 +336,11 @@
    ts_triangle **tria=tlist->tria -1;

    for(i=1;i<=vesicle->vlist->n;i++){
        for(j=1;j<=vtx[i]->data->neigh_no;j++){
            k1=vtx[i]->data->neigh[j-1];
        for(j=1;j<=vtx[i]->neigh_no;j++){
            k1=vtx[i]->neigh[j-1];
            jp=j+1;
            if(j == vtx[i]->data->neigh_no) jp=1;
            k2=vtx[i]->data->neigh[jp-1];
            if(j == vtx[i]->neigh_no) jp=1;
            k2=vtx[i]->neigh[jp-1];
            for(k=1;k<=tlist->n;k++){        // VERY NON-OPTIMAL!!! too many loops (vlist.n * vtx.neigh * tlist.n )!
                k3=tria[k]->vertex[0];
                k4=tria[k]->vertex[1];

 src/io.c

@@ -15,7 +15,7 @@
    printf("Number of vertices: %u\n",vlist->n);
    for(i=0;i<vlist->n;i++){
        printf("%u: %f %f %f\n",
vlist->vtx[i]->idx,vlist->vtx[i]->data->x, vlist->vtx[i]->data->y, vlist->vtx[i]->data->z);
vlist->vtx[i]->idx,vlist->vtx[i]->x, vlist->vtx[i]->y, vlist->vtx[i]->z);
    }
    return TS_SUCCESS;
}
@@ -25,10 +25,10 @@
    ts_vertex **vtx=vlist->vtx;
    printf("Vertex id(neigh no): (neighvertex coord) (neighvertex coord) ...\n");
    for(i=0;i<vlist->n;i++){
        printf("%u(%u): ",vtx[i]->idx,vtx[i]->data->neigh_no);
        for(j=0;j<vtx[i]->data->neigh_no;j++){
            printf("(%f,%f,%f)",vtx[i]->data->neigh[j]->data->x,
vtx[i]->data->neigh[j]->data->y,vtx[i]->data->neigh[j]->data->z);
        printf("%u(%u): ",vtx[i]->idx,vtx[i]->neigh_no);
        for(j=0;j<vtx[i]->neigh_no;j++){
            printf("(%f,%f,%f)",vtx[i]->neigh[j]->x,
vtx[i]->neigh[j]->y,vtx[i]->neigh[j]->z);
        }
        printf("\n");
    }
@@ -47,7 +47,7 @@
        return TS_FAIL;
    }
    for(i=0;i<vlist->n;i++)
        fprintf(fh," %E\t%E\t%E\n",vtx[i]->data->x,vtx[i]->data->y, vtx[i]->data->z);
        fprintf(fh," %E\t%E\t%E\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);

    fclose(fh);
return TS_SUCCESS;
@@ -57,11 +57,11 @@
ts_bool fprint_vertex_list(FILE *fh,ts_vertex_list *vlist){
    ts_uint i,j;
    for(i=0;i<vlist->n;i++){
        fprintf(fh," %.17E\t%.17E\t%.17E\t%u\n",vlist->vtx[i]->data->x,
            vlist->vtx[i]->data->y, vlist->vtx[i]->data->z,
            vlist->vtx[i]->data->neigh_no);
        for(j=0;j<vlist->vtx[i]->data->neigh_no;j++){
            fprintf(fh,"\t%u",(ts_uint)(vlist->vtx[i]->data->neigh[j]->idx));
        fprintf(fh," %.17E\t%.17E\t%.17E\t%u\n",vlist->vtx[i]->x,
            vlist->vtx[i]->y, vlist->vtx[i]->z,
            vlist->vtx[i]->neigh_no);
        for(j=0;j<vlist->vtx[i]->neigh_no;j++){
            fprintf(fh,"\t%u",(ts_uint)(vlist->vtx[i]->neigh[j]->idx));
        //-vlist->vtx+1));
        }
        fprintf(fh,"\n");
@@ -72,9 +72,9 @@
ts_bool fprint_tristar(FILE *fh, ts_vesicle *vesicle){
    ts_uint i,j;
    for(i=0;i<vesicle->vlist->n;i++){
        fprintf(fh,"\t%u",vesicle->vlist->vtx[i]->data->tristar_no);
        for(j=0;j<vesicle->vlist->vtx[i]->data->tristar_no;j++){
            fprintf(fh,"\t%u",(ts_uint)(vesicle->vlist->vtx[i]->data->tristar[j]->idx));//-vesicle->tlist->tria+1));
        fprintf(fh,"\t%u",vesicle->vlist->vtx[i]->tristar_no);
        for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
            fprintf(fh,"\t%u",(ts_uint)(vesicle->vlist->vtx[i]->tristar[j]->idx));//-vesicle->tlist->tria+1));
        }
        fprintf(fh,"\n");
    }
@@ -105,14 +105,14 @@
    ts_uint i,j;
    for(i=0;i<vlist->n;i++){
        fprintf(fh," %.17E\t%.17E\t%.17E\t%.17E\t%.17E\t%u\n",
        vlist->vtx[i]->data->xk,vlist->vtx[i]->data->c,vlist->vtx[i]->data->energy,
        vlist->vtx[i]->data->energy_h, vlist->vtx[i]->data->curvature, 0);
        for(j=0;j<vlist->vtx[i]->data->neigh_no;j++){
            fprintf(fh," %.17E", vlist->vtx[i]->data->bond[j]->bond_length_dual);
        vlist->vtx[i]->xk,vlist->vtx[i]->c,vlist->vtx[i]->energy,
        vlist->vtx[i]->energy_h, vlist->vtx[i]->curvature, 0);
        for(j=0;j<vlist->vtx[i]->neigh_no;j++){
            fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length_dual);
        }
            fprintf(fh,"\n");
        for(j=0;j<vlist->vtx[i]->data->neigh_no;j++){
            fprintf(fh," %.17E", vlist->vtx[i]->data->bond[j]->bond_length);
        for(j=0;j<vlist->vtx[i]->neigh_no;j++){
            fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length);
        }
            fprintf(fh,"\n");
    }
@@ -226,7 +226,7 @@

    fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
    for(i=0;i<vlist->n;i++){
        fprintf(fh,"%e %e %e\n",vtx[i]->data->x,vtx[i]->data->y, vtx[i]->data->z);
        fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
    }

    fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
@@ -270,7 +270,7 @@
    fprintf(fh,"DATASET UNSTRUCTURED_GRID\n");
    fprintf(fh,"POINTS %u double\n", vlist->n);
    for(i=0;i<vlist->n;i++){
        fprintf(fh,"%e %e %e\n",vtx[i]->data->x,vtx[i]->data->y, vtx[i]->data->z);
        fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
    }
    
    fprintf(fh,"CELLS %u %u\n",blist->n,3*blist->n);
@@ -286,7 +286,7 @@
    fprintf(fh,"LOOKUP_TABLE default\n");

    for(i=0;i<vlist->n;i++)
        fprintf(fh,"%u\n",vtx[i]->data->idx);
        fprintf(fh,"%u\n",vtx[i]->idx);

    fclose(fh);
    return TS_SUCCESS;
@@ -351,11 +351,11 @@
    vesicle->vlist=init_vertex_list(nvtx);
    vlist=vesicle->vlist;
    for(i=0;i<nvtx;i++){
   //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->data->x,&vlist->vtx[i]->data->y,&vlist->vtx[i]->data->z);
   //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->x,&vlist->vtx[i]->y,&vlist->vtx[i]->z);
       retval=fscanf(fh,"%F %F %F",&x,&y,&z);
        vlist->vtx[i]->data->x=x;
        vlist->vtx[i]->data->y=y;
        vlist->vtx[i]->data->z=z;
        vlist->vtx[i]->x=x;
        vlist->vtx[i]->y=y;
        vlist->vtx[i]->z=z;
    }
    for(i=0;i<nedges;i++){
        retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
@@ -368,8 +368,8 @@
    /*
    for(i=0;i<ntria;i++){
        retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
        vtxi1=vesicle->blist->data->vertex1->idx;
        vtxi2=vesicle->blist->data->vertex1->idx;
        vtxi1=vesicle->blist->vertex1->idx;
        vtxi2=vesicle->blist->vertex1->idx;
        
    }
    */

 src/main.c

@@ -36,8 +36,8 @@
retval=vtx_remove_neighbour(vlist->vtx[0],vlist->vtx[1]);
vtx_add_neighbour(vlist->vtx[0],vlist->vtx[1]);

vlist->vtx[0]->data->x=1.0;
vlist->vtx[0]->data->x=1.1;
vlist->vtx[0]->x=1.0;
vlist->vtx[0]->x=1.1;
vlist1=vertex_list_copy(vlist);
bond_add(blist, vlist->vtx[1],vlist->vtx[0]);
triangle_add(tlist,vlist->vtx[1],vlist->vtx[2],vlist->vtx[3]);
@@ -47,7 +47,7 @@
printf("Cell idx=1 has vertices=%u\n",clist->cell[0]->nvertex);
cell_add_vertex(clist->cell[0], vlist->vtx[0]);
printf("Cell idx=1 has vertices=%u\n",clist->cell[0]->nvertex);
printf("Cell idx=1 has vertex[0] has x coordinate=%e\n",clist->cell[0]->vertex[0]->data->x);
printf("Cell idx=1 has vertex[0] has x coordinate=%e\n",clist->cell[0]->vertex[0]->x);
cell_list_cell_occupation_clear(clist);
printf("Cell idx=1 has vertices=%u\n",clist->cell[0]->nvertex);
cell_add_vertex(clist->cell[0], vlist->vtx[0]);

 src/shdiscover.c

@@ -42,17 +42,17 @@

/*we calculate new position of each vertex of vesicle */
for(i=0;i<n;i++){
    fi=atan2(vlist->vtx[i]->data->y, vlist->vtx[i]->data->x);
    fi=atan2(vlist->vtx[i]->y, vlist->vtx[i]->x);
/*    theta=atan2(
        sqrt(vlist->vtx[i]->data->x*vlist->vtx[i]->data->x + 
        vlist->vtx[i]->data->y*vlist->vtx[i]->data->y),
        vlist->vtx[i]->data->z 
        ); */
    theta=acos(
        vlist->vtx[i]->data->z /
        sqrt(vlist->vtx[i]->data->x*vlist->vtx[i]->data->x + 
        vlist->vtx[i]->data->y*vlist->vtx[i]->data->y+
        vlist->vtx[i]->data->z*vlist->vtx[i]->data->z)
        vlist->vtx[i]->z /
        sqrt(vlist->vtx[i]->x*vlist->vtx[i]->x + 
        vlist->vtx[i]->y*vlist->vtx[i]->y+
        vlist->vtx[i]->z*vlist->vtx[i]->z)

        );

@@ -67,9 +67,9 @@
        /*ts_fprintf(stdout,"l=%d, m=%d, u=%s\n",l,m,argv[j]);*/
    }

    vlist->vtx[i]->data->z=fabs(r)*cos(theta);
    vlist->vtx[i]->data->x=fabs(r)*sin(theta)*cos(fi);
    vlist->vtx[i]->data->y=fabs(r)*sin(theta)*sin(fi);
    vlist->vtx[i]->z=fabs(r)*cos(theta);
    vlist->vtx[i]->x=fabs(r)*sin(theta)*cos(fi);
    vlist->vtx[i]->y=fabs(r)*sin(theta)*sin(fi);
}

write_vertex_xml_file(vesicle,0);

 src/triangle.c

@@ -194,12 +194,12 @@
  */
ts_bool triangle_normal_vector(ts_triangle *tria){
    ts_double x21,x31,y21,y31,z21,z31,xden;
    x21=tria->vertex[1]->data->x - tria->vertex[0]->data->x;
    x31=tria->vertex[2]->data->x - tria->vertex[0]->data->x;
    y21=tria->vertex[1]->data->y - tria->vertex[0]->data->y;
    y31=tria->vertex[2]->data->y - tria->vertex[0]->data->y;
    z21=tria->vertex[1]->data->z - tria->vertex[0]->data->z;
    z31=tria->vertex[2]->data->z - tria->vertex[0]->data->z;
    x21=tria->vertex[1]->x - tria->vertex[0]->x;
    x31=tria->vertex[2]->x - tria->vertex[0]->x;
    y21=tria->vertex[1]->y - tria->vertex[0]->y;
    y31=tria->vertex[2]->y - tria->vertex[0]->y;
    z21=tria->vertex[1]->z - tria->vertex[0]->z;
    z31=tria->vertex[2]->z - tria->vertex[0]->z;

    tria->xnorm=y21*z31 - z21*y31;
    tria->ynorm=z21*x31 - x21*z31;

 src/vertex.c

@@ -8,7 +8,6 @@

ts_vertex_list *init_vertex_list(ts_uint N){    
    ts_int i;
    ts_vertex *tlist;
    ts_vertex_list *vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
    
    if(N==0){
@@ -18,27 +17,16 @@
        return vlist;
    }
    
    vlist->vtx=(ts_vertex **)malloc(N*sizeof(ts_vertex *));
    tlist=(ts_vertex *)malloc(N*sizeof(ts_vertex));
    if(vlist->vtx==NULL || tlist==NULL)
    vlist->vtx=(ts_vertex **)calloc(N,sizeof(ts_vertex *));
    if(vlist->vtx==NULL)
        fatal("Fatal error reserving memory space for vertex list! Could number of requsted vertices be too large?", 100);
    for(i=0;i<N;i++) {
        vlist->vtx[i]=&tlist[i];
        vlist->vtx[i]->data=init_vertex_data();
        vlist->vtx[i]=(ts_vertex *)calloc(1,sizeof(ts_vertex));
        vlist->vtx[i]->idx=i;
    }
    vlist->n=N;
    return vlist;
}

ts_vertex_data *init_vertex_data(){
    ts_vertex_data *data;
    data=(ts_vertex_data *)calloc(1,sizeof(ts_vertex_data));
    if(data==NULL)
        fatal("Fatal error reserving memory space for ts_vertex! Memory full?", 100);
    return data;
}


ts_bool vtx_add_neighbour(ts_vertex *vtx, ts_vertex *nvtx){
    ts_uint i;
@@ -46,12 +34,12 @@
    if(vtx==NULL || nvtx==NULL) return TS_FAIL;
    
    /*if it is already a neighbour don't add it to the list */
    for(i=0; i<vtx->data->neigh_no;i++){
        if(vtx->data->neigh[i]==nvtx) return TS_FAIL;
    for(i=0; i<vtx->neigh_no;i++){
        if(vtx->neigh[i]==nvtx) return TS_FAIL;
    }
    ts_uint nn=++vtx->data->neigh_no;
    vtx->data->neigh=(ts_vertex **)realloc(vtx->data->neigh, nn*sizeof(ts_vertex *));
    vtx->data->neigh[nn-1]=nvtx;
    ts_uint nn=++vtx->neigh_no;
    vtx->neigh=(ts_vertex **)realloc(vtx->neigh, nn*sizeof(ts_vertex *));
    vtx->neigh[nn-1]=nvtx;
/* This was a bug in creating DIPYRAMID (the neighbours were not in right
 * order).
 */
@@ -74,31 +62,31 @@
/* find a neighbour */
/* remove it from the list while shifting remaining neighbours up */
    ts_uint i,j=0;
    for(i=0;i<vtx->data->neigh_no;i++){
        if(vtx->data->neigh[i]!=nvtx){
            vtx->data->neigh[j]=vtx->data->neigh[i];
    for(i=0;i<vtx->neigh_no;i++){
        if(vtx->neigh[i]!=nvtx){
            vtx->neigh[j]=vtx->neigh[i];
            j++;
        }
    }
/* resize memory. potentionally time consuming */
    vtx->data->neigh_no--;
    vtx->data->neigh=(ts_vertex **)realloc(vtx->data->neigh,vtx->data->neigh_no*sizeof(ts_vertex *));
    if(vtx->data->neigh == NULL && vtx->data->neigh_no!=0)
    vtx->neigh_no--;
    vtx->neigh=(ts_vertex **)realloc(vtx->neigh,vtx->neigh_no*sizeof(ts_vertex *));
    if(vtx->neigh == NULL && vtx->neigh_no!=0)
        fatal("Reallocation of memory failed during removal of vertex neighbour in vtx_remove_neighbour",100);

/* repeat for the neighbour */
/* find a neighbour */
/* remove it from the list while shifting remaining neighbours up */
    for(i=0;i<nvtx->data->neigh_no;i++){
        if(nvtx->data->neigh[i]!=vtx){
            nvtx->data->neigh[j]=nvtx->data->neigh[i];
    for(i=0;i<nvtx->neigh_no;i++){
        if(nvtx->neigh[i]!=vtx){
            nvtx->neigh[j]=nvtx->neigh[i];
            j++;
        }
    }
/* resize memory. potentionally time consuming. */
    nvtx->data->neigh_no--;
    nvtx->data->neigh=(ts_vertex **)realloc(nvtx->data->neigh,nvtx->data->neigh_no*sizeof(ts_vertex *));
    if(nvtx->data->neigh == NULL && nvtx->data->neigh_no!=0)
    nvtx->neigh_no--;
    nvtx->neigh=(ts_vertex **)realloc(nvtx->neigh,nvtx->neigh_no*sizeof(ts_vertex *));
    if(nvtx->neigh == NULL && nvtx->neigh_no!=0)
        fatal("Reallocation of memory failed during removal of vertex neighbour in vtx_remove_neighbour",100);

    return TS_SUCCESS;
@@ -110,13 +98,13 @@
    ts_bond *bond;
    bond=bond_add(blist,vtx1,vtx2);
    if(bond==NULL) return TS_FAIL;
    vtx1->data->bond_no++;
    vtx1->bond_no++;
   // vtx2->data->bond_no++;

    vtx1->data->bond=(ts_bond **)realloc(vtx1->data->bond, vtx1->data->bond_no*sizeof(ts_bond *)); 
    vtx1->bond=(ts_bond **)realloc(vtx1->bond, vtx1->bond_no*sizeof(ts_bond *)); 
   // vtx2->data->bond=(ts_bond **)realloc(vtx2->data->bond, vtx2->data->bond_no*sizeof(ts_bond *)); 
    vtx1->data->bond[vtx1->data->bond_no-1]=bond;
   // vtx2->data->bond[vtx2->data->bond_no-1]=bond;
    vtx1->bond[vtx1->bond_no-1]=bond;
   // vtx2->ata->bond[vtx2->data->bond_no-1]=bond;
    return TS_SUCCESS;
}

@@ -139,20 +127,11 @@
}



ts_bool vtx_data_free(ts_vertex_data *data){
    if(data->neigh!=NULL)   free(data->neigh);
    if(data->tristar!=NULL) free(data->tristar);
    if(data->bond!=NULL)    free(data->bond);
//Cells are freed separately.
 //   if(data->cell!=NULL)    free(data->cell);
    free(data);
    return TS_SUCCESS;
}

/*not usable. can be deleted */
ts_bool vtx_free(ts_vertex  *vtx){
    vtx_data_free(vtx->data);
    if(vtx->neigh!=NULL)   free(vtx->neigh);
    if(vtx->tristar!=NULL) free(vtx->tristar);
    if(vtx->bond!=NULL)    free(vtx->bond);

    free(vtx);
    return TS_SUCCESS;
}
@@ -160,9 +139,9 @@
ts_bool vtx_list_free(ts_vertex_list *vlist){
    int i;
    for(i=0;i<vlist->n;i++){
        if(vlist->vtx[i]->data!=NULL) vtx_data_free(vlist->vtx[i]->data);
        if(vlist->vtx[i]!=NULL) vtx_free(vlist->vtx[i]);
    }
    free(*(vlist->vtx));
    //free(*(vlist->vtx));
    free(vlist->vtx);
    free(vlist);
    return TS_SUCCESS;
@@ -170,15 +149,14 @@

inline ts_double vtx_distance_sq(ts_vertex *vtx1, ts_vertex *vtx2){
    ts_double dist;
    ts_vertex_data *vd1=vtx1->data, *vd2=vtx2->data;
#ifdef TS_DOUBLE_DOUBLE
    dist=pow(vd1->x-vd2->x,2) + pow(vd1->y-vd2->y,2) + pow(vd1->z-vd2->z,2);
    dist=pow(vtx1->x-vtx2->x,2) + pow(vtx1->y-vtx2->y,2) + pow(vtx1->z-vtx2->z,2);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
    dist=powl(vd1->x-vd2->x,2) + powl(vd1->y-vd2->y,2) + powl(vd1->z-vd2->z,2);
    dist=powl(vtx1->x-vtx2->x,2) + powl(vtx1->y-vtx2->y,2) + powl(vtx1->z-vtx2->z,2);
#endif
#ifdef TS_DOUBLE_FLOAT
    dist=powf(vd1->x-vd2->x,2) + powf(vd1->y-vd2->y,2) + powf(vd1->z-vd2->z,2);
    dist=powf(vtx1->x-vtx2->x,2) + powf(vtx1->y-vtx2->y,2) + powf(vtx1->z-vtx2->z,2);
#endif
    return(dist);
}
@@ -189,32 +167,32 @@
    ts_double xk=vesicle->bending_rigidity;
    ts_uint i; 
    for(i=0;i<vesicle->vlist->n;i++){
        vesicle->vlist->vtx[i]->data->xk=xk;
        vesicle->vlist->vtx[i]->xk=xk;
    }
    return TS_SUCCESS;
}

inline ts_double vtx_direct(ts_vertex *vtx1, ts_vertex *vtx2, ts_vertex *vtx3){
    ts_double dX2=vtx2->data->x-vtx1->data->x;
    ts_double dY2=vtx2->data->y-vtx1->data->y;
    ts_double dZ2=vtx2->data->z-vtx1->data->z;
    ts_double dX3=vtx3->data->x-vtx1->data->x;
    ts_double dY3=vtx3->data->y-vtx1->data->y;
    ts_double dZ3=vtx3->data->z-vtx1->data->z;
    ts_double direct=vtx1->data->x*(dY2*dZ3 -dZ2*dY3)+ 
        vtx1->data->y*(dZ2*dX3-dX2*dZ3)+
        vtx1->data->z*(dX2*dY3-dY2*dX3);
    ts_double dX2=vtx2->x-vtx1->x;
    ts_double dY2=vtx2->y-vtx1->y;
    ts_double dZ2=vtx2->z-vtx1->z;
    ts_double dX3=vtx3->x-vtx1->x;
    ts_double dY3=vtx3->y-vtx1->y;
    ts_double dZ3=vtx3->z-vtx1->z;
    ts_double direct=vtx1->x*(dY2*dZ3 -dZ2*dY3)+ 
        vtx1->y*(dZ2*dX3-dX2*dZ3)+
        vtx1->z*(dX2*dY3-dY2*dX3);
    return(direct);    
}


inline ts_bool vertex_add_tristar(ts_vertex *vtx, ts_triangle *tristarmem){
    vtx->data->tristar_no++;
    vtx->data->tristar=(ts_triangle **)realloc(vtx->data->tristar,vtx->data->tristar_no*sizeof(ts_triangle *));
    if(vtx->data->tristar==NULL){
    vtx->tristar_no++;
    vtx->tristar=(ts_triangle **)realloc(vtx->tristar,vtx->tristar_no*sizeof(ts_triangle *));
    if(vtx->tristar==NULL){
            fatal("Reallocation of memory while adding tristar failed.",3);
    }
    vtx->data->tristar[vtx->data->tristar_no-1]=tristarmem;
    vtx->tristar[vtx->tristar_no-1]=tristarmem;
    return TS_SUCCESS;
}

@@ -225,22 +203,18 @@

ts_bool vtx_copy(ts_vertex *cvtx, ts_vertex *ovtx){
    memcpy((void *)cvtx,(void *)ovtx,sizeof(ts_vertex));
    cvtx->data=(ts_vertex_data *)malloc(sizeof(ts_vertex_data));
    memcpy((void *)cvtx->data,(void *)ovtx->data,sizeof(ts_vertex_data));
    cvtx->data->neigh=NULL;
    cvtx->data->neigh_no=0;
    cvtx->data->tristar_no=0;
    cvtx->data->bond_no=0;
    cvtx->data->tristar=NULL;
    cvtx->data->bond=NULL;
    cvtx->data->cell=NULL;
    cvtx->neigh=NULL;
    cvtx->neigh_no=0;
    cvtx->tristar_no=0;
    cvtx->bond_no=0;
    cvtx->tristar=NULL;
    cvtx->bond=NULL;
    cvtx->cell=NULL;
    return TS_SUCCESS;
}

ts_bool vtx_duplicate(ts_vertex *cvtx, ts_vertex *ovtx){
    memcpy((void *)cvtx,(void *)ovtx,sizeof(ts_vertex));
    cvtx->data=(ts_vertex_data *)malloc(sizeof(ts_vertex_data));
    memcpy((void *)cvtx->data,(void *)ovtx->data,sizeof(ts_vertex_data));
    return TS_SUCCESS;
}

@@ -256,12 +230,11 @@
    ts_vertex_list *vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
    vlist=memcpy((void *)vlist, (void *)ovlist, sizeof(ts_vertex_list));
    ts_vertex **vtx=(ts_vertex **)malloc(vlist->n*sizeof(ts_vertex *));
    ts_vertex *tlist=(ts_vertex *)calloc(vlist->n,sizeof(ts_vertex));
    vlist->vtx=vtx;
    if(vlist->vtx==NULL || tlist==NULL)
    if(vlist->vtx==NULL)
        fatal("Fatal error reserving memory space for vertex list! Could number of requsted vertices be too large?", 100);
    for(i=0;i<vlist->n;i++) {
        vlist->vtx[i]=&tlist[i];
        vlist->vtx[i]=(ts_vertex *)calloc(1,sizeof(ts_vertex));
        vlist->vtx[i]->idx=i;
        vtx_copy(vlist->vtx[i],ovlist->vtx[i]);
    }

 src/vertex.h

@@ -14,12 +14,10 @@
 *    @returns ts_bool value 1 on success, 0 otherwise
*/
ts_vertex_list *init_vertex_list(ts_uint N);
ts_vertex_data *init_vertex_data(void);
ts_bool vtx_add_neighbour(ts_vertex *vtx, ts_vertex *nvtx);
ts_bool vtx_add_cneighbour(ts_bond_list *blist,ts_vertex *vtx1,ts_vertex *vtx2);
ts_bool vtx_add_bond(ts_bond_list *blist,ts_vertex *vtx1,ts_vertex *vtx2);
ts_bool vtx_remove_neighbour(ts_vertex *vtx, ts_vertex *nvtx);
ts_bool vtx_data_free(ts_vertex_data *data);
ts_bool vtx_free(ts_vertex *vtx);
ts_bool vtx_list_free(ts_vertex_list *vlist);
inline ts_double vtx_distance_sq(ts_vertex *vtx1, ts_vertex *vtx2);

 src/vertexmove.c

@@ -16,8 +16,8 @@
*rn){
    ts_uint i;
    ts_double dist;
    ts_vertex *tvtx=(ts_vertex *)malloc(sizeof(ts_vertex));
    tvtx->data=init_vertex_data();
    ts_vertex *tvtx=(ts_vertex *)calloc(1,sizeof(ts_vertex));
//    tvtx->data=init_vertex_data();
    ts_bool retval; 
    ts_uint cellidx; 
    ts_double xold,yold,zold;
@@ -25,46 +25,50 @@
    ts_vertex *ovtx;

    //randomly we move the temporary vertex
    tvtx->data->x=vtx->data->x+vesicle->stepsize*(2.0*rn[0]-1.0);
    tvtx->data->y=vtx->data->y+vesicle->stepsize*(2.0*rn[1]-1.0);
    tvtx->data->z=vtx->data->z+vesicle->stepsize*(2.0*rn[2]-1.0);
    tvtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
    tvtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
    tvtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
    //check we if some length to neighbours are too much
    for(i=0;i<vtx->data->neigh_no;i++){
        dist=vtx_distance_sq(tvtx,vtx->data->neigh[i]);
        if(dist<1.0 || dist>vesicle->dmax) return TS_FAIL;
    for(i=0;i<vtx->neigh_no;i++){
        dist=vtx_distance_sq(tvtx,vtx->neigh[i]);
        if(dist<1.0 || dist>vesicle->dmax) {
        vtx_free(tvtx);
        return TS_FAIL;
        }
    }
fprintf(stderr,"Was here!\n");
//fprintf(stderr,"Was here!\n");
    //self avoidance check with distant vertices
     cellidx=vertex_self_avoidance(vesicle, tvtx);
    //check occupation number
     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx,tvtx);
    if(retval==TS_FAIL){
    vtx_free(tvtx);
        return TS_FAIL;
    } 
    
    //if all the tests are successful, then we update the vertex position
    xold=vtx->data->x;
    yold=vtx->data->y;
    zold=vtx->data->z;
    xold=vtx->x;
    yold=vtx->y;
    zold=vtx->z;
    ovtx=malloc(sizeof(ts_vertex));
    vtx_copy(ovtx,vtx);
    vtx->data->x=tvtx->data->x;
    vtx->data->y=tvtx->data->y;
    vtx->data->z=tvtx->data->z;
    vtx->x=tvtx->x;
    vtx->y=tvtx->y;
    vtx->z=tvtx->z;

    delta_energy=0;
    //update the normals of triangles that share bead i.
    for(i=0;i<vtx->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]);
    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
    //energy and curvature
    energy_vertex(vtx);
    delta_energy=vtx->data->xk*(vtx->data->energy - ovtx->data->energy);
    delta_energy=vtx->xk*(vtx->energy - ovtx->energy);
    //the same is done for neighbouring vertices
    for(i=0;i<vtx->data->neigh_no;i++){
        oenergy=vtx->data->neigh[i]->data->energy;
        energy_vertex(vtx->data->neigh[i]);
        delta_energy+=vtx->data->neigh[i]->data->xk*(vtx->data->neigh[i]->data->energy-oenergy);
    for(i=0;i<vtx->neigh_no;i++){
        oenergy=vtx->neigh[i]->energy;
        energy_vertex(vtx->neigh[i]);
        delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
    }
    fprintf(stderr, "DE=%f\n",delta_energy);
//    fprintf(stderr, "DE=%f\n",delta_energy);
    //MONTE CARLOOOOOOOO
    if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
@@ -78,18 +82,19 @@
#endif
    {
    //not accepted, reverting changes
    vtx->data->x=xold;
    vtx->data->y=yold;
    vtx->data->z=zold;
    vtx->x=xold;
    vtx->y=yold;
    vtx->z=zold;
    //update the normals of triangles that share bead i.
    for(i=0;i<vtx->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]);
    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
    //energy and curvature
    energy_vertex(vtx);
    //the same is done for neighbouring vertices
    for(i=0;i<vtx->data->neigh_no;i++) energy_vertex(vtx->data->neigh[i]);
  free(ovtx->data->bond_length);
    free(ovtx->data->bond_length_dual);
    for(i=0;i<vtx->neigh_no;i++) energy_vertex(vtx->neigh[i]);
    free(ovtx->bond_length);
    free(ovtx->bond_length_dual);
    free(ovtx);
    vtx_free(tvtx);
    return TS_FAIL; 
    }
}
@@ -97,9 +102,10 @@

    //TODO: change cell occupation if necessary!

    free(ovtx->data->bond_length);
    free(ovtx->data->bond_length_dual);
    free(ovtx->bond_length);
    free(ovtx->bond_length_dual);
    free(ovtx);
    vtx_free(tvtx);
    return TS_SUCCESS;
}


 src/vesicle.c

@@ -21,9 +21,9 @@
    ts_vertex **vtx=vesicle->vlist->vtx;
    ts_uint nn=vesicle->vlist->n;
    for(i=0;i<nn;i++){
        vtx[i]->data->x+=x;
        vtx[i]->data->y+=y;
        vtx[i]->data->z+=z;
        vtx[i]->x+=x;
        vtx[i]->y+=y;
        vtx[i]->z+=z;
    }
    return TS_SUCCESS;
}