From c0ae90d68622063a44af443acc01f7e4925fbc56 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Thu, 18 Sep 2014 10:17:22 +0000
Subject: [PATCH] Added constant area. Unfortunately dump has changed sue to change in datastructure
---
src/timestep.c | 11 ++-
src/io.c | 1
src/vesicle.h | 2
src/tape | 4 +
src/vesicle.c | 17 +++++
src/bondflip.c | 47 +++++++++++++--
src/general.h | 6 +
src/vertexmove.c | 39 ++++++++++--
8 files changed, 106 insertions(+), 21 deletions(-)
diff --git a/src/bondflip.c b/src/bondflip.c
index 4c8b50c..bb6bd33 100644
--- a/src/bondflip.c
+++ b/src/bondflip.c
@@ -31,7 +31,7 @@
ts_vertex *k=bond->vtx2;
ts_uint nei,neip,neim;
ts_uint i,j;
- ts_double oldenergy, delta_energy, dvol=0.0;
+ ts_double oldenergy, delta_energy, dvol=0.0, darea=0.0;
ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL;
ts_vertex *kp,*km;
@@ -170,7 +170,7 @@
//Neigbours of k, it, km, kp don't change its energy.
if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){dvol = -lm->volume - lp->volume;}
-
+ if(vesicle->tape->constareaswitch==2){darea=-lm->area-lp->area;}
/* vesicle_volume(vesicle);
fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
*/
@@ -192,6 +192,39 @@
dvol = dvol + lm->volume + lp->volume;
if(vesicle->pswitch==1) delta_energy-= vesicle->pressure*dvol;
}
+ if(vesicle->tape->constareaswitch==2){
+ darea=darea+lm->area+lp->area;
+/*check whether the dvol is gt than epsvol */
+ if(fabs(vesicle->area+darea-A0)>epsarea){
+ //restore old state.
+ /* restoration procedure copied from few lines below */
+ for(i=0;i<4;i++){
+ // fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ free(orig_vtx[i]->neigh);
+ free(orig_vtx[i]->tristar);
+ free(orig_vtx[i]->bond);
+ free(orig_tria[i]->neigh);
+ memcpy((void *)orig_vtx[i],(void *)bck_vtx[i],sizeof(ts_vertex));
+ memcpy((void *)orig_tria[i],(void *)bck_tria[i],sizeof(ts_triangle));
+ // fprintf(stderr,"Restored vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ /* level 2 pointers are redirected*/
+ }
+ memcpy(bond,bck_bond,sizeof(ts_bond));
+ for(i=0;i<4;i++){
+ free(bck_vtx[i]);
+ free(bck_tria[i]);
+ /* fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d =",i, orig_vtx[i]->neigh_no );
+ for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx);
+ fprintf(stderr,"\n"); */
+ }
+ free(bck_bond);
+ return TS_FAIL;
+
+ }
+ }
+
+
+
if(vesicle->tape->constvolswitch == 2){
/*check whether the dvol is gt than epsvol */
@@ -305,12 +338,14 @@
/* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
// fprintf(stderr,"SUCCESS!!!\n");
- if(vesicle->tape->constvolswitch == 2){
- vesicle->volume+=dvol;
- } else
- if(vesicle->tape->constvolswitch == 1){
+ if(vesicle->tape->constvolswitch == 2){
+ vesicle->volume+=dvol;
+ } else if(vesicle->tape->constvolswitch == 1){
constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
}
+ if(vesicle->tape->constareaswitch==2){
+ vesicle->area+=darea;
+ }
// delete all backups
for(i=0;i<4;i++){
free(bck_vtx[i]->neigh);
diff --git a/src/general.h b/src/general.h
index c768c0f..393eafe 100644
--- a/src/general.h
+++ b/src/general.h
@@ -261,6 +261,7 @@
long int R_nucleus;
long int pswitch;
long int constvolswitch;
+ long int constareaswitch;
ts_double constvolprecision;
char *multiprocessing;
long int brezveze0;
@@ -305,7 +306,7 @@
ts_int pswitch;
ts_tape *tape;
ts_double R_nucleus;
-
+ ts_double area;
} ts_vesicle;
@@ -314,8 +315,9 @@
int quiet;
ts_double V0;
+ts_double A0;
ts_double epsvol;
-
+ts_double epsarea;
/* FUNCTIONS */
/** Non-fatal error function handler:
diff --git a/src/io.c b/src/io.c
index f9b0164..acc61c3 100644
--- a/src/io.c
+++ b/src/io.c
@@ -1008,6 +1008,7 @@
CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
CFG_SIMPLE_INT("pswitch",&tape->pswitch),
CFG_SIMPLE_INT("constvolswitch",&tape->constvolswitch),
+ CFG_SIMPLE_INT("constareaswitch",&tape->constareaswitch),
CFG_SIMPLE_FLOAT("constvolprecision",&tape->constvolprecision),
CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
diff --git a/src/tape b/src/tape
index bb4b1db..9c7c574 100644
--- a/src/tape
+++ b/src/tape
@@ -20,10 +20,12 @@
constvolswitch=2
constvolprecision=1e-14
+#Constant area constraint (0 disable constant area, 2 enable constant area with epsarea)
+constareaswitch=2
####### Polymer (brush) definitions ###########
# npoly is a number of polymers attached to npoly distinct vertices on vesicle
-npoly=5
+npoly=0
# nmono is a number of monomers in each polymer
nmono=10
# Spring constant between monomers of the polymer
diff --git a/src/timestep.c b/src/timestep.c
index fcb0ff1..315c1d1 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -19,7 +19,7 @@
ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){
ts_uint i, j,k,l,m;
ts_double r0,kc1,kc2,kc3,kc4;
- ts_double l1,l2,l3,volume=0.0,area=0.0,vmsr,bfsr, vmsrt, bfsrt;
+ ts_double l1,l2,l3,vmsr,bfsr, vmsrt, bfsrt;
ts_ulong epochtime;
FILE *fd1,*fd2=NULL;
char filename[10000];
@@ -49,8 +49,11 @@
centermass(vesicle);
cell_occupation(vesicle);
vesicle_volume(vesicle); //needed for constant volume at this moment
+ vesicle_area(vesicle); //needed for constant area at this moment
V0=vesicle->volume;
+ A0=vesicle->area;
epsvol=4.0*sqrt(2.0*M_PI)/pow(3.0,3.0/4.0)*V0/pow(vesicle->tlist->n,3.0/2.0);
+ epsarea=A0/(ts_double)vesicle->tlist->n;
// fprintf(stderr, "DVol=%1.16f (%1.16f), V0=%1.16f\n", epsvol,0.003e-2*V0,V0);
if(start_iteration<inititer) ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
for(i=start_iteration;i<inititer+iterations;i++){
@@ -77,8 +80,8 @@
write_master_xml_file(command_line_args.output_fullfilename);
epochtime=get_epoch();
gyration_eigen(vesicle, &l1, &l2, &l3);
- vesicle_volume(vesicle); //calculates just volume. Area is not added to ts_vesicle yet!
- get_area_volume(vesicle, &area,&volume); //that's why I must recalculate area (and volume for no particular reason).
+ vesicle_volume(vesicle); //calculates just volume.
+ vesicle_area(vesicle); //calculates area.
r0=getR0(vesicle);
if(vesicle->sphHarmonics!=NULL){
preparationSh(vesicle,r0);
@@ -117,7 +120,7 @@
}
- fprintf(fd, "%lu %u %e %e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e\n",epochtime,i,vmsr,bfsr,volume, area,l1,l2,l3,kc1, kc2, kc3,kc4);
+ fprintf(fd, "%lu %u %e %e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e\n",epochtime,i,vmsr,bfsr,vesicle->volume, vesicle->area,l1,l2,l3,kc1, kc2, kc3,kc4);
fflush(fd);
// sprintf(filename,"timestep-%05d.pov",i-inititer);
diff --git a/src/vertexmove.c b/src/vertexmove.c
index 28a8880..c552a87 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -20,7 +20,7 @@
ts_double dist;
ts_bool retval;
ts_uint cellidx;
- ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0;
+ ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0;
ts_double costheta,sintheta,phi,r;
//This will hold all the information of vtx and its neighbours
ts_vertex backupvtx[20], *constvol_vtx_moved=NULL, *constvol_vtx_backup=NULL;
@@ -41,7 +41,7 @@
// 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:
+//random move in a sphere with radius stepsize:
r=vesicle->stepsize*rn[0];
phi=rn[1]*2*M_PI;
costheta=2*rn[2]-1;
@@ -51,12 +51,12 @@
vtx->z=vtx->z+r*costheta;
- //distance with neighbours check
+//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;
+ vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
+ return TS_FAIL;
}
}
@@ -87,14 +87,19 @@
}
- //if all the tests are successful, then energy for vtx and neighbours is calculated
+//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));
}
if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){
for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
- };
+ }
+
+ if(vesicle->tape->constareaswitch==2){
+ for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area;
+
+ }
delta_energy=0;
@@ -118,6 +123,22 @@
if(vesicle->pswitch==1) delta_energy-=vesicle->pressure*dvol;
};
+ if(vesicle->tape->constareaswitch==2){
+ /* check whether the darea is gt epsarea */
+ for(i=0;i<vtx->tristar_no;i++) darea+=vtx->tristar[i]->area;
+ if(fabs(vesicle->area+darea-A0)>epsarea){
+ //restore old state.
+ 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));
+ }
+ for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
+ //fprintf(stderr,"fajlam!\n");
+ return TS_FAIL;
+ }
+
+
+ }
if(vesicle->tape->constvolswitch==2){
/*check whether the dvol is gt than epsvol */
@@ -211,6 +232,10 @@
if(vesicle->tape->constvolswitch == 1){
constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
}
+
+ if(vesicle->tape->constareaswitch==2){
+ vesicle->area+=darea;
+ }
// if(oldcellidx);
//END MONTE CARLOOOOOOO
// vesicle_volume(vesicle);
diff --git a/src/vesicle.c b/src/vesicle.c
index f63eb12..5fc51cc 100644
--- a/src/vesicle.c
+++ b/src/vesicle.c
@@ -60,3 +60,20 @@
vesicle->volume=volume;
return TS_SUCCESS;
}
+
+/* @brief Function makes a sum of partial areas of each triangle.
+ *
+ *
+ *
+ */
+ts_bool vesicle_area(ts_vesicle *vesicle){
+ ts_double area;
+ ts_uint i;
+ ts_triangle **tria=vesicle->tlist->tria;
+ area=0;
+ for(i=0;i<vesicle->tlist->n;i++){
+ area=area+tria[i]->area;
+ }
+ vesicle->area=area;
+ return TS_SUCCESS;
+}
diff --git a/src/vesicle.h b/src/vesicle.h
index 2992491..50a9e3e 100644
--- a/src/vesicle.h
+++ b/src/vesicle.h
@@ -5,5 +5,5 @@
ts_bool vesicle_translate(ts_vesicle *vesicle,ts_double x, ts_double y, ts_double z);
ts_bool vesicle_free(ts_vesicle *vesicle);
ts_bool vesicle_volume(ts_vesicle *vesicle);
-
+ts_bool vesicle_area(ts_vesicle *vesicle);
#endif
--
Gitblit v1.9.3