From d06a9476717d788c426e88b6365d7a863928b5c9 Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@gmail.com> Date: Mon, 09 Dec 2013 00:29:28 +0000 Subject: [PATCH] Debugging bondflip --- src/bondflip.h | 2 src/bondflip.c | 105 ++++++++++++++++++++++------------------------------ src/bond.c | 6 ++- 3 files changed, 50 insertions(+), 63 deletions(-) diff --git a/src/bond.c b/src/bond.c index b70a871..7e98110 100644 --- a/src/bond.c +++ b/src/bond.c @@ -89,8 +89,10 @@ else i3=i; } bond->adjvtx[im]=bond->tria[1]->vertex[i3]; - - +#ifdef DEBUG + if(bond->adjvtx[0]==NULL) fatal("AdjVTX[0] is NULL",1); + if(bond->adjvtx[1]==NULL) fatal("AdjVTX[1] is NULL",1); +#endif return TS_SUCCESS; } diff --git a/src/bondflip.c b/src/bondflip.c index 06479e2..cc29021 100644 --- a/src/bondflip.c +++ b/src/bondflip.c @@ -27,63 +27,28 @@ */ ts_vertex *it=bond->vtx1; ts_vertex *k=bond->vtx2; - ts_uint nei,neip,neim; - ts_uint i; //j; + // ts_uint nei,neip,neim; + ts_uint i; ts_double oldenergy, delta_energy; - // ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lp2=NULL, *lm1=NULL, *lm2=NULL; - ts_vertex *kp=NULL,*km=NULL; - ts_triangle *lm=NULL, *lp=NULL; - - if(it->neigh->n< 3) return TS_FAIL; - if(k->neigh->n< 3) return TS_FAIL; - - if(k==NULL || it==NULL){ - fatal("In bondflip, number of neighbours of k or it is less than 3!",999); - } - - nei=0; - for(i=0;i<it->neigh->n;i++){ // Finds the nn of it, that is k - if(it->neigh->vtx[i]==k){ - nei=i; - break; - } - } - neip=nei+1; // I don't like it.. Smells like I must have it in correct order - neim=nei-1; - if(neip>=it->neigh->n) neip=0; - if((ts_int)neim<0) neim=it->neigh->n-1; /* casting is essential... If not -there the neim is never <0 !!! */ - // fprintf(stderr,"The numbers are: %u %u\n",neip, neim); - km=it->neigh->vtx[neim]; // We located km and kp - kp=it->neigh->vtx[neip]; - - //find lm and lp (triangles that is common with k and it) - for(i=0;i<it->tristar_no;i++){ - if(it->tristar[i]->vertex[0]==k || it->tristar[i]->vertex[1]==k || it->tristar[i]->vertex[2]==k){ - if(lm==NULL) lm=it->tristar[i]; - else lp=it->tristar[i]; - } - } - - if(lm==NULL || lp==NULL) { - fatal("In bondflip, cannot determine lp and lm",999); - } - - if(km==NULL || kp==NULL){ - fatal("In bondflip, cannot determine km and kp!",999); - } - - // fprintf(stderr,"I WAS HERE! after the 4 vertices are known!\n"); +#ifdef DEBUG + if(bond->adjvtx[0]==NULL || bond->adjvtx[1]==NULL) fatal ("single_bondflip_timestep: Adjvertex for called bond is null. This should not happen!",999); +#endif + ts_vertex *kp=bond->adjvtx[1],*km=bond->adjvtx[0]; +// ts_triangle *lm=NULL, *lp=NULL; /* test if the membrane is wrapped too much, so that kp is nearest neighbour of * km. If it is true, then don't flip! */ - for(i=0;i<km->neigh->n;i++){ +/*TODO: dec2013, i think this never happens! + for(i=0;i<km->neigh->n;i++){ if(km->neigh->vtx[i] == kp) return TS_FAIL; } +*/ + + // fprintf(stderr,"Membrane didn't wrap too much.. Continue.\n"); /* if bond would be too long, return... */ - if(vtx_distance_sq(km,kp) > vesicle->dmax ) return TS_FAIL; - // fprintf(stderr,"Bond will not be too long.. Continue.\n"); + if(vtx_distance_sq(bond->adjvtx[0],bond->adjvtx[1]) > vesicle->dmax ) return TS_FAIL; + // fprintf(stderr,"Bond will not be too long.. Continue.\n"); /* we make a bond flip. this is different than in original fortran */ // 0. step. Get memory prior the flip @@ -114,12 +79,13 @@ */ -// fprintf(stderr,"I WAS HERE! Before bondflip!\n"); - ts_flip_bond(k,it,km,kp, bond); + // fprintf(stderr,"I WAS HERE! Before bondflip!\n"); + ts_flip_bond(bond); // fprintf(stderr,"I WAS HERE! Bondflip successful!\n"); /* Calculating the new energy */ delta_energy=0; +/* TODO: why this neighs are recalculated? Not necessary! */ for(i=0;i<k->neigh->n;i++) energy_vertex(k->neigh->vtx[i]); for(i=0;i<kp->neigh->n;i++) energy_vertex(kp->neigh->vtx[i]); for(i=0;i<km->neigh->n;i++) energy_vertex(km->neigh->vtx[i]); @@ -150,7 +116,7 @@ // fprintf(stderr,"Failed to move, due to MC\n"); // ts_flip_bond(km,kp,it,k, bond); - ts_flip_bond(kp,km,k,it, bond); + ts_flip_bond(bond); /* @@ -184,16 +150,15 @@ } -ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp, -ts_bond *bond){ - +ts_bool ts_flip_bond(ts_bond *bond){ +fprintf(stderr,"Called!\n"); ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL; - ts_uint i,j; //lmidx, lpidx; -if(k==NULL || it==NULL || km==NULL || kp==NULL){ + ts_uint i; //lmidx, lpidx; +/*This is no longer necessary! if(k==NULL || it==NULL || km==NULL || kp==NULL){ fatal("ts_flip_bond: You called me with invalid pointers to vertices",999); -} +}*/ // 1. step. We find lm and lp from k->tristar ! - for(i=0;i<it->tristar_no;i++){ +/* for(i=0;i<it->tristar_no;i++){ for(j=0;j<k->tristar_no;j++){ if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1] @@ -229,6 +194,14 @@ } } } +*/ +//Upper block no longer neccessary. +ts_vertex *it=bond->vtx1; +ts_vertex *k=bond->vtx2; +ts_vertex *km=bond->adjvtx[0]; +ts_vertex *kp=bond->adjvtx[1]; +lm=bond->tria[0]; +lp=bond->tria[1]; /* // DEBUG TESTING! fprintf(stderr,"*** Naslov k=%d\n",k); @@ -244,8 +217,17 @@ // END DEBUG TESTING! */ +//find bonds between k and kp (adj[1] vtx) +//find bonds between it and km (adj[0] vtx) +for(i=0;i<it->tristar_no;i++){ + if((it->tristar[i]->vertex[0]==km || it->tristar[i]->vertex[1]==km || it->tristar[i]->vertex[2]==km) && (it->tristar[i]->vertex[0]==it || it->tristar[i]->vertex[1]==it || it->tristar[i]->vertex[2]==it)){ + lm2=it->tristar[i]; + } + else if ((it->tristar[i]->vertex[0]==kp || it->tristar[i]->vertex[1]==kp || it->tristar[i]->vertex[2]==kp) && (it->tristar[i]->vertex[0]==k || it->tristar[i]->vertex[1]==k || it->tristar[i]->vertex[2]==k)){ + lp1=it->tristar[i]; + } +} if(lm2==NULL || lp1==NULL) fatal("ts_flip_bond: Cannot find triangles lm2 and lp1!",999); - /* //DEBUG TESTING fprintf(stderr,"1. step: lm, lm2, lp1 and lp found!\n"); @@ -364,6 +346,9 @@ energy_vertex(km); energy_vertex(it); +//Extra steps for new bond data structure +bond->adjvtx[0]=k; +bond->adjvtx[1]=it; // END modifications to data structure! diff --git a/src/bondflip.h b/src/bondflip.h index 8d7ee37..f2f443d 100644 --- a/src/bondflip.h +++ b/src/bondflip.h @@ -3,5 +3,5 @@ ts_bool single_bondflip_timestep(ts_vesicle *vesicle, ts_bond *bond, ts_double *rn); -ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp, ts_bond *bond); +ts_bool ts_flip_bond(ts_bond *bond); #endif -- Gitblit v1.9.3