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 | 106 ++++++++++++++++++++++++++++-------------------------
1 files changed, 56 insertions(+), 50 deletions(-)
diff --git a/src/constvol.c b/src/constvol.c
index 33922b4..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>
@@ -12,76 +13,81 @@
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;
-/* for(k=0;k<vtx_moved->neigh[j]->neigh_no;k++){
- if(vtx_moved->neigh[j]->neigh[k]==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. */
+ 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);
-// fprintf(stderr,"Prej (x,y,z)=(%e,%e,%e).\n",vtx_moved->x,vtx_moved->y,vtx_moved->z);
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);
-// fprintf(stderr,"Potem (x,y,z)=(%e,%e,%e). Vol=%e\n",vtx_moved->x,vtx_moved->y,vtx_moved->z,Vol);
+
+//SKIPPING FIRST CHECK of CONSTRAINTS. This is not a final move.
//check for constraints
- if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
+/* if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
- continue;
- }
-// fprintf(stderr,"Sprejet.\n");
-
+// continue;
+ break;
+ } */
// All checks OK!
- fprintf(stderr, "Step 1 success\n");
- 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_moved_retval=vtx_moved;
fprintf(stderr, "Preliminary success\n");
return TS_SUCCESS;
- }
- fprintf(stderr, "Step 2 success\n");
+ } */
+
+// 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));
@@ -90,45 +96,42 @@
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=0.0;
+ dvol=volFirst;
for(j=0;j<vtx_moved->tristar_no;j++){
- dvol-=vtx_moved->tristar[j]->volume;
triangle_normal_vector(vtx_moved->tristar[j]);
dvol+=vtx_moved->tristar[j]->volume;
}
-
- fprintf(stderr, "Step 3a success voldiff=%e\n",voldiff);
voldiff=dvol-Vol;
- fprintf(stderr, "Step 3b success voldiff=%e\n",voldiff);
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_moved_retval=vtx_moved;
- fprintf(stderr, "DVOL=%e\n",voldiff);
return TS_SUCCESS;
}
-
-
}
free(backupvtx);
- fprintf(stderr, "fail\n");
return TS_FAIL;
}
@@ -168,9 +171,12 @@
ts_bool constvolumerestore(ts_vertex *vtx_moved,ts_vertex *vtx_backup){
ts_uint j;
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));
+ // 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;
}
--
Gitblit v1.9.3