From 4e6285b2342038e3740ff9cea7c31364b4caa5e5 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Tue, 11 Feb 2014 12:03:29 +0000
Subject: [PATCH] Merge branch 'master' of bitbucket.org:samop/trisurf-ng
---
src/cell.c | 53 +++++++-
src/main.c | 25 +---
src/io.c | 33 +++-
src/tape | 9 +
src/cell.h | 6
src/general.h | 1
src/vertexmove.c | 109 ++++++++++-------
src/io.h | 3
src/timestep.c | 31 ++++
src/spherical_trisurf.c | 41 ++++++
src/timestep.h | 1
src/bond.c | 2
src/energy.c | 6
aclocal.m4 | 10
src/frame.c | 7 +
15 files changed, 226 insertions(+), 111 deletions(-)
diff --git a/aclocal.m4 b/aclocal.m4
index da06da5..3094178 100644
--- a/aclocal.m4
+++ b/aclocal.m4
@@ -1,4 +1,4 @@
-# generated automatically by aclocal 1.11.3 -*- Autoconf -*-
+# generated automatically by aclocal 1.11.6 -*- Autoconf -*-
# Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
# 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation,
@@ -14,8 +14,8 @@
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
-m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.68],,
-[m4_warning([this file was generated for autoconf 2.68.
+m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.69],,
+[m4_warning([this file was generated for autoconf 2.69.
You have another version of autoconf. It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically `autoreconf'.])])
@@ -38,7 +38,7 @@
[am__api_version='1.11'
dnl Some users find AM_AUTOMAKE_VERSION and mistake it for a way to
dnl require some minimum version. Point them to the right macro.
-m4_if([$1], [1.11.3], [],
+m4_if([$1], [1.11.6], [],
[AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl
])
@@ -54,7 +54,7 @@
# Call AM_AUTOMAKE_VERSION and AM_AUTOMAKE_VERSION so they can be traced.
# This function is AC_REQUIREd by AM_INIT_AUTOMAKE.
AC_DEFUN([AM_SET_CURRENT_AUTOMAKE_VERSION],
-[AM_AUTOMAKE_VERSION([1.11.3])dnl
+[AM_AUTOMAKE_VERSION([1.11.6])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])
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/io.c b/src/io.c
index f993e72..9e278ec 100644
--- a/src/io.c
+++ b/src/io.c
@@ -8,7 +8,7 @@
#include<stdlib.h>
#include <sys/types.h>
#include <dirent.h>
-
+#include "initial_distribution.h"
ts_bool print_vertex_list(ts_vertex_list *vlist){
ts_uint i;
@@ -294,12 +294,14 @@
-ts_bool parsetape(ts_vesicle *vesicle,ts_uint *iterations){
+ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations){
long int nshell=17,ncxmax=60, ncymax=60, nczmax=60; // THIS IS DUE TO CONFUSE BUG!
- char buf[255];
- long int brezveze=1;
+ char *buf=malloc(255*sizeof(char));
+ long int brezveze0=1;
+ long int brezveze1=1;
+ long int brezveze2=1;
ts_double xk0=25.0, dmax=1.67,stepsize=0.15;
- *iterations=1000;
+ long int iter=1000, init=1000, mcsw=1000;
cfg_opt_t opts[] = {
CFG_SIMPLE_INT("nshell", &nshell),
CFG_SIMPLE_FLOAT("dmax", &dmax),
@@ -308,12 +310,14 @@
CFG_SIMPLE_INT("nxmax", &ncxmax),
CFG_SIMPLE_INT("nymax", &ncymax),
CFG_SIMPLE_INT("nzmax", &nczmax),
- CFG_SIMPLE_INT("iterations",iterations),
+ CFG_SIMPLE_INT("iterations",&iter),
+ CFG_SIMPLE_INT("mcsweeps",&mcsw),
+ CFG_SIMPLE_INT("inititer", &init),
CFG_SIMPLE_BOOL("quiet",&quiet),
CFG_SIMPLE_STR("multiprocessing",buf),
- CFG_SIMPLE_INT("smp_cores",&brezveze),
- CFG_SIMPLE_INT("cluster_nodes",&brezveze),
- CFG_SIMPLE_INT("distributed_processes",&brezveze),
+ CFG_SIMPLE_INT("smp_cores",&brezveze0),
+ CFG_SIMPLE_INT("cluster_nodes",&brezveze1),
+ CFG_SIMPLE_INT("distributed_processes",&brezveze2),
CFG_END()
};
cfg_t *cfg;
@@ -326,6 +330,11 @@
else if(retval==CFG_PARSE_ERROR){
fatal("Invalid tape!",100);
}
+ ts_vesicle *vesicle;
+ *iterations=iter;
+ *inititer=init;
+ *mcsweeps=mcsw;
+ vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize);
vesicle->nshell=nshell;
vesicle->dmax=dmax*dmax;
vesicle->bending_rigidity=xk0;
@@ -334,9 +343,11 @@
vesicle->clist->ncmax[1]=ncymax;
vesicle->clist->ncmax[2]=nczmax;
vesicle->clist->max_occupancy=8;
+
cfg_free(cfg);
-// fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
- return TS_SUCCESS;
+ free(buf);
+ // fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
+ return vesicle;
}
diff --git a/src/io.h b/src/io.h
index 0dab974..20eab38 100644
--- a/src/io.h
+++ b/src/io.h
@@ -48,7 +48,6 @@
ts_bool write_vertex_vtk_file(ts_vesicle *vesicle,ts_char *filename, ts_char *text);
ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno);
ts_bool write_master_xml_file(ts_char *filename);
-ts_bool parsetape(ts_vesicle *vesicle,ts_uint *iterations);
-
+ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations);
#endif
diff --git a/src/main.c b/src/main.c
index d84b611..7e12c1c 100644
--- a/src/main.c
+++ b/src/main.c
@@ -18,14 +18,15 @@
*/
int main(int argv, char *argc[]){
+ts_uint inititer,mcsweeps, iterations;
+ts_vesicle *vesicle;
+/* THIS SHOULD GO INTO UNIT TEST
ts_bool retval;
-ts_uint i;
ts_vertex_list *vlist=init_vertex_list(5);
ts_vertex_list *vlist1;
ts_bond_list *blist=init_bond_list();
ts_triangle_list *tlist=init_triangle_list();
ts_cell_list *clist=init_cell_list(3,3,3,0.3);
-ts_vesicle *vesicle;
retval=vtx_add_cneighbour(blist,vlist->vtx[1],vlist->vtx[0]);
if(retval==TS_FAIL) printf("1. already a member or vertex is null!\n");
@@ -61,24 +62,10 @@
vtx_list_free(vlist1);
printf("Tests complete.\n");
-vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
-//parsetape(vesicle,&i);
+*/
+vesicle=parsetape(&mcsweeps, &inititer, &iterations);
+run_simulation(vesicle, mcsweeps, inititer, iterations);
-//these four must come from parsetype!
-vesicle->dmax=1.67*1.67;
-vesicle->stepsize=0.15;
-vesicle->clist->max_occupancy=8;
-vesicle->bending_rigidity=25.0;
-fprintf(stderr,"xk=%f",vesicle->bending_rigidity);
-
-centermass(vesicle);
-cell_occupation(vesicle);
-for(i=0;i<100;i++){
-single_timestep(vesicle);
-if(i%100==0){
-write_vertex_xml_file(vesicle,i/100);
-}
-}
write_master_xml_file("test.pvd");
write_dout_fcompat_file(vesicle,"dout");
vesicle_free(vesicle);
diff --git a/src/spherical_trisurf.c b/src/spherical_trisurf.c
index 82e500a..7fa40f9 100644
--- a/src/spherical_trisurf.c
+++ b/src/spherical_trisurf.c
@@ -19,11 +19,14 @@
*/
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);
//parsetape(vesicle,&i);
+
+//similar to nmax in fortran code
+ts_uint nmax;
//these four must come from parsetype!
vesicle->dmax=1.67*1.67;
@@ -33,6 +36,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);
@@ -42,14 +60,25 @@
calculateYlmi(vesicle);
calculateUlm(vesicle);
-
-
-for(i=0;i<10;i++){
+//preloop:
+for(i=0;i<1000;i++){
cell_occupation(vesicle);
for(j=0;j<1000;j++){
single_timestep(vesicle);
}
centermass(vesicle);
+ fprintf(stderr, "Preloop %d completed.\n",i+1);
+}
+
+nmax=1000;
+for(i=0;i<nmax;i++){
+ for(j=0;j<200;j++){
+ cell_occupation(vesicle);
+ for(k=0;k<5;k++){
+ single_timestep(vesicle);
+ }
+ centermass(vesicle);
+ }
vesicle_volume(vesicle);
r0=getR0(vesicle);
@@ -61,8 +90,10 @@
saveAvgUlm2(vesicle);
write_vertex_xml_file(vesicle,i);
- fprintf(stderr, "Loop %d completed.\n",i+1);
+ fprintf(stderr, "Loop %d out of %d completed.\n",i+1,nmax);
+
}
+
write_master_xml_file("test.pvd");
write_dout_fcompat_file(vesicle,"dout");
vesicle_free(vesicle);
diff --git a/src/tape b/src/tape
index 38518d3..0351c83 100644
--- a/src/tape
+++ b/src/tape
@@ -16,8 +16,13 @@
####### Program Control ############
-#how many iteration are there in a run?
-iterations=20000
+#how many MC sweeps between subsequent records of states to disk
+mcsweeps=100
+#how many initial mcsweeps*inititer MC sweeps before recording to disk?
+inititer=100
+#how many records do you want on the disk iteration are there in a run?
+iterations=1000
+
#shut up if we are using cluster!!!
quiet=false
diff --git a/src/timestep.c b/src/timestep.c
index 01cf777..9e5e61c 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -6,11 +6,31 @@
#include "timestep.h"
#include "vertexmove.h"
#include "bondflip.h"
+#include "frame.h"
+#include "io.h"
+ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations){
+ ts_uint i, j;
+
+ centermass(vesicle);
+ cell_occupation(vesicle);
+ ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
+ for(i=0;i<inititer+iterations;i++){
+ for(j=0;j<mcsweeps;j++){
+ single_timestep(vesicle);
+ }
+ centermass(vesicle);
+ cell_occupation(vesicle);
+ if(i>inititer){
+ write_vertex_xml_file(vesicle,i-inititer);
+ }
+ }
+ return TS_SUCCESS;
+}
ts_bool single_timestep(ts_vesicle *vesicle){
ts_bool retval;
ts_double rnvec[3];
- ts_uint i;
+ ts_uint i, b;
for(i=0;i<vesicle->vlist->n;i++){
rnvec[0]=drand48();
rnvec[1]=drand48();
@@ -19,13 +39,16 @@
}
// ts_int cnt=0;
- for(i=0;i<vesicle->blist->n;i++){
- rnvec[0]=drand48();
+ for(i=0;i<vesicle->vlist->n;i++){
+//why is rnvec needed in bondflip?
+/* rnvec[0]=drand48();
rnvec[1]=drand48();
rnvec[2]=drand48();
+*/
+ b=rand() % vesicle->blist->n;
//find a bond and return a pointer to a bond...
//call single_bondflip_timestep...
- retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[i],rnvec);
+ retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec);
// if(retval==TS_SUCCESS) cnt++;
}
// printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n);
diff --git a/src/timestep.h b/src/timestep.h
index c4015ae..e25234d 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -1,4 +1,5 @@
#ifndef _TIMESTEP_H
#define _TIMESTEP_H
ts_bool single_timestep(ts_vesicle *vesicle);
+ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations);
#endif
diff --git a/src/vertexmove.c b/src/vertexmove.c
index c62aeda..a79b2f0 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -11,6 +11,7 @@
//#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){
@@ -18,50 +19,69 @@
ts_double dist;
ts_bool retval;
ts_uint cellidx;
- ts_double xold,yold,zold;
ts_double delta_energy,oenergy;
- ts_vertex *ovtx;
- ts_vertex *tvtx=(ts_vertex *)calloc(1,sizeof(ts_vertex));
+ 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->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
- tvtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
- tvtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
- //check we if some length to neighbours are too much
+ //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(tvtx,vtx->neigh[i]);
+ dist=vtx_distance_sq(vtx,vtx->neigh[i]);
if(dist<1.0 || dist>vesicle->dmax) {
- vtx_free(tvtx);
-// fprintf(stderr,"Fail 1, dist=%f, vesicle->dmax=%f\n", dist, vesicle->dmax);
+ vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
return TS_FAIL;
}
}
//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_free(tvtx);
-// fprintf(stderr,"Fail 2\n");
+ 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->x;
- yold=vtx->y;
- zold=vtx->z;
- ovtx=malloc(sizeof(ts_vertex));
- vtx_copy(ovtx,vtx);
- vtx->x=tvtx->x;
- vtx->y=tvtx->y;
- vtx->z=tvtx->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->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
- //energy and curvature
+ oenergy=vtx->energy;
energy_vertex(vtx);
- delta_energy=vtx->xk*(vtx->energy - ovtx->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;
@@ -82,30 +102,27 @@
#endif
{
//not accepted, reverting changes
- vtx->x=xold;
- vtx->y=yold;
- vtx->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->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
- //energy and curvature
- energy_vertex(vtx);
- //the same is done for neighbouring vertices
- for(i=0;i<vtx->neigh_no;i++) energy_vertex(vtx->neigh[i]);
- free(ovtx->bond_length);
- free(ovtx->bond_length_dual);
- free(ovtx);
- vtx_free(tvtx);
+ 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!
-// fprintf(stderr,"Success!!\n");
- free(ovtx->bond_length);
- free(ovtx->bond_length_dual);
- free(ovtx);
- vtx_free(tvtx);
return TS_SUCCESS;
}
--
Gitblit v1.9.3