From a63f1719d2c7fd2c69accc0eb3eb038af50e555e Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Fri, 13 Jul 2012 13:10:01 +0000
Subject: [PATCH] Fixed frame to make changes in cell.
---
src/cell.c | 53 ++++++++++++++---
src/cell.h | 6 +-
src/spherical_trisurf.c | 28 +++++++-
src/bond.c | 2
src/energy.c | 6 +
src/general.h | 1
src/frame.c | 7 ++
src/vertexmove.c | 33 +++++++++-
8 files changed, 109 insertions(+), 27 deletions(-)
diff --git a/src/bond.c b/src/bond.c
index 9d82ea2..af39532 100644
--- a/src/bond.c
+++ b/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];
diff --git a/src/cell.c b/src/cell.c
index 47d7fab..246866b 100644
--- a/src/cell.c
+++ b/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){
+
+ //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],tvtx);
+ dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],vtx);
// fprintf(stderr,"dist was %f\n",dist);
- if(dist<1) return TS_FAIL;
+ if(dist<=1.0) return TS_FAIL;
}
}
}
}
}
- }
-
+ }
return TS_SUCCESS;
}
diff --git a/src/cell.h b/src/cell.h
index 76ba105..20a58e4 100644
--- a/src/cell.h
+++ b/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
diff --git a/src/energy.c b/src/energy.c
index 4d93753..abbb8e7 100644
--- a/src/energy.c
+++ b/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);
diff --git a/src/frame.c b/src/frame.c
index b01c2d2..5cfede9 100644
--- a/src/frame.c
+++ b/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]);
diff --git a/src/general.h b/src/general.h
index 6263862..f3328f6 100644
--- a/src/general.h
+++ b/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;
diff --git a/src/spherical_trisurf.c b/src/spherical_trisurf.c
index e3e1986..950ecf8 100644
--- a/src/spherical_trisurf.c
+++ b/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++){
- cell_occupation(vesicle);
- for(j=0;j<1000;j++){
+ for(j=0;j<20;j++){
+ cell_occupation(vesicle);
+ for(k=0;k<5;k++){
single_timestep(vesicle);
+ }
+ centermass(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);
diff --git a/src/vertexmove.c b/src/vertexmove.c
index b98835c..f38e3c7 100644
--- a/src/vertexmove.c
+++ b/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;
@@ -77,9 +92,8 @@
//not accepted, reverting changes
vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
for(i=0;i<vtx->neigh_no;i++){
- vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
+ 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;
}
--
Gitblit v1.9.3