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 | 203 ++++++++++++++++++++++++++++++++++++--------------
1 files changed, 144 insertions(+), 59 deletions(-)
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 0f59aec..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>
@@ -9,28 +10,110 @@
#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);
- vesicle->vlist = init_sort_neighbours(vesicle->blist,vesicle->vlist);
+ 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(stderr,"initial_distribution finished!\n");
+ 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){
@@ -64,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;
}
}
}
@@ -105,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++){
@@ -166,15 +249,16 @@
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_cneighbour(blist,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_cneighbour(blist,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;
@@ -206,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;
}
@@ -232,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]);
}
}
}
@@ -254,7 +339,7 @@
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)){
@@ -334,11 +419,11 @@
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]->vertex[0];
k4=tria[k]->vertex[1];
--
Gitblit v1.9.3