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/initial_distribution.c | 287 +++++++++++++++++++++++++++++++++++++--------------------
1 files changed, 187 insertions(+), 100 deletions(-)
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index d5f1e94..7dfd23d 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -1,3 +1,4 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
#include<stdlib.h>
#include<math.h>
#include<stdio.h>
@@ -8,25 +9,111 @@
#include "vertex.h"
#include "triangle.h"
#include "initial_distribution.h"
+#include "energy.h"
+#include "poly.h"
+#include "io.h"
+#include "sh.h"
+#include "shcomplex.h"
ts_vesicle *initial_distribution_dipyramid(ts_uint nshell, ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
- ts_fprintf(stderr,"Starting initial_distribution on vesicle with %u shells!...\n",nshell);
+ ts_fprintf(stdout,"Starting initial_distribution on vesicle with %u shells!...\n",nshell);
ts_bool retval;
- ts_uint no_vertices=5*nshell*nshell+2;
-
- ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize);
- vesicle->nshell=nshell;
- retval = vtx_set_global_values(vesicle);
- retval = pentagonal_dipyramid_vertex_distribution(vesicle->vlist);
- retval = init_vertex_neighbours(vesicle->vlist);
- retval = init_sort_neighbours(vesicle->vlist);
- retval = init_vesicle_bonds(vesicle);
- retval = init_triangles(vesicle);
- retval = init_triangle_neighbours(vesicle);
- retval = init_common_vertex_triangle_neighbours(vesicle);
- ts_fprintf(stderr,"initial_distribution finished!\n");
+ ts_uint no_vertices=5*nshell*nshell+2;
+ ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize);
+ vesicle->nshell=nshell;
+ //retval = vtx_set_global_values(vesicle);
+ retval = pentagonal_dipyramid_vertex_distribution(vesicle->vlist);
+ retval = init_vertex_neighbours(vesicle->vlist);
+ vesicle->vlist = init_sort_neighbours(vesicle->blist,vesicle->vlist);
+ // retval = init_vesicle_bonds(vesicle); // bonds are created in sort_neigh
+ retval = init_triangles(vesicle);
+ retval = init_triangle_neighbours(vesicle);
+ retval = init_common_vertex_triangle_neighbours(vesicle);
+ retval = init_normal_vectors(vesicle->tlist);
+ retval = mean_curvature_and_energy(vesicle);
+ ts_fprintf(stdout,"initial_distribution finished!\n");
+ if(retval);
return vesicle;
}
+
+
+
+ts_vesicle *create_vesicle_from_tape(ts_tape *tape){
+ ts_vesicle *vesicle;
+
+ vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
+ vesicle->tape=tape;
+ set_vesicle_values_from_tape(vesicle);
+ return vesicle;
+}
+
+ts_bool set_vesicle_values_from_tape(ts_vesicle *vesicle){
+ // Nucleus:
+ ts_vertex *vtx;
+ ts_tape *tape=vesicle->tape;
+ vesicle->R_nucleus=tape->R_nucleus*tape->R_nucleus;
+
+ vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
+
+ //Initialize grafted polymers (brush):
+ vesicle->poly_list=init_poly_list(tape->npoly,tape->nmono, vesicle->vlist, vesicle);
+ vesicle->spring_constant=tape->kspring;
+ poly_assign_spring_const(vesicle);
+
+ //Initialize filaments (polymers inside the vesicle):
+ vesicle->filament_list=init_poly_list(tape->nfil,tape->nfono, NULL, vesicle);
+ poly_assign_filament_xi(vesicle,tape);
+
+ ts_uint i,j;
+ for(i=0;i<vesicle->filament_list->n;i++){
+ for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
+ bond_vector(vesicle->filament_list->poly[i]->blist->bond[j]);
+ vesicle->filament_list->poly[i]->blist->bond[j]->bond_length = sqrt(vtx_distance_sq(vesicle->filament_list->poly[i]->blist->bond[j]->vtx1,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2));
+ }
+ }
+
+ for(i=0;i<vesicle->filament_list->n;i++){
+ for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
+ vtx = vesicle->filament_list->poly[i]->vlist->vtx[j];
+ if(vtx->bond_no == 2){
+ vtx->energy = -(vtx->bond[0]->x*vtx->bond[1]->x + vtx->bond[0]->y*vtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length;
+ }
+ }
+ }
+
+ for(i=0;i<vesicle->filament_list->n;i++){
+ vertex_list_assign_id(vesicle->filament_list->poly[i]->vlist,TS_ID_FILAMENT);
+ }
+
+// vesicle->spring_constant=tape->kspring;
+// poly_assign_spring_const(vesicle);
+
+
+ vesicle->nshell=tape->nshell;
+ vesicle->dmax=tape->dmax*tape->dmax; /* dmax^2 in the vesicle dmax variable */
+ vesicle->bending_rigidity=tape->xk0;
+ vtx_set_global_values(vesicle); /* make xk0 default value for every vertex */
+ ts_fprintf(stdout, "Tape setting: xk0=%e\n",tape->xk0);
+ vesicle->stepsize=tape->stepsize;
+ vesicle->clist->ncmax[0]=tape->ncxmax;
+ vesicle->clist->ncmax[1]=tape->ncymax;
+ vesicle->clist->ncmax[2]=tape->nczmax;
+ vesicle->clist->max_occupancy=8; /* hard coded max occupancy? */
+
+ vesicle->pressure= tape->pressure;
+ vesicle->pswitch=tape->pswitch;
+ if(tape->shc>0){
+ vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc);
+ }
+ else {
+ vesicle->sphHarmonics=NULL;
+ }
+ return TS_SUCCESS;
+
+}
+
+
+
ts_bool pentagonal_dipyramid_vertex_distribution(ts_vertex_list *vlist){
@@ -60,40 +147,40 @@
ts_double dx,dy; // end loop prereq
/* topmost vertex */
- vtx[1]->data->x=0.0;
- vtx[1]->data->y=0.0;
- vtx[1]->data->z=z0*(ts_double)nshell;
+ vtx[1]->x=0.0;
+ vtx[1]->y=0.0;
+ vtx[1]->z=z0*(ts_double)nshell;
/* starting from to in circular order on pentagrams */
for(i=1;i<=nshell;i++){
n0=2+5*i*(i-1)/2; //-1 would be for the reason that C index starts from 0
- vtx[n0]->data->x=0.0;
- vtx[n0]->data->y=(ts_double)i*xl0;
- vtx[n0+i]->data->x=vtx[n0]->data->y*s1;
- vtx[n0+i]->data->y=vtx[n0]->data->y*c1;
- vtx[n0+2*i]->data->x=vtx[n0]->data->y*s2;
- vtx[n0+2*i]->data->y=vtx[n0]->data->y*c2;
- vtx[n0+3*i]->data->x=-vtx[n0+2*i]->data->x;
- vtx[n0+3*i]->data->y=vtx[n0+2*i]->data->y;
- vtx[n0+4*i]->data->x=-vtx[n0+i]->data->x;
- vtx[n0+4*i]->data->y=vtx[n0+i]->data->y;
+ vtx[n0]->x=0.0;
+ vtx[n0]->y=(ts_double)i*xl0;
+ vtx[n0+i]->x=vtx[n0]->y*s1;
+ vtx[n0+i]->y=vtx[n0]->y*c1;
+ vtx[n0+2*i]->x=vtx[n0]->y*s2;
+ vtx[n0+2*i]->y=vtx[n0]->y*c2;
+ vtx[n0+3*i]->x=-vtx[n0+2*i]->x;
+ vtx[n0+3*i]->y=vtx[n0+2*i]->y;
+ vtx[n0+4*i]->x=-vtx[n0+i]->x;
+ vtx[n0+4*i]->y=vtx[n0+i]->y;
}
/* vertexes on the faces of the dipyramid */
for(i=1;i<=nshell;i++){
n0=2+5*i*(i-1)/2; // -1 would be because of C!
for(j=1;j<=i-1;j++){
- dx=(vtx[n0]->data->x-vtx[n0+4*i]->data->x)/(ts_double)i;
- dy=(vtx[n0]->data->y-vtx[n0+4*i]->data->y)/(ts_double)i;
- vtx[n0+4*i+j]->data->x=(ts_double)j*dx+vtx[n0+4*i]->data->x;
- vtx[n0+4*i+j]->data->y=(ts_double)j*dy+vtx[n0+4*i]->data->y;
+ dx=(vtx[n0]->x-vtx[n0+4*i]->x)/(ts_double)i;
+ dy=(vtx[n0]->y-vtx[n0+4*i]->y)/(ts_double)i;
+ vtx[n0+4*i+j]->x=(ts_double)j*dx+vtx[n0+4*i]->x;
+ vtx[n0+4*i+j]->y=(ts_double)j*dy+vtx[n0+4*i]->y;
}
for(k=0;k<=3;k++){ // I would be worried about zero starting of for
- dx=(vtx[n0+(k+1)*i]->data->x - vtx[n0+k*i]->data->x)/(ts_double) i;
- dy=(vtx[n0+(k+1)*i]->data->y - vtx[n0+k*i]->data->y)/(ts_double) i;
+ dx=(vtx[n0+(k+1)*i]->x - vtx[n0+k*i]->x)/(ts_double) i;
+ dy=(vtx[n0+(k+1)*i]->y - vtx[n0+k*i]->y)/(ts_double) i;
for(j=1; j<=i-1;j++){
- vtx[n0+k*i+j]->data->x= (ts_double)j*dx+vtx[n0+k*i]->data->x;
- vtx[n0+k*i+j]->data->y= (ts_double)j*dy+vtx[n0+k*i]->data->y;
+ vtx[n0+k*i+j]->x= (ts_double)j*dx+vtx[n0+k*i]->x;
+ vtx[n0+k*i+j]->y= (ts_double)j*dy+vtx[n0+k*i]->y;
}
}
}
@@ -101,15 +188,15 @@
for(i=1;i<=nshell;i++){
n0= 2+ 5*i*(i-1)/2;
for(j=0;j<=5*i-1;j++){
- vtx[n0+j]->data->z= z0*(ts_double)(nshell-i); // I would be worried about zero starting of for
+ vtx[n0+j]->z= z0*(ts_double)(nshell-i); // I would be worried about zero starting of for
}
}
/* for botom part of dipyramide we calculate the positions of vertices */
for(i=2+5*nshell*(nshell+1)/2;i<=vlist->n;i++){
- vtx[i]->data->x=vtx[vlist->n - i +1]->data->x;
- vtx[i]->data->y=vtx[vlist->n - i +1]->data->y;
- vtx[i]->data->z=-vtx[vlist->n - i +1]->data->z;
+ vtx[i]->x=vtx[vlist->n - i +1]->x;
+ vtx[i]->y=vtx[vlist->n - i +1]->y;
+ vtx[i]->z=-vtx[vlist->n - i +1]->z;
}
for(i=1;i<=vlist->n;i++){
@@ -148,29 +235,30 @@
return TS_SUCCESS;
}
-// TODO: with new datastructure can be rewritten.
-ts_bool init_sort_neighbours(ts_vertex_list *vlist){
+// TODO: with new datastructure can be rewritten. Partially it is done, but it is complicated.
+ts_vertex_list *init_sort_neighbours(ts_bond_list *blist,ts_vertex_list *vlist){
ts_vertex **vtx=vlist->vtx -1; // take a look at dipyramid function for comment.
ts_uint i,l,j,jj,jjj,k=0;
ts_double eps=0.001; // Take a look if EPS from math.h can be used
/*lets initialize memory for temporary vertex_list. Should we write a function instead */
- ts_vertex_list *tvlist=init_vertex_list(vlist->n);
+ ts_vertex_list *tvlist=vertex_list_copy(vlist);
ts_vertex **tvtx=tvlist->vtx -1; /* again to compensate for 0-indexing */
ts_double dist2; // Square of distance of neighbours
ts_double direct; // Something, dont know what, but could be normal of some kind
for(i=1;i<=vlist->n;i++){
k++; // WHY i IS NOT GOOD??
- vtx_add_neighbour(tvtx[k], tvtx[vtx[i]->data->neigh[0]->idx+1]); //always add 1st
+ vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->neigh[0]->idx+1]); //always add 1st
jjj=1;
jj=1;
- for(l=2;l<=vtx[i]->data->neigh_no;l++){
- for(j=2;j<=vtx[i]->data->neigh_no;j++){
- dist2=vtx_distance_sq(vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
- direct=vtx_direct(vtx[i],vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
+ for(l=2;l<=vtx[i]->neigh_no;l++){
+ for(j=2;j<=vtx[i]->neigh_no;j++){
+ dist2=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
+ direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
+// TODO: check if fabs can be used with all floating point types!!
if( (fabs(dist2-A0*A0)<=eps) && (direct>0.0) && (j!=jjj) ){
- vtx_add_neighbour(tvtx[k],tvtx[vtx[i]->data->neigh[j-1]->idx+1]);
+ vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->neigh[j-1]->idx+1]);
jjj=jj;
jj=j;
break;
@@ -178,24 +266,19 @@
}
}
}
-
- for(i=1;i<=vlist->n;i++){
- for(j=1;j<=vtx[i]->data->neigh_no;j++){
- if(vtx[i]->data->neigh_no!=tvtx[i]->data->neigh_no){ //doesn't work with nshell=1!
-// fprintf(stderr,"data1=%u data2=%u\n",vtx[i]->data->neigh_no,tvtx[i]->data->neigh_no);
- fatal("Number of neighbours not the same in init_sort_neighbours.",4);
- }
- //we must correct the pointers in original to point to their
- //neighbours according to indexes. Must be sure not to do it any
- //other way! Also, we need to repair the collection of bonds...
- vtx[i]->data->neigh[j-1]=vtx[tvtx[i]->data->neigh[j-1]->idx+1];
- }
+/* We use the temporary vertex for our main vertices and we abandon main
+ * vertices, because their neighbours are not correctly ordered */
+ // tvtx=vlist->vtx;
+ // vlist->vtx=tvtx;
+ // tvlist->vtx=vtx;
+ vtx_list_free(vlist);
+/* Let's make a check if the number of bonds is correct */
+ if((blist->n)!=3*(tvlist->n-2)){
+ ts_fprintf(stderr,"Number of bonds is %u should be %u!\n", blist->n, 3*(tvlist->n-2));
+ fatal("Number of bonds is not 3*(no_vertex-2).",4);
}
- // Must free memory for temporary vertex array to avoid memory leak! HERE! NOW!
-// free_vertex(tvlist.vertex,tvlist.n);
- vtx_list_free(tvlist);
- return TS_SUCCESS;
+ return tvlist;
}
@@ -207,8 +290,8 @@
ts_uint i,j,k;
for(i=1;i<=vlist->n;i++){
for(j=i+1;j<=vlist->n;j++){
- for(k=0;k<vtx[i]->data->neigh_no;k++){ // has changed 0 to < instead of 1 and <=
- if(vtx[i]->data->neigh[k]==vtx[j]){ //if addresses matches it is the same
+ for(k=0;k<vtx[i]->neigh_no;k++){ // has changed 0 to < instead of 1 and <=
+ if(vtx[i]->neigh[k]==vtx[j]){ //if addresses matches it is the same
bond_add(blist,vtx[i],vtx[j]);
break;
}
@@ -233,14 +316,15 @@
ts_double eps=0.001; // can we use EPS from math.h?
k=0;
for(i=1;i<=vesicle->vlist->n;i++){
- for(j=1;j<=vtx[i]->data->neigh_no;j++){
- for(jj=1;jj<=vtx[i]->data->neigh_no;jj++){
+ for(j=1;j<=vtx[i]->neigh_no;j++){
+ for(jj=1;jj<=vtx[i]->neigh_no;jj++){
// ts_fprintf(stderr,"%u: (%u,%u) neigh_no=%u ",i,j,jj,vtx[i].neigh_no);
// ts_fprintf(stderr,"%e, %e",vtx[i].neigh[j-1]->x,vtx[i].neigh[jj-1]->x);
- dist=vtx_distance_sq(vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
- direct=vtx_direct(vtx[i],vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
- if(fabs(dist-A0*A0)<=eps && direct < 0.0 && vtx[i]->data->neigh[j-1]->idx+1 > i && vtx[i]->data->neigh[jj-1]->idx+1 >i){
- triangle_add(tlist,vtx[i],vtx[i]->data->neigh[j-1],vtx[i]->data->neigh[jj-1]);
+ dist=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
+ direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
+// TODO: same as above
+ if(fabs(dist-A0*A0)<=eps && direct < 0.0 && vtx[i]->neigh[j-1]->idx+1 > i && vtx[i]->neigh[jj-1]->idx+1 >i){
+ triangle_add(tlist,vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
}
}
}
@@ -251,11 +335,11 @@
for(i=0;i<tlist->n;i++){
k=0;
for(j=0;j<3;j++){
- if(tlist->tria[i]->data->vertex[j]!=NULL)
+ if(tlist->tria[i]->vertex[j]!=NULL)
k++;
}
if(k!=3){
- fatal("Some triangles has less than 3 vertices..",4);
+ fatal("Some triangles have less than 3 vertices..",4);
}
}
if(tlist->n!=2*(vesicle->vlist->n -2)){
@@ -275,14 +359,14 @@
ts_triangle **tria=tlist->tria -1;
nobo=0;
for(i=1;i<=tlist->n;i++){
- i1=tria[i]->data->vertex[0];
- i2=tria[i]->data->vertex[1];
- i3=tria[i]->data->vertex[2];
+ i1=tria[i]->vertex[0];
+ i2=tria[i]->vertex[1];
+ i3=tria[i]->vertex[2];
for(j=1;j<=tlist->n;j++){
if(j==i) continue;
- j1=tria[j]->data->vertex[0];
- j2=tria[j]->data->vertex[1];
- j3=tria[j]->data->vertex[2];
+ j1=tria[j]->vertex[0];
+ j2=tria[j]->vertex[1];
+ j3=tria[j]->vertex[2];
if((i1==j1 && i3==j2) || (i1==j2 && i3==j3) || (i1==j3 && i3==j1)){
triangle_add_neighbour(tria[i],tria[j]);
nobo++;
@@ -290,14 +374,14 @@
}
}
for(i=1;i<=tlist->n;i++){
- i1=tria[i]->data->vertex[0];
- i2=tria[i]->data->vertex[1];
- i3=tria[i]->data->vertex[2];
+ i1=tria[i]->vertex[0];
+ i2=tria[i]->vertex[1];
+ i3=tria[i]->vertex[2];
for(j=1;j<=tlist->n;j++){
if(j==i) continue;
- j1=tria[j]->data->vertex[0];
- j2=tria[j]->data->vertex[1];
- j3=tria[j]->data->vertex[2];
+ j1=tria[j]->vertex[0];
+ j2=tria[j]->vertex[1];
+ j3=tria[j]->vertex[2];
if((i1==j1 && i2==j3) || (i1==j3 && i2==j2) || (i1==j2 && i2==j1)){
triangle_add_neighbour(tria[i],tria[j]);
nobo++;
@@ -305,14 +389,14 @@
}
}
for(i=1;i<=tlist->n;i++){
- i1=tria[i]->data->vertex[0];
- i2=tria[i]->data->vertex[1];
- i3=tria[i]->data->vertex[2];
+ i1=tria[i]->vertex[0];
+ i2=tria[i]->vertex[1];
+ i3=tria[i]->vertex[2];
for(j=1;j<=tlist->n;j++){
if(j==i) continue;
- j1=tria[j]->data->vertex[0];
- j2=tria[j]->data->vertex[1];
- j3=tria[j]->data->vertex[2];
+ j1=tria[j]->vertex[0];
+ j2=tria[j]->vertex[1];
+ j3=tria[j]->vertex[2];
if((i2==j1 && i3==j3) || (i2==j3 && i3==j2) || (i2==j2 && i3==j1)){
triangle_add_neighbour(tria[i],tria[j]);
nobo++;
@@ -335,20 +419,23 @@
ts_triangle **tria=tlist->tria -1;
for(i=1;i<=vesicle->vlist->n;i++){
- for(j=1;j<=vtx[i]->data->neigh_no;j++){
- k1=vtx[i]->data->neigh[j-1];
+ for(j=1;j<=vtx[i]->neigh_no;j++){
+ k1=vtx[i]->neigh[j-1];
jp=j+1;
- if(j == vtx[i]->data->neigh_no) jp=1;
- k2=vtx[i]->data->neigh[jp-1];
+ if(j == vtx[i]->neigh_no) jp=1;
+ k2=vtx[i]->neigh[jp-1];
for(k=1;k<=tlist->n;k++){ // VERY NON-OPTIMAL!!! too many loops (vlist.n * vtx.neigh * tlist.n )!
- k3=tria[k]->data->vertex[0];
- k4=tria[k]->data->vertex[1];
- k5=tria[k]->data->vertex[2];
+ k3=tria[k]->vertex[0];
+ k4=tria[k]->vertex[1];
+ k5=tria[k]->vertex[2];
// ts_fprintf(stderr,"%u %u: k=(%u %u %u)\n",k1,k2,k3,k4,k5);
if((vtx[i]==k3 && k1==k4 && k2==k5) ||
(vtx[i]==k4 && k1==k5 && k2==k3) ||
(vtx[i]==k5 && k1==k3 && k2==k4)){
- // ts_fprintf(stderr, "Added to tristar! ");
+
+//TODO: probably something wrong with neighbour distribution.
+// if(vtx[i]==k3 || vtx[i]==k4 || vtx[i]==k5){
+ // if(i==6) ts_fprintf(stdout, "Vtx[%u] > Added to tristar!\n",i);
vertex_add_tristar(vtx[i],tria[k]);
}
}
--
Gitblit v1.9.3