Fixed frame to make changes in cell.

Samo Penic

2012-07-13 a63f1719d2c7fd2c69accc0eb3eb038af50e555e

Fixed frame to make changes in cell.

 src/bond.c

@@ -28,7 +28,7 @@
    //NOW insert vertices into data!    
    blist->bond[blist->n - 1]->vtx1=vtx1;    
    blist->bond[blist->n - 1]->vtx2=vtx2;

    blist->bond[blist->n - 1]->tainted=0;
    //Should we calculate bond length NOW?
    
    return blist->bond[blist->n-1];

 src/cell.c

@@ -53,13 +53,13 @@
    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){
    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 == 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 == 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)*clist->ncmax[2] + (ncx-1)*clist->ncmax[2]* 
@@ -69,16 +69,48 @@


//TODO: looks ok, but debug anyway in the future
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=(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;
}

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;
@@ -93,7 +125,7 @@
}

// TODO: compiles ok, but it is completely untested and undebugged. It was debugged before rewrite, but this was long time ago.
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;
    ts_double dist;
@@ -117,17 +149,18 @@
// cell!
                if(cell_occupation>1){
                    for(l=0;l<cell_occupation;l++){
                        if(clist->cell[neigh_cidx]->vertex[l]!=vtx){
                    //        fprintf(stderr,"calling dist on vertex %i\n",l);
                           dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],tvtx);
                    //        fprintf(stderr,"dist was %f\n",dist);
                            if(dist<1) return TS_FAIL;
                        }
                    }
                }
            }
        }
    }
    
                //carefull with this checks!
                        if(clist->cell[neigh_cidx]->vertex[l]->idx!=vtx->idx){
                    //        fprintf(stderr,"calling dist on vertex %i\n",l);
                           dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],vtx);
                    //        fprintf(stderr,"dist was %f\n",dist);
                            if(dist<=1.0) return TS_FAIL;
                        }
                    }
                }
            }
        }
    } 
    return TS_SUCCESS;
}

 src/cell.h

@@ -4,8 +4,8 @@
ts_bool cell_free(ts_cell* cell);
ts_bool cell_list_free(ts_cell_list *clist);
inline ts_uint vertex_self_avoidance(ts_vesicle *vesicle, ts_vertex *vtx);
ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx);
inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx);
inline ts_bool cell_remove_vertex(ts_cell *cell, ts_vertex *vtx);
ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist);

ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx, ts_vertex *tvtx);
inline ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx);
#endif

 src/energy.c

@@ -28,7 +28,7 @@
    ts_uint jjp,jjm;
    ts_vertex *j,*jp, *jm;
    ts_triangle *jt;
    ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0;
    ts_double s=0.0,xh=0.0,yh=0.0,zh=0.0,txn=0.0,tyn=0.0,tzn=0.0;
    ts_double x1,x2,x3,ctp,ctm,tot,xlen;
    ts_double h,ht;
    for(jj=1; jj<=vtx->neigh_no;jj++){
@@ -72,7 +72,9 @@
#endif
        tot=ctp+ctm;
        tot=0.5*tot;

        xlen=vtx_distance_sq(j,vtx);
/*
#ifdef  TS_DOUBLE_DOUBLE 
        vtx->bond[jj-1]->bond_length=sqrt(xlen); 
#endif
@@ -84,7 +86,7 @@
#endif

        vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length;

*/
        s+=tot*xlen;
        xh+=tot*(j->x - vtx->x);
        yh+=tot*(j->y - vtx->y);

 src/frame.c

@@ -23,6 +23,10 @@
        vtx[i]->z-=vesicle->cm[2]; 
    } 

    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;

    return TS_SUCCESS;
}

@@ -37,7 +41,8 @@
    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]->cell=vesicle->clist->cell[cellidx];
//    already done in cell_add_vertex
//    vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx];

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


 src/general.h

@@ -169,6 +169,7 @@
    ts_vertex *vtx2;
    ts_double bond_length;
    ts_double bond_length_dual;
    ts_bool tainted;
};
typedef struct ts_bond ts_bond;


 src/spherical_trisurf.c

@@ -19,7 +19,7 @@
*/
ts_bool saveAvgUlm2(ts_vesicle *vesicle);
int main(int argv, char *argc[]){
ts_uint i,j;
ts_uint i,j,k;
ts_vesicle *vesicle;
ts_double r0;
vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
@@ -33,6 +33,21 @@
//fprintf(stderr,"xk=%f",vesicle->bending_rigidity);

    centermass(vesicle);
cell_occupation(vesicle);

//test if the structure is internally organized into cells correctly 
ts_uint cind;
for(i=0;i<vesicle->vlist->n;i++){
    cind=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]);

    if(vesicle->clist->cell[cind]==vesicle->vlist->vtx[i]->cell){
        //fprintf(stdout,"(T) Idx match!\n");
    } else {
        fprintf(stderr,"(T) ***** Idx don't match!\n");

    }
}
//end test
vesicle->sphHarmonics=sph_init(vesicle->vlist, 21);

vesicle_volume(vesicle);
@@ -43,13 +58,14 @@
calculateUlm(vesicle);



for(i=0;i<1;i++){
    for(j=0;j<20;j++){
    cell_occupation(vesicle);
    for(j=0;j<1000;j++){
        for(k=0;k<5;k++){
        single_timestep(vesicle);
    }    
    centermass(vesicle);
    }	
    vesicle_volume(vesicle);
    r0=getR0(vesicle);

@@ -62,7 +78,9 @@

    write_vertex_xml_file(vesicle,i);
    fprintf(stderr, "Loop %d completed.\n",i+1);

}

write_master_xml_file("test.pvd");
write_dout_fcompat_file(vesicle,"dout");
vesicle_free(vesicle);

 src/vertexmove.c

@@ -23,11 +23,22 @@
    //This will hold all the information of vtx and its neighbours
    ts_vertex backupvtx[20];
    memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));

    //Some stupid tests for debugging cell occupation!
/*         cellidx=vertex_self_avoidance(vesicle, vtx);
    if(vesicle->clist->cell[cellidx]==vtx->cell){
        fprintf(stderr,"Idx match!\n");
    } else {
        fprintf(stderr,"***** Idx don't match!\n");
        fatal("ENding.",1);
    }
*/

        //temporarly moving the vertex
    vtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
        vtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
        vtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
        //check we if some length to neighbours are too much
        //distance with neighbours check
    for(i=0;i<vtx->neigh_no;i++){
        dist=vtx_distance_sq(vtx,vtx->neigh[i]);
        if(dist<1.0 || dist>vesicle->dmax) {
@@ -38,7 +49,8 @@
    //self avoidance check with distant vertices
     cellidx=vertex_self_avoidance(vesicle, vtx);
    //check occupation number
     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,&backupvtx[0],vtx);
     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	
    if(retval==TS_FAIL){
        vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
        return TS_FAIL;
@@ -50,11 +62,14 @@
    memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
    }



    delta_energy=0;
    //update the normals of triangles that share bead i.
    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
    oenergy=vtx->energy;
    energy_vertex(vtx);
    delta_energy=vtx->xk*(vtx->energy - (&backupvtx[0])->energy);
    delta_energy=vtx->xk*(vtx->energy - oenergy);
    //the same is done for neighbouring vertices
    for(i=0;i<vtx->neigh_no;i++){
        oenergy=vtx->neigh[i]->energy;
@@ -79,7 +94,6 @@
    for(i=0;i<vtx->neigh_no;i++){
    vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
    }
//    fprintf(stderr,"Reverted\n");
    
    //update the normals of triangles that share bead i.
   for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
@@ -87,6 +101,15 @@
    return TS_FAIL; 
    }
}
		
//    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
    if(vtx->cell!=vesicle->clist->cell[cellidx]){
        retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
//        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
        if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
		
    }
//    if(oldcellidx);
    //END MONTE CARLOOOOOOO
    return TS_SUCCESS;
}