finished rewritting. start debugging

Samo Penic

2013-11-30 6fa0c93a432b95a6b888882494fcafb11c41515c

 src/initial_distribution.c

@@ -20,7 +20,7 @@
    ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize);

//TODO: debugging only. Please remove ASAP!
   vesicle->bending_rigidity=25.0*25.0;
   vesicle->bending_rigidity=25.0;

    vesicle->nshell=nshell;
    retval = vtx_set_global_values(vesicle);
@@ -172,16 +172,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]->neigh[0]->idx+1]); //always add 1st
          vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->neigh->vtx[0]->idx+1]); //always add 1st
          jjj=1;
          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]);
          for(l=2;l<=vtx[i]->neigh->n;l++){
              for(j=2;j<=vtx[i]->neigh->n;j++){
                  dist2=vtx_distance_sq(vtx[i]->neigh->vtx[j-1],vtx[i]->neigh->vtx[jj-1]);
                  direct=vtx_direct(vtx[i],vtx[i]->neigh->vtx[j-1],vtx[i]->neigh->vtx[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]->neigh[j-1]->idx+1]);
                    vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->neigh->vtx[j-1]->idx+1]);
                    jjj=jj;
                    jj=j;
                    break;
@@ -213,8 +213,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]->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
         for(k=0;k<vtx[i]->neigh->n;k++){ // has changed 0 to < instead of 1 and <=
            if(vtx[i]->neigh->vtx[k]==vtx[j]){  //if addresses matches it is the same
               bond_add(blist,vtx[i],vtx[j]);
               break;
            }
@@ -239,15 +239,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]->neigh_no;j++){
         for(jj=1;jj<=vtx[i]->neigh_no;jj++){
      for(j=1;j<=vtx[i]->neigh->n;j++){
         for(jj=1;jj<=vtx[i]->neigh->n;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]->neigh[j-1],vtx[i]->neigh[jj-1]);
            direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);				
            dist=vtx_distance_sq(vtx[i]->neigh->vtx[j-1],vtx[i]->neigh->vtx[jj-1]);
            direct=vtx_direct(vtx[i],vtx[i]->neigh->vtx[j-1],vtx[i]->neigh->vtx[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]);
            if(fabs(dist-A0*A0)<=eps && direct < 0.0 && vtx[i]->neigh->vtx[j-1]->idx+1 > i && vtx[i]->neigh->vtx[jj-1]->idx+1 >i){
               triangle_add(tlist,vtx[i],vtx[i]->neigh->vtx[j-1],vtx[i]->neigh->vtx[jj-1]);
            }   
         }   
      }
@@ -342,11 +342,11 @@
    ts_triangle **tria=tlist->tria -1;

    for(i=1;i<=vesicle->vlist->n;i++){
        for(j=1;j<=vtx[i]->neigh_no;j++){
            k1=vtx[i]->neigh[j-1];
        for(j=1;j<=vtx[i]->neigh->n;j++){
            k1=vtx[i]->neigh->vtx[j-1];
            jp=j+1;
            if(j == vtx[i]->neigh_no) jp=1;
            k2=vtx[i]->neigh[jp-1];
            if(j == vtx[i]->neigh->n) jp=1;
            k2=vtx[i]->neigh->vtx[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];