Merge branch 'nirgov' of bitbucket.org:samop/trisurf-ng into nirgov

Samo Penic

2016-07-05 dd4e182bc818974b3ad529422e631aa6a0764f2f

Merge branch 'nirgov' of bitbucket.org:samop/trisurf-ng into nirgov

 src/Makefile.am

@@ -1,4 +1,4 @@
bin_PROGRAMS = trisurf tsmeasure
bin_PROGRAMS = trisurf tsmeasure tspoststat
trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c cluster.c
GITVERSION:=$(shell git --no-pager describe --tags --always --dirty)
AM_CFLAGS = -Wall -Werror -DTS_VERSION=\"$(GITVERSION)\" -fgnu89-inline
@@ -20,5 +20,8 @@

tsmeasure_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c tsmeasure.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c cluster.c
tsmeasure_LDADD = ${libcurl_LIBS} ${libxml2_LIBS}

tspoststat_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c tspoststat.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c cluster.c
tspoststat_LDADD = ${libcurl_LIBS} ${libxml2_LIBS}
#gitversion.c: .git/HEAD .git/index
#    echo "const char *gitversion = \"$(shell git rev-parse HEAD)\";" > $@

 src/initial_distribution.c

@@ -44,6 +44,7 @@
    vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
        vesicle->tape=tape;
    set_vesicle_values_from_tape(vesicle);
        initial_population_with_c0(vesicle,tape);
    return vesicle;
}

@@ -111,7 +112,14 @@
        vesicle->sphHarmonics=NULL;
    }

    int rndvtx;
    
    return TS_SUCCESS;

}


ts_bool initial_population_with_c0(ts_vesicle *vesicle, ts_tape *tape){
    int rndvtx,i,j;
    if(tape->number_of_vertices_with_c0>0){
        ts_fprintf(stderr,"Setting values for spontaneous curvature as defined in tape\n");
        j=0;
@@ -133,13 +141,8 @@
            sweep_attraction_bond_energy(vesicle);
        }
    }
    
    return TS_SUCCESS;

}





ts_bool pentagonal_dipyramid_vertex_distribution(ts_vertex_list *vlist){

 src/initial_distribution.h

@@ -15,6 +15,8 @@

ts_vesicle *create_vesicle_from_tape(ts_tape *tape);
ts_bool set_vesicle_values_from_tape(ts_vesicle *vesicle);

ts_bool initial_population_with_c0(ts_vesicle *vesicle, ts_tape *tape);
/** Sets the initial position of the vertexes to dipyramid
 *
 *      @param *vlist is a pointer to list of vertices

 src/main.c

@@ -127,6 +127,7 @@
    }
            //printf("nucleus coords: %.17e %.17e %.17e\n",vesicle->nucleus_center[0], vesicle->nucleus_center[1], vesicle->nucleus_center[2]);

            //write_vertex_xml_file(vesicle,1000);
    run_simulation(vesicle, tape->mcsweeps, tape->inititer, tape->iterations, start_iteration);
    write_master_xml_file(command_line_args.output_fullfilename);
    write_dout_fcompat_file(vesicle,"dout");

 src/restore.c

@@ -54,6 +54,10 @@
                if ((!xmlStrcmp(cur1->name, (const xmlChar *)"Piece"))){
                    cur2=cur1->xmlChildrenNode;
                    while(cur2!=NULL){
                        if ((!xmlStrcmp(cur2->name, (const xmlChar *)"PointData"))){
                            if(vesicle!=NULL)
                                parseXMLPointData(vesicle,doc,cur2);
                        }
                        if ((!xmlStrcmp(cur2->name, (const xmlChar *)"Points"))){
                            //fprintf(stderr,"Found point data\n");
                            if(vesicle!=NULL)
@@ -81,6 +85,7 @@

    init_normal_vectors(vesicle->tlist);
    mean_curvature_and_energy(vesicle);
    sweep_attraction_bond_energy(vesicle);

/* TODO: filaments */

@@ -182,7 +187,6 @@

    vesicle->tape=tape;
    set_vesicle_values_from_tape(vesicle);

    return vesicle;
}

@@ -318,7 +322,48 @@
    return TS_SUCCESS;
}

/* this parses the data for vertices (like spontaneous curvature, etc.) */
ts_bool parseXMLPointData(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
    xmlNodePtr child = cur->xmlChildrenNode;
    xmlChar *property_name;
    xmlChar *values;
    char *vals;
    char *token;
    int idx, polyidx, monoidx, filidx, fonoidx;
    while (child != NULL) {
        if ((!xmlStrcmp(child->name, (const xmlChar *)"DataArray"))){
            property_name=xmlGetProp(child, (xmlChar *)"Name");
        //    fprintf(stderr,"Name: %s\n", property_name);
            if(!xmlStrcmp(property_name,(const xmlChar *)"spontaneous_curvature")){
                values=xmlNodeListGetString(doc,child->xmlChildrenNode,1);
                vals=(char *)values;
                token=strtok(vals," ");
                idx=0;
                while(token!=NULL){
                    if(idx<vesicle->vlist->n){
                        vesicle->vlist->vtx[idx]->c=atof(token);
                    } else if(vesicle->tape->nmono && vesicle->tape->npoly && idx<vesicle->vlist->n+vesicle->tape->nmono*vesicle->tape->npoly) {
                        polyidx=(idx-vesicle->vlist->n)/vesicle->tape->nmono;
                        monoidx=(idx-vesicle->vlist->n)%vesicle->tape->nmono;
                        vesicle->poly_list->poly[polyidx]->vlist->vtx[monoidx]->c=atof(token);
                    } else {
                        filidx=(idx-vesicle->vlist->n-vesicle->tape->nmono*vesicle->tape->npoly)/vesicle->tape->nfono;
                        fonoidx=(idx-vesicle->vlist->n-vesicle->tape->nmono*vesicle->tape->npoly)%vesicle->tape->nfono;
                        //fprintf(stderr,"filidx=%d, fonoidx=%d, coord=%s,%s,%s\n",filidx,fonoidx,token[0],token[1],token[2]);
                        vesicle->filament_list->poly[filidx]->vlist->vtx[fonoidx]->c=atof(token);
                    }
                    idx++;
                    token=strtok(NULL," ");
                }
                xmlFree(values);		
            }
            xmlFree(property_name);
        }

        child=child->next;
    }
    return TS_SUCCESS;
}
/* this is a parser of vertex positions and bonds from main xml data */
ts_bool parseXMLVertexPosition(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
    xmlNodePtr child = cur->xmlChildrenNode;

 src/restore.h

@@ -9,6 +9,7 @@
ts_bool parseTrisurfTria(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);
ts_bool parseTrisurfTriaNeigh(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);
ts_bool parseTrisurfTristar(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);
ts_bool parseXMLPointData(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur);
ts_bool parseXMLVertexPosition(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur);
ts_bool parseXMLBonds(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur);
ts_bool parseTrisurfNucleus(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);

 src/tsmeasure.c

@@ -84,7 +84,6 @@
            tape_free(vesicle->tape);
            vesicle_free(vesicle);
                }
        free(ent);  
        }
    for (n = 0; n < count; n++)
      {

 src/tspoststat.c

New file
@@ -0,0 +1,119 @@
/* vim: set ts=4 sts=4 sw=4 noet : */
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include "general.h"
//#include "vertex.h"
//#include "bond.h"
//#include "triangle.h"
//#include "cell.h"
#include "vesicle.h"
#include "io.h"
//#include "initial_distribution.h"
//#include "frame.h"
//#include "timestep.h"
//#include "poly.h"
#include "stats.h"
#include "sh.h"
#include "shcomplex.h"
#include "dumpstate.h"
#include "restore.h"
#include <string.h>
#include <getopt.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <dirent.h>
#include <errno.h>
#include <snapshot.h>
#include<gsl/gsl_complex.h>
#include<gsl/gsl_complex_math.h>


ts_vesicle *restoreVesicle(char *filename){
    ts_vesicle *vesicle = parseDump(filename);
    return vesicle;
}

void vesicle_calculate_ulm2(ts_vesicle *vesicle){
    //complex_sph_free(vesicle->sphHarmonics);

    //vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,21);
    vesicle_volume(vesicle);
    preparationSh(vesicle,getR0(vesicle));
    calculateUlmComplex(vesicle);
    ts_int i,j;
    for(i=0;i<vesicle->sphHarmonics->l;i++){
            for(j=i;j<2*i+1;j++){
            printf("%e ", gsl_complex_abs2(vesicle->sphHarmonics->ulmComplex[i][j]));
            }
    }
        printf("\n");

}



int count_bonds_with_energy(ts_bond_list *blist){

    unsigned int i, cnt;
    cnt=0;
    for(i=0;i<blist->n;i++){
        if(fabs(blist->bond[i]->energy)>1e-16) cnt++;
    }
    return cnt;
}

int main(){
    ts_vesicle *vesicle;
    ts_char *i,*j;
    ts_uint tstep,n;
        ts_char *number;
    struct dirent **list;
    ts_double l1,l2,l3;
    int count;
    ts_fprintf(stderr,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);

    fprintf(stdout, "OuterLoop Volume Area lamdba1 lambda2 lambda3 Nbw/Nb\n");


    count=scandir(".",&list,0,alphasort);
    if(count<0){
        fatal("Error, cannot open directory.",1);
    }
        tstep=0;
    for(n=0;n<count;n++){
        struct dirent *ent;
        ent=list[n];	
                i=rindex(ent->d_name,'.');
                if(i==NULL) {
                continue;
        }
                if(strcmp(i+1,"vtu")==0){
                    j=rindex(ent->d_name,'_');
                    if(j==NULL) continue;
                    number=strndup(j+1,j-i); 
            quiet=1;
                    ts_fprintf(stdout,"timestep: %u filename: %s\n",atoi(number),ent->d_name);
//            printf("%u ",atoi(number));
            vesicle=restoreVesicle(ent->d_name);
//            vesicle_calculate_ulm2(vesicle);
            vesicle_volume(vesicle);
            vesicle_area(vesicle);
            gyration_eigen(vesicle,&l1,&l2,&l3);
            fprintf(stdout,"%d %.17e %.17e %.17e %.17e %.17e %.17e\n",atoi(number),vesicle->volume, vesicle->area,l1,l2,l3, (ts_double)count_bonds_with_energy(vesicle->blist)/(ts_double)vesicle->blist->n),
                        tstep++;

                    free(number);
            tape_free(vesicle->tape);
            vesicle_free(vesicle);
                }
        }
    for (n = 0; n < count; n++)
      {
          free(list[n]);
      }
	
    free(list);
    return 0;
}