From bd826de2f539f2e48c8c01d2d7f9f34c7e97104a Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Fri, 13 May 2016 07:43:27 +0000
Subject: [PATCH] Fix in trisurf output, inhibiting print of successful reconstruction. Multiple fixes and improvements in python module. Added symlinking of tapes into the running directories and dumping tapes from snapshots into tape files.
---
src/constvol.c | 135 ++++++++++++++++++++++++++++----------------
1 files changed, 85 insertions(+), 50 deletions(-)
diff --git a/src/constvol.c b/src/constvol.c
index da07add..63c07e0 100644
--- a/src/constvol.c
+++ b/src/constvol.c
@@ -1,3 +1,4 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
#include<stdlib.h>
#include<stdio.h>
#include<string.h>
@@ -9,69 +10,84 @@
#include "vertex.h"
#include "cell.h"
-ts_bool constvolume(ts_vesicle *vesicle, ts_vertex *vtx_avoid, ts_double Vol, ts_double *retEnergy, ts_vertex *vtx_moved, ts_vertex *vtx_backup){
-
+ts_bool constvolume(ts_vesicle *vesicle, ts_vertex *vtx_avoid, ts_double Vol, ts_double *retEnergy, ts_vertex **vtx_moved_retval, ts_vertex **vtx_backup){
+ ts_vertex *vtx_moved;
ts_uint vtxind,i,j;
- ts_uint Ntries=20;
+ ts_uint Ntries=100;
ts_vertex *backupvtx;
- ts_double Rv, dh, dvol, voldiff, oenergy,delta_energy;
-
+ ts_double Rv, dh, dvol, volFirst, voldiff, oenergy,delta_energy;
backupvtx=(ts_vertex *)calloc(sizeof(ts_vertex),10);
-
+ ts_double l0 = (1.0 + sqrt(vesicle->dmax))/2.0; //make this a global constant if necessary
for(i=0;i<Ntries;i++){
vtxind=rand() % vesicle->vlist->n;
vtx_moved=vesicle->vlist->vtx[vtxind];
- if(vtx_moved==vtx_avoid) continue;
-
+ /* chosen vertex must not be a nearest neighbour. WASTODO: probably must.
+ * DONE: solved below, when using different algorithm.
+ * extend search in case of bondflip */
+/* if(vtx_moved==vtx_avoid) continue;
for(j=0;j<vtx_moved->neigh_no;j++){
if(vtx_moved->neigh[j]==vtx_avoid) continue;
}
-
- memcpy((void *)&backupvtx[0],(void *)vtx_moved,sizeof(ts_vertex));
- //move vertex in specified direction. first try, test move!
+*/
+/* different check of nearest neighbour (and second nearest neighbour) check.
+ * Checking the distance between vertex and vertex to be moved to assure
+ * constant volume. Solves upper todo problem. */
- Rv=sqrt(pow(vtx_moved->x,2)+pow(vtx_moved->y,2)+pow(vtx_moved->z,2));
- dh=2*Rv*vesicle->dmax/sqrt(3);
- vtx_moved->x=vtx_moved->x*(1-dh/Rv);
- vtx_moved->y=vtx_moved->y*(1-dh/Rv);
- vtx_moved->z=vtx_moved->z*(1-dh/Rv);
-
- //check for constraints
- if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
- vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
+ if(vtx_distance_sq(vtx_moved,vtx_avoid)<16.0*vesicle->dmax){
continue;
}
+
+ memcpy((void *)&backupvtx[0],(void *)vtx_moved,sizeof(ts_vertex));
+ //move vertex in specified direction. first try, test move!
+ Rv=sqrt(pow(vtx_moved->x,2)+pow(vtx_moved->y,2)+pow(vtx_moved->z,2));
+ dh=2.0*Vol/(sqrt(3.0)*l0*l0);
+ vtx_moved->x=vtx_moved->x*(1.0-dh/Rv);
+ vtx_moved->y=vtx_moved->y*(1.0-dh/Rv);
+ vtx_moved->z=vtx_moved->z*(1.0-dh/Rv);
+
+
+//SKIPPING FIRST CHECK of CONSTRAINTS. This is not a final move.
+ //check for constraints
+/* if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
+ vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
+// continue;
+ break;
+ } */
// All checks OK!
- // doing second and final move.
- for(j=0;j<vtx_moved->neigh_no;j++){
- memcpy((void *)&backupvtx[j+1],(void *)vtx_moved->neigh[j],sizeof(ts_vertex));
- }
dvol=0.0;
for(j=0;j<vtx_moved->tristar_no;j++){
dvol-=vtx_moved->tristar[j]->volume;
+ }
+ volFirst=dvol;
+ for(j=0;j<vtx_moved->tristar_no;j++){
triangle_normal_vector(vtx_moved->tristar[j]);
dvol+=vtx_moved->tristar[j]->volume;
}
-
- voldiff=dvol-Vol;
-
+//TODO: here there is a bug. Don't know where, but preliminary success can
+//happen sometimes. And when I do this checks, constant value is not achieved
+//anymore
+/* voldiff=dvol-Vol;
if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
//calculate energy, return change in energy...
oenergy=vtx_moved->energy;
energy_vertex(vtx_moved);
delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
//the same is done for neighbouring vertices
- for(i=0;i<vtx_moved->neigh_no;i++){
- oenergy=vtx_moved->neigh[i]->energy;
- energy_vertex(vtx_moved->neigh[i]);
- delta_energy+=vtx_moved->neigh[i]->xk*(vtx_moved->neigh[i]->energy-oenergy);
+ for(j=0;j<vtx_moved->neigh_no;j++){
+ oenergy=vtx_moved->neigh[j]->energy;
+ energy_vertex(vtx_moved->neigh[j]);
+ delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy);
}
*retEnergy=delta_energy;
- vtx_backup=backupvtx;
+ *vtx_backup=backupvtx;
+ *vtx_moved_retval=vtx_moved;
+ fprintf(stderr, "Preliminary success\n");
return TS_SUCCESS;
- }
+ } */
+
+// fprintf(stderr, "Step 2 success\n");
//do it again ;)
dh=Vol*dh/dvol;
vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
@@ -80,31 +96,40 @@
vtx_moved->z=vtx_moved->z*(1-dh/Rv);
//check for constraints
if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
- for(j=0;j<vtx_moved->neigh_no;j++){
- memcpy((void *)vtx_moved->neigh[j],(void *)&backupvtx[j+1],sizeof(ts_vertex));
- }
vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
- continue;
+ //also, restore normals
+ for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]);
+// continue;
+ break;
}
+ dvol=volFirst;
+ for(j=0;j<vtx_moved->tristar_no;j++){
+ triangle_normal_vector(vtx_moved->tristar[j]);
+ dvol+=vtx_moved->tristar[j]->volume;
+ }
voldiff=dvol-Vol;
if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
//calculate energy, return change in energy...
+// fprintf(stderr, "Constvol success! %e\n",voldiff);
+ for(j=0;j<vtx_moved->neigh_no;j++){
+ memcpy((void *)&backupvtx[j+1],(void *)vtx_moved->neigh[j],sizeof(ts_vertex));
+ }
+
oenergy=vtx_moved->energy;
energy_vertex(vtx_moved);
delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
//the same is done for neighbouring vertices
- for(i=0;i<vtx_moved->neigh_no;i++){
- oenergy=vtx_moved->neigh[i]->energy;
- energy_vertex(vtx_moved->neigh[i]);
- delta_energy+=vtx_moved->neigh[i]->xk*(vtx_moved->neigh[i]->energy-oenergy);
+ for(j=0;j<vtx_moved->neigh_no;j++){
+ oenergy=vtx_moved->neigh[j]->energy;
+ energy_vertex(vtx_moved->neigh[j]);
+ delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy);
}
*retEnergy=delta_energy;
- vtx_backup=backupvtx;
+ *vtx_backup=backupvtx;
+ *vtx_moved_retval=vtx_moved;
return TS_SUCCESS;
}
-
-
}
free(backupvtx);
return TS_FAIL;
@@ -145,16 +170,26 @@
ts_bool constvolumerestore(ts_vertex *vtx_moved,ts_vertex *vtx_backup){
ts_uint j;
- for(j=0;j<vtx_moved->neigh_no;j++){
- memcpy((void *)vtx_moved->neigh[j],(void *)&vtx_backup[j+1],sizeof(ts_vertex));
- }
- vtx_moved=memcpy((void *)vtx_moved,(void *)&vtx_backup[0],sizeof(ts_vertex));
+ memcpy((void *)vtx_moved,(void *)&vtx_backup[0],sizeof(ts_vertex));
+ for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]);
+ for(j=0;j<vtx_moved->neigh_no;j++){
+ // memcpy((void *)vtx_moved->neigh[j],(void *)&vtx_backup[j+1],sizeof(ts_vertex));
+ energy_vertex(vtx_moved->neigh[j]);
+ }
+ free(vtx_backup);
return TS_SUCCESS;
}
-ts_bool constvolumeaccept(ts_vertex *vtx_moved, ts_vertex *vtx_backup){
-
+ts_bool constvolumeaccept(ts_vesicle *vesicle,ts_vertex *vtx_moved, ts_vertex *vtx_backup){
+ ts_bool retval;
+ ts_uint cellidx=vertex_self_avoidance(vesicle, vtx_moved);
+ if(vtx_moved->cell!=vesicle->clist->cell[cellidx]){
+ retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx_moved);
+ if(retval==TS_SUCCESS) cell_remove_vertex(vtx_backup[0].cell,vtx_moved);
+
+ }
+ free(vtx_backup);
return TS_SUCCESS;
}
--
Gitblit v1.9.3