From 091f7d56b288b83a32a2f7286ce532e5471ed49b Mon Sep 17 00:00:00 2001
From: mihaf <miha.fosnaric@gmail.com>
Date: Fri, 11 Apr 2014 08:10:45 +0000
Subject: [PATCH] Merge branch 'trout-rbc' of https://bitbucket.org/samop/trisurf-ng into trout-rbc
---
src/main.c | 8 ++
src/io.c | 1
src/tape | 3
src/vesicle.c | 2
src/initial_distribution.c | 8 +
src/spherical_trisurf_ff.c | 24 ------
src/io.h | 1
src/Makefile.am | 2
src/timestep.c | 21 ++++-
src/sh.h | 1
src/spherical_trisurf.c | 23 -----
src/sh.c | 28 ++++++
src/triangle.c | 96 ++++++++++++-----------
13 files changed, 119 insertions(+), 99 deletions(-)
diff --git a/src/Makefile.am b/src/Makefile.am
index 204da92..970dec2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,6 +1,6 @@
trisurfdir=../
trisurf_PROGRAMS = trisurf
-trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c
+trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c
#trisurf_LDFLAGS = -lm -lconfuse
shdiscoverdir=../
shdiscover_PROGRAMS= shdiscover
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 0587d18..5006a64 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -11,6 +11,7 @@
#include "energy.h"
#include "poly.h"
#include "io.h"
+#include "sh.h"
ts_vesicle *initial_distribution_dipyramid(ts_uint nshell, ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
ts_fprintf(stdout,"Starting initial_distribution on vesicle with %u shells!...\n",nshell);
@@ -92,7 +93,12 @@
vesicle->pressure= tape->pressure;
vesicle->pswitch=tape->pswitch;
-
+ if(tape->shc>0){
+ vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
+ }
+ else {
+ vesicle->sphHarmonics=NULL;
+ }
return vesicle;
}
diff --git a/src/io.c b/src/io.c
index a25b16c..d001b7e 100644
--- a/src/io.c
+++ b/src/io.c
@@ -988,6 +988,7 @@
CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
+ CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
CFG_END()
};
cfg_t *cfg;
diff --git a/src/io.h b/src/io.h
index 117c9ff..699c73e 100644
--- a/src/io.h
+++ b/src/io.h
@@ -35,6 +35,7 @@
long int inititer;
long int mcsweeps;
long int quiet;
+ long int shc;
} ts_tape;
typedef struct{
diff --git a/src/main.c b/src/main.c
index eef5835..9bb3430 100644
--- a/src/main.c
+++ b/src/main.c
@@ -11,6 +11,7 @@
#include "frame.h"
#include "timestep.h"
#include "poly.h"
+#include "sh.h"
/** Entrance function to the program
* @param argv is a number of parameters used in program call (including the program name
@@ -50,6 +51,13 @@
vesicle->dmax=tape->dmax*tape->dmax;
poly_assign_filament_xi(vesicle,tape);
vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
+ /* spherical harmonics */
+ if(tape->shc>0){
+ vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
+ }
+ else {
+ vesicle->sphHarmonics=NULL;
+ }
if(command_line_args.reset_iteration_count) start_iteration=tape->inititer;
else start_iteration++;
diff --git a/src/sh.c b/src/sh.c
index 0c4b06d..5ab56eb 100644
--- a/src/sh.c
+++ b/src/sh.c
@@ -49,6 +49,7 @@
ts_bool sph_free(ts_spharm *sph){
int i,j;
+ if(sph==NULL) return TS_FAIL;
for(i=0;i<sph->l;i++){
if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
@@ -189,7 +190,7 @@
K=-sqrt(1.0/(M_PI))*cos(m*fi);
}
- return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));
+ return K*sqrt((2.0*l+1.0)/2.0*(ts_double)(fac1/fac2))*plgndr(l,abs(m),cos(theta));
}
@@ -379,3 +380,28 @@
sph->N++;
return TS_SUCCESS;
}
+
+
+ts_bool saveAvgUlm2(ts_vesicle *vesicle){
+
+ FILE *fh;
+
+ fh=fopen("sph2out.dat", "w");
+ if(fh==NULL){
+ err("Cannot open file %s for writing");
+ return TS_FAIL;
+ }
+
+ ts_spharm *sph=vesicle->sphHarmonics;
+ ts_int i,j;
+ fprintf(fh,"l,\tm,\tulm^2avg\n");
+ for(i=0;i<sph->l;i++){
+ for(j=0;j<2*i+1;j++){
+ fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
+
+ }
+ fprintf(fh,"\n");
+ }
+ fclose(fh);
+ return TS_SUCCESS;
+}
diff --git a/src/sh.h b/src/sh.h
index 7853139..32df46b 100644
--- a/src/sh.h
+++ b/src/sh.h
@@ -13,4 +13,5 @@
ts_bool preparationSh(ts_vesicle *vesicle, ts_double r0);
ts_bool calculateYlmi(ts_vesicle *vesicle);
ts_bool calculateUlm(ts_vesicle *vesicle);
+ts_bool saveAvgUlm2(ts_vesicle *vesicle);
#endif
diff --git a/src/spherical_trisurf.c b/src/spherical_trisurf.c
index 6c5946f..0669981 100644
--- a/src/spherical_trisurf.c
+++ b/src/spherical_trisurf.c
@@ -104,26 +104,3 @@
-ts_bool saveAvgUlm2(ts_vesicle *vesicle){
-
- FILE *fh;
-
- fh=fopen("sph2out.dat", "w");
- if(fh==NULL){
- err("Cannot open file %s for writing");
- return TS_FAIL;
- }
-
- ts_spharm *sph=vesicle->sphHarmonics;
- ts_int i,j;
- fprintf(fh,"l,\tm,\tulm^2avg\n");
- for(i=0;i<sph->l;i++){
- for(j=0;j<2*i+1;j++){
- fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
-
- }
- fprintf(fh,"\n");
- }
- fclose(fh);
- return TS_SUCCESS;
-}
diff --git a/src/spherical_trisurf_ff.c b/src/spherical_trisurf_ff.c
index 00a4cc0..157b0e4 100644
--- a/src/spherical_trisurf_ff.c
+++ b/src/spherical_trisurf_ff.c
@@ -87,27 +87,3 @@
}
-
-ts_bool saveAvgUlm2(ts_vesicle *vesicle){
-
- FILE *fh;
-
- fh=fopen("sph2out.dat", "w");
- if(fh==NULL){
- err("Cannot open file %s for writing");
- return TS_FAIL;
- }
-
- ts_spharm *sph=vesicle->sphHarmonics;
- ts_int i,j;
- fprintf(fh,"l,\tm,\tulm^2avg\n");
- for(i=0;i<sph->l;i++){
- for(j=0;j<2*i+1;j++){
- fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
-
- }
- fprintf(fh,"\n");
- }
- fclose(fh);
- return TS_SUCCESS;
-}
diff --git a/src/tape b/src/tape
index 713e2b7..c1f23e1 100644
--- a/src/tape
+++ b/src/tape
@@ -51,6 +51,9 @@
iterations=500
+###### Spherical harmonics ###########
+spherical_harmonics_coefficients=21
+
#shut up if we are using cluster!!!
quiet=false
diff --git a/src/timestep.c b/src/timestep.c
index 9ab60fd..b7857dd 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -9,9 +9,12 @@
#include "frame.h"
#include "io.h"
#include "stats.h"
+#include "sh.h"
+#include "vesicle.h"
ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){
ts_uint i, j;
+ ts_double r0;
ts_double l1,l2,l3,volume=0.0,area=0.0,vmsr,bfsr, vmsrt, bfsrt;
ts_ulong epochtime;
// char filename[255];
@@ -19,7 +22,7 @@
if(fd==NULL){
fatal("Cannot open statistics.csv file for writing",1);
}
- fprintf(fd, "Epoch OuterLoop VertexMoveSucessRate BondFlipSuccessRate Volume Area lamdba1 lambda2 lmabda3\n");
+ fprintf(fd, "Epoch OuterLoop VertexMoveSucessRate BondFlipSuccessRate Volume Area lamdba1 lambda2 lambda3\n");
centermass(vesicle);
cell_occupation(vesicle);
if(start_iteration<inititer) ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
@@ -42,9 +45,19 @@
write_master_xml_file("test.pvd");
epochtime=get_epoch();
gyration_eigen(vesicle, &l1, &l2, &l3);
- get_area_volume(vesicle, &area,&volume);
- fprintf(fd, "%lu %u %e %e %e %e %e %e %e\n",epochtime,i,vmsr,bfsr,volume, area,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).
+ r0=getR0(vesicle);
+ if(vesicle->sphHarmonics!=NULL){
+ preparationSh(vesicle,r0);
+ calculateYlmi(vesicle);
+ calculateUlm(vesicle);
+ storeUlm2(vesicle);
+ saveAvgUlm2(vesicle);
+ }
+
+ fprintf(fd, "%lu %u %e %e %1.16e %1.16e %1.16e %1.16e %1.16e\n",epochtime,i,vmsr,bfsr,volume, area,l1,l2,l3);
+ fflush(fd);
// sprintf(filename,"timestep-%05d.pov",i-inititer);
// write_pov_file(vesicle,filename);
}
diff --git a/src/triangle.c b/src/triangle.c
index 311e138..7dd5a6e 100644
--- a/src/triangle.c
+++ b/src/triangle.c
@@ -5,6 +5,7 @@
#include<math.h>
/** @brief Prepares the list for triangles.
+ * @returns pointer to empty data structure for maintaining triangle list.
*
* Create empty list for holding the information on triangles. Triangles are
* added later on with triangle_add().
@@ -27,14 +28,17 @@
}
/** @brief Add the triangle to the triangle list and create necessary data
- * structures.
+ * structures.
+ * @param *tlist is a pointer to triangle list where triangle should be created
+ * @param *vtx1, *vtx2, *vtx3 are the three vertices defining the triangle
+ * @returns pointer to the newly created triangle on success and NULL if
+ * triangle could not be created. It breaks program execution if memory
+ * allocation of triangle list can't be done.
*
- * Add the triangle ts_triangle with ts_triangle_data to the ts_triangle_list.
+ * Add the triangle ts_triangle to the ts_triangle_list.
* The triangle list is resized, the ts_triangle is allocated and
- * ts_triangle_data is allocated and zeroed. The function receives 4 arguments:
- * ts_triangle_list *tlist as list of triangles and 3 ts_vertex *vtx as
- * vertices that are used to form a triangle. Returns a pointer to newly
- * created triangle. This pointer doesn't need assigning, since it is
+ * triangle data is zeroed. Returned pointer to newly
+ * created triangle doesn't need assigning, since it is
* referenced by triangle list.
*
* WARNING: Function can be accelerated a bit by removing the NULL checks.
@@ -57,7 +61,6 @@
tlist->tria[tlist->n-1]=(ts_triangle *)calloc(1,sizeof(ts_triangle));
if(tlist->tria[tlist->n-1]==NULL) fatal("Cannot reallocate memory for additional ts_triangle.",5);
- // tlist->tria[tlist->n-1]->data=(ts_triangle_data *)calloc(1,sizeof(ts_triangle_data));
//NOW insert vertices!
tlist->tria[tlist->n - 1]->idx=tlist->n-1;
@@ -68,9 +71,14 @@
}
/** @brief Add the neigbour to triangles.
+ * @param *tria is a first triangle.
+ * @param *ntria is a second triangle.
+ * @returns TS_SUCCES on sucessful adition to the list, TS_FAIL if triangles
+ * are NULL and breaks execution FATALY if memory allocation error occurs.
*
* Add the neigbour to the list of neighbouring triangles. The
- * neighbouring triangles are those, who share two vertices. Function resizes
+ * neighbouring triangles are those, who share two vertices and corresponding
+ * bond. Function resizes
* the list and adds the pointer to neighbour. It receives two arguments of
* ts_triangle type. It then adds second triangle to the list of first
* triangle, but not the opposite. Upon
@@ -84,42 +92,40 @@
* debugging stupid NULL pointers.
*
* Example of usage:
- * triangle_remove_neighbour(tlist->tria[3], tlist->tria[4]);
+ * triangle_add_neighbour(tlist->tria[3], tlist->tria[4]);
*
- * Triangles 3 and 4 are not neighbours anymore.
+ * Triangle 4 is a neighbour of triangle 3, but (strangely) not the
+ * oposite. The function should be called again with the changed order of
+ * triangles to make neighbourship mutual.
*
*/
ts_bool triangle_add_neighbour(ts_triangle *tria, ts_triangle *ntria){
if(tria==NULL || ntria==NULL) return TS_FAIL;
-/*TODO: check if the neighbour already exists! Now there is no such check
- * because of the performance issue. */
tria->neigh_no++;
tria->neigh=realloc(tria->neigh,tria->neigh_no*sizeof(ts_triangle *));
if(tria->neigh == NULL)
fatal("Reallocation of memory failed during insertion of triangle neighbour in triangle_add_neighbour",3);
- tria->neigh[tria->neigh_no-1]=ntria;
-
-
-/* we repeat the procedure for the neighbour */
-/* ntria->data->neigh_no++;
- ntria->data->neigh=realloc(ntria->data->neigh,ntria->data->neigh_no*sizeof(ts_triangle *));
- if(ntria->data->neigh == NULL)
- fatal("Reallocation of memory failed during insertion of triangle neighbour in triangle_add_neighbour",3);
- ntria->data->neigh[ntria->data->neigh_no-1]=tria;
-*/
+ tria->neigh[tria->neigh_no-1]=ntria;
return TS_SUCCESS;
}
/** @brief Remove the neigbours from triangle.
+ * @param *tria is a first triangle.
+ * @param *ntria is neighbouring triangle.
+ * @returns TS_SUCCESS on successful removal, TS_FAIL if triangles are not
+ * neighbours and it breaks program execution FATALY if memory allocation
+ * problem occurs.
*
* Removes the neigbour from the list of neighbouring triangles. The
- * neighbouring triangles are those, who share two vertices. Function resizes
+ * neighbouring triangles are those, who share two vertices and corresponding
+ * bond. Function resizes
* the list and deletes the pointer to neighbour. It receives two arguments of
- * ts_triangle type. It then removes eachother form eachother's list. Upon
+ * ts_triangle type. It then mutually removes triangles from eachouther
+ * neighbour list. Upon
* success it returns TS_SUCCESS, upon failure to find the triangle in the
- * neighbour list returns TS_FAIL and it FATALY ends when the datastructure
- * cannot be resized.
+ * neighbour list returns TS_FAIL. It FATALY breaks program execution when the datastructure
+ * cannot be resized due to memory constrain problems.
*
* WARNING: The function doesn't check whether the pointer is NULL or invalid. It is the
* job of programmer to make sure the pointer is valid.
@@ -144,10 +150,8 @@
}
if(j==i) {
return TS_FAIL;
- //fatal("In triangle_remove_neighbour: Specified neighbour does not exist for given triangle",3);
}
tria->neigh_no--;
-// fprintf(stderr,"*** tria_number=%d\n",tria->neigh_no);
tria->neigh=(ts_triangle **)realloc(tria->neigh,tria->neigh_no*sizeof(ts_triangle *));
if(tria->neigh == NULL){
fprintf(stderr,"Ooops: tria->neigh_no=%d\n",tria->neigh_no);
@@ -163,10 +167,8 @@
}
if(j==i) {
return TS_FAIL;
- //fatal("In triangle_remove_neighbour: Specified neighbour does not exist for given triangle",3);
}
ntria->neigh_no--;
-// fprintf(stderr,"*** ntria_number=%d\n",ntria->neigh_no);
ntria->neigh=(ts_triangle **)realloc(ntria->neigh,ntria->neigh_no*sizeof(ts_triangle *));
if(ntria->neigh == NULL){
fprintf(stderr,"Ooops: ntria->neigh_no=%d\n",ntria->neigh_no);
@@ -176,25 +178,33 @@
}
-/** @brief Calculates normal vector of the triangle.
-
+/** @brief Calculates normal vector of the triangle, its corresponding area and volume.
+ * @param *tria is a triangle pointer for which normal, area and volume is
+ * to be calculated.
+ * @returns TS_SUCCESS on success. (always)
*
* Calculate normal vector of the triangle (xnorm, ynorm and znorm) and stores
- * information in underlying ts_triangle_data data_structure.
+ * information. At the same time
+ * triangle area is determined, since we already have the normal and volume of
+ * triangular pyramid with given triangle as a base and vesicle centroid as a
+ * tip.
*
* Function receives one argument of type ts_triangle. It should be corectly
- * initialized with underlying data structure of type ts_triangle_data. the
- * result is stored in triangle->data->xnorm, triangle->data->ynorm,
- * triangle->data->znorm. Returns TS_SUCCESS on completion.
+ * initialized. The
+ * result is stored in triangle->xnorm, triangle->ynorm, triangle->znorm.
+ * Area and volume are stored into triangle->area and triangle->volume.
+ * Returns TS_SUCCESS on completion.
*
- * NOTE: Function uses math.h library. pow function implementation is selected
- * accordind to the setting in genreal.h
+ * NOTE: Function uses math.h library. Function pow implementation is selected
+ * accordind to the used TS_DOUBLE_* definition set in general.h, so it should
+ * be compatible with any type of floating point precision.
*
* Example of usage:
* triangle_normal_vector(tlist->tria[3]);
*
* Computes normals and stores information into tlist->tria[3]->xnorm,
- * tlist->tria[3]->ynorm, tlist->tria[3]->znorm.
+ * tlist->tria[3]->ynorm, tlist->tria[3]->znorm tlist->tria[3]->area and
+ * tlist->tria[3]->volume.
*
*/
ts_bool triangle_normal_vector(ts_triangle *tria){
@@ -239,13 +249,9 @@
return TS_SUCCESS;
}
-
-
-
-
-
/** @brief Frees the memory allocated for data structure of triangle list
- * (ts_triangle_list)
+ * @param *tlist is a pointer to datastructure triangle list to be freed.
+ * @returns TS_SUCCESS on success (always).
*
* Function frees the memory of ts_triangle_list previously allocated. It
* accepts one argument, the address of data structure. It destroys all
diff --git a/src/vesicle.c b/src/vesicle.c
index 43b16b5..0a1b84b 100644
--- a/src/vesicle.c
+++ b/src/vesicle.c
@@ -6,6 +6,7 @@
#include "cell.h"
#include "stdlib.h"
#include "poly.h"
+#include "sh.h"
ts_vesicle *init_vesicle(ts_uint N, ts_uint ncmax1, ts_uint ncmax2, ts_uint
ncmax3, ts_double stepsize){
@@ -35,6 +36,7 @@
triangle_list_free(vesicle->tlist);
cell_list_free(vesicle->clist);
poly_list_free(vesicle->poly_list);
+ sph_free(vesicle->sphHarmonics);
free(vesicle);
return TS_SUCCESS;
}
--
Gitblit v1.9.3