From 632960fbf1ebde72860966058354f7a6b61f5267 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Thu, 10 Apr 2014 09:35:04 +0000
Subject: [PATCH] Added possibility to skip calculations of spherical harmonics. Fixed documentation in triangle.c.
---
src/main.c | 8 ++
src/timestep.c | 22 ++++---
src/vesicle.c | 2
src/initial_distribution.c | 7 ++
src/sh.c | 3
src/triangle.c | 96 +++++++++++++++++---------------
6 files changed, 80 insertions(+), 58 deletions(-)
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 7d5647c..5006a64 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -93,7 +93,12 @@
vesicle->pressure= tape->pressure;
vesicle->pswitch=tape->pswitch;
- vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
+ if(tape->shc>0){
+ vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
+ }
+ else {
+ vesicle->sphHarmonics=NULL;
+ }
return vesicle;
}
diff --git a/src/main.c b/src/main.c
index 63b72cf..9bb3430 100644
--- a/src/main.c
+++ b/src/main.c
@@ -51,7 +51,13 @@
vesicle->dmax=tape->dmax*tape->dmax;
poly_assign_filament_xi(vesicle,tape);
vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
- vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
+ /* 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 9f28d2f..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));
}
diff --git a/src/timestep.c b/src/timestep.c
index c61099b..b7857dd 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -22,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);
@@ -45,17 +45,19 @@
write_master_xml_file("test.pvd");
epochtime=get_epoch();
gyration_eigen(vesicle, &l1, &l2, &l3);
- get_area_volume(vesicle, &area,&volume);
- vesicle_volume(vesicle);
+ 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);
- preparationSh(vesicle,r0);
- calculateYlmi(vesicle);
- calculateUlm(vesicle);
- storeUlm2(vesicle);
- saveAvgUlm2(vesicle);
+ if(vesicle->sphHarmonics!=NULL){
+ preparationSh(vesicle,r0);
+ calculateYlmi(vesicle);
+ calculateUlm(vesicle);
+ storeUlm2(vesicle);
+ saveAvgUlm2(vesicle);
+ }
- fprintf(fd, "%lu %u %e %e %e %e %e %e %e\n",epochtime,i,vmsr,bfsr,volume, area,l1,l2,l3);
-
+ 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