From 142a67fe82b830e5c7816914afa62445959c87ca Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Tue, 05 Nov 2013 14:04:21 +0000
Subject: [PATCH] changes in bondflip call. No need to bondflip all the bonds, but only as many bonds as there are vertices. Also, rnvec seems to be not needed for bondflip, so it is commented out
---
src/vertexmove.c | 123 ++++++++++++++++++++++++----------------
1 files changed, 73 insertions(+), 50 deletions(-)
diff --git a/src/vertexmove.c b/src/vertexmove.c
index 243463f..a79b2f0 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -11,60 +11,84 @@
//#include "io.h"
#include<stdio.h>
#include "vertexmove.h"
+#include <string.h>
ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double
*rn){
ts_uint i;
ts_double dist;
- ts_vertex *tvtx=(ts_vertex *)malloc(sizeof(ts_vertex));
- tvtx->data=init_vertex_data();
ts_bool retval;
ts_uint cellidx;
- ts_double xold,yold,zold;
ts_double delta_energy,oenergy;
- ts_vertex *ovtx;
+ ts_double costheta,sintheta,phi,r;
+ //This will hold all the information of vtx and its neighbours
+ ts_vertex backupvtx[20];
+ memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));
- //randomly we move the temporary vertex
- tvtx->data->x=vtx->data->x+vesicle->stepsize*(2.0*rn[0]-1.0);
- tvtx->data->y=vtx->data->y+vesicle->stepsize*(2.0*rn[1]-1.0);
- tvtx->data->z=vtx->data->z+vesicle->stepsize*(2.0*rn[2]-1.0);
- //check we if some length to neighbours are too much
- for(i=0;i<vtx->data->neigh_no;i++){
- dist=vtx_distance_sq(tvtx,vtx->data->neigh[i]);
- if(dist<1.0 || dist>vesicle->dmax) return TS_FAIL;
+ //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);
+
+ //random move in a sphere with radius stepsize:
+ r=vesicle->stepsize*rn[0];
+ phi=rn[1]*2*M_PI;
+ costheta=2*rn[2]-1;
+ sintheta=sqrt(1-pow(costheta,2));
+ vtx->x=vtx->x+r*sintheta*cos(phi);
+ vtx->y=vtx->y+r*sintheta*sin(phi);
+ vtx->z=vtx->z+r*costheta;
+
+
+ //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) {
+ vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
+ return TS_FAIL;
+ }
}
-fprintf(stderr,"Was here!\n");
//self avoidance check with distant vertices
- cellidx=vertex_self_avoidance(vesicle, tvtx);
+ cellidx=vertex_self_avoidance(vesicle, vtx);
//check occupation number
- retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx,tvtx);
+ 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;
}
-
- //if all the tests are successful, then we update the vertex position
- xold=vtx->data->x;
- yold=vtx->data->y;
- zold=vtx->data->z;
- ovtx=malloc(sizeof(ts_vertex));
- vtx_copy(ovtx,vtx);
- vtx->data->x=tvtx->data->x;
- vtx->data->y=tvtx->data->y;
- vtx->data->z=tvtx->data->z;
+
+
+ //if all the tests are successful, then energy for vtx and neighbours is calculated
+ for(i=0;i<vtx->neigh_no;i++){
+ 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->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]);
- //energy and curvature
+ for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
+ oenergy=vtx->energy;
energy_vertex(vtx);
- delta_energy=vtx->data->xk*(vtx->data->energy - ovtx->data->energy);
+ delta_energy=vtx->xk*(vtx->energy - oenergy);
//the same is done for neighbouring vertices
- for(i=0;i<vtx->data->neigh_no;i++){
- oenergy=vtx->data->neigh[i]->data->energy;
- energy_vertex(vtx->data->neigh[i]);
- delta_energy+=vtx->data->neigh[i]->data->xk*(vtx->data->neigh[i]->data->energy-oenergy);
+ for(i=0;i<vtx->neigh_no;i++){
+ oenergy=vtx->neigh[i]->energy;
+ energy_vertex(vtx->neigh[i]);
+ delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
}
- fprintf(stderr, "DE=%f\n",delta_energy);
+// fprintf(stderr, "DE=%f\n",delta_energy);
//MONTE CARLOOOOOOOO
if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
@@ -78,28 +102,27 @@
#endif
{
//not accepted, reverting changes
- vtx->data->x=xold;
- vtx->data->y=yold;
- vtx->data->z=zold;
+ 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));
+ }
+
//update the normals of triangles that share bead i.
- for(i=0;i<vtx->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]);
- //energy and curvature
- energy_vertex(vtx);
- //the same is done for neighbouring vertices
- for(i=0;i<vtx->data->neigh_no;i++) energy_vertex(vtx->data->neigh[i]);
- free(ovtx->data->bond_length);
- free(ovtx->data->bond_length_dual);
- free(ovtx);
+ for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
+
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
-
- //TODO: change cell occupation if necessary!
-
- free(ovtx->data->bond_length);
- free(ovtx->data->bond_length_dual);
- free(ovtx);
return TS_SUCCESS;
}
--
Gitblit v1.9.3