Fixed x-y bug for force

Miha

2016-07-15 c372c14d1f552a3219a733394141179c9e800683

 src/cell.c

@@ -1,87 +1,134 @@
/* vim: set ts=4 sts=4 sw=4 noet : */
#include<stdlib.h>
#include "general.h"
#include "vertex.h"
#include <stdio.h>
ts_bool init_cell_list(ts_cell_list *clist, ts_double stepsize){

ts_cell_list  *init_cell_list(ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
    ts_uint i;
    ts_uint nocells=clist->ncmax[0]*clist->ncmax[1]*clist->ncmax[2];
    clist->cell=malloc(nocells*sizeof(ts_cell));
    ts_uint nocells=ncmax1*ncmax2*ncmax3;
    ts_cell_list *clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
    if(clist==NULL) fatal("Error while allocating memory for cell list!",100);

    clist->ncmax[0]=ncmax1;
    clist->ncmax[1]=ncmax2;
    clist->ncmax[2]=ncmax3;
    clist->cellno=nocells;
    clist->dcell=1.0/(1.0 + stepsize);
    clist->shift=(ts_double) clist->ncmax[0]/2;
    clist->cellno=nocells;

    clist->cell=(ts_cell **)malloc(nocells*sizeof(ts_cell *));
    if(clist->cell==NULL) fatal("Error while allocating memory for cell list! ncmax too large?",101);

    for(i=0;i<nocells;i++){
        clist->cell[i].idx=i+1; // We enumerate cells! Probably never required!
        clist->cell[i].nvertex=0;
        clist->cell[i].vertex=NULL;
        clist->cell[i]=(ts_cell *)calloc(1,sizeof(ts_cell));
        if(clist->cell[i]==NULL) fatal("Error while allocating memory for cell list! ncmax too large?",102);
        clist->cell[i]->idx=i+1; // We enumerate cells! Probably never required!
    }
    return clist;
}

ts_bool cell_free(ts_cell* cell){
    if(cell->vertex!=NULL) free(cell->vertex);
    free(cell);
    return TS_SUCCESS;
}

ts_bool cell_list_free(ts_cell_list *clist){
    ts_uint i;
     ts_uint nocells=clist->ncmax[0]*clist->ncmax[1]*clist->ncmax[2];

    if(clist==NULL) return TS_FAIL;
    ts_uint nocells=clist->cellno;
    for(i=0;i<nocells;i++)
         if(clist->cell->vertex != NULL) free(clist->cell->vertex);
         if(clist->cell[i] != NULL) cell_free(clist->cell[i]);
    free(clist->cell);

    free(clist);
    return TS_SUCCESS;
}


inline ts_uint vertex_self_avoidance(ts_vesicle *vesicle, ts_vertex *vtx){
    ts_uint i,cellidx;
    ts_uint cellidx;
    ts_uint ncx, ncy,ncz;
    ts_cell_list *clist=vesicle->clist;
    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);

    ncx=(ts_uint)((vtx->x-vesicle->cm[0])*vesicle->clist.dcell+vesicle->clist.shift);
    ncy=(ts_uint)((vtx->y-vesicle->cm[1])*vesicle->clist.dcell+vesicle->clist.shift);
    ncz=(ts_uint)((vtx->z-vesicle->cm[2])*vesicle->clist.dcell+vesicle->clist.shift);
//    fprintf(stderr,"(ncx,ncy,ncz)=(%i %i %i)\t",ncx,ncy,ncz);
//    fprintf(stderr,"(ncxmax,ncymax,nczmax)=(%i %i %i)\n",vesicle->clist.ncmax[0], vesicle->clist.ncmax[1], vesicle->clist.ncmax[2]);
    if(ncx == vesicle->clist.ncmax[0]-1 || ncx == 2){
    if(ncx >= clist->ncmax[0]-1 || ncx <= 2){
        fatal("Vesicle is positioned outside the cell covered area. Coordinate x is the problem.",1500);
    }
    if(ncy == vesicle->clist.ncmax[1]-1 || ncy == 2){
    if(ncy >= clist->ncmax[1]-1 || ncy <= 2){
        fatal("Vesicle is positioned outside the cell covered area. Coordinate y is the problem.",1500);
    }
    if(ncz == vesicle->clist.ncmax[2]-1 || ncz == 2){
    if(ncz >= clist->ncmax[2]-1 || ncz <= 2){
        fatal("Vesicle is positioned outside the cell covered area. Coordinate z is the problem.",1500);
    }
    cellidx=ncz+(ncy-1)*vesicle->clist.ncmax[2] +
(ncx-1)*vesicle->clist.ncmax[2]* vesicle->clist.ncmax[1] - 1; // -1 is because of 0 based indexing
    return cellidx;
    cellidx=ncz+(ncy-1)*clist->ncmax[2] + (ncx-1)*clist->ncmax[2]* 
                                    clist->ncmax[1] - 1; // -1 is because of 0 based indexing
        return cellidx;
}



ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){
    
inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){
   ts_uint i;
   for(i=0;i<cell->nvertex;i++){
      if(cell->vertex[i]==vtx){
   //vertex is already in the cell!
         //fprintf(stderr,"VTX in the cell!\n");
         return TS_FAIL;
      }
   }
         //fprintf(stderr,"VTX added to the cell!\n");
    cell->nvertex++;
   cell->vertex=realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *));
      if(vtx->neigh == NULL){
         fatal("Reallocation of memory failed during insertion of vertex neighbour in vertex_add_neighbour",3);
   cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *));
      if(cell->vertex == NULL){
         fatal("Reallocation of memory failed during insertion of vertex in cell_add_vertex",3);
        }
    cell->vertex[cell->nvertex-1]=vtx;

   vtx->cell=cell;
    return TS_SUCCESS;
}

ts_bool cell_list_cell_ocupation_clear(ts_cell_list *clist){
inline ts_bool cell_remove_vertex(ts_cell *cell, ts_vertex *vtx){
   ts_uint i,j=0;
    for(i=0;i<cell->nvertex;i++){
        if(cell->vertex[i]!=vtx){
            cell->vertex[j]=cell->vertex[i];
            j++;
        }
    }
   if(j==i){
   fatal("Vertex was not in the cell!",3);
   } 
   //fprintf(stderr, "Vertex deleted from the cell!\n");

/* resize memory. potentionally time consuming */
    cell->nvertex--;
   cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *));
    if(vtx->neigh == NULL && vtx->neigh_no!=0)
      if(cell->vertex == NULL){
         fatal("Reallocation of memory failed during removal of vertex in cell_remove_vertex",3);
        }
   return TS_SUCCESS;
}

ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist){
    ts_uint i;
    for(i=0;i<clist->cellno;i++){
        if(clist->cell[i].vertex != NULL){
            free(clist->cell[i].vertex);
            clist->cell[i].vertex=NULL;
        if(clist->cell[i]->vertex != NULL){
            free(clist->cell[i]->vertex);
            clist->cell[i]->vertex=NULL;
        }
        clist->cell[i].nvertex=0;
        clist->cell[i]->nvertex=0;
    }
    return TS_SUCCESS;
}

ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx, ts_vertex *tvtx){

ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx){
    ts_uint ncx,ncy,ncz,remainder,cell_occupation;
    ts_uint i,j,k,l,neigh_cidx,mcell;
    ts_uint i,j,k,l,neigh_cidx;
    ts_double dist;
    ncx=(cellidx+1)/(clist->ncmax[2]*clist->ncmax[1])+1;
    ncx=(cellidx+1)/(clist->ncmax[2]*clist->ncmax[1])+1; //+1 because of zero indexing.
    remainder=(cellidx+1)%(clist->ncmax[2]*clist->ncmax[1]);
    ncy=remainder/clist->ncmax[2]+1;
    ncz=remainder%clist->ncmax[2];
@@ -92,26 +139,32 @@
            for(k=ncz-1;k<=ncz+1;k++){
                neigh_cidx=k+(j-1)*clist->ncmax[2]+(i-1)*clist->ncmax[2]*clist->ncmax[1] -1;
          //      fprintf(stderr,"neigh_cell_index=%i\n",neigh_cidx);
                cell_occupation=clist->cell[neigh_cidx].nvertex;
                cell_occupation=clist->cell[neigh_cidx]->nvertex;
          //      fprintf(stderr, "cell_occupation=%i\n",cell_occupation);
                if(cell_occupation>clist->max_occupancy){
                    fatal("Neighbouring cell occupation more than set max_occupancy value.",2000);
                }
// Now we check whether we didn't come close to some other vertices in the same
// cell!
                if(cell_occupation>1){
                if(cell_occupation>0){
                    for(l=0;l<cell_occupation;l++){
                        if(clist->cell[neigh_cidx].vertex[l]!=vtx){

            //carefull with this checks!
                        if(clist->cell[neigh_cidx]->vertex[l]!=vtx){
                    //        fprintf(stderr,"calling dist on vertex %i\n",l);
                           dist=vertex_distance_sq(clist->cell[neigh_cidx].vertex[l],tvtx);
                    //        fprintf(stderr,"dist was %f\n",dist);
                            if(dist<1) return TS_FAIL;
                           dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],vtx);

//            if(vtx->idx==1)
//            fprintf(stderr,"VTX(0) ima bliznji vertex z indeksom, %d, tipa %d \n", clist->cell[neigh_cidx]->vertex[l]->idx, clist->cell[neigh_cidx]->vertex[l]->id);
//            if(vtx->idx==0 && clist->cell[neigh_cidx]->vertex[l]->idx==0)
//                            fprintf(stderr,"*** dist was %f\n",dist);
				
                            if(dist<=1.0 || (dist<=clist->dmin_interspecies && (clist->cell[neigh_cidx]->vertex[l]->id != vtx->id))) return TS_FAIL;
                        }
                    }
                }
            }
        }
    }
    
    } 
    return TS_SUCCESS;
}