Merge branch 'trout-rbc' of https://bitbucket.org/samop/trisurf-ng into tro...

mihaf

2014-04-11 091f7d56b288b83a32a2f7286ce532e5471ed49b

Merge branch 'trout-rbc' of https://bitbucket.org/samop/trisurf-ng into trout-rbc

Conflicts:
	src/tape

 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

 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;

}

 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;    

 src/io.h

@@ -35,6 +35,7 @@
    long int inititer;
    long int mcsweeps;
    long int quiet;
    long int shc;
} ts_tape;

typedef struct{

 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++;

 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;
}

 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

 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;
}

 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;
}

 src/tape

@@ -51,6 +51,9 @@
iterations=500


###### Spherical harmonics ###########
spherical_harmonics_coefficients=21

#shut up if we are using cluster!!!
quiet=false


 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);
        }

 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

 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;
}