Some rewritting done.

Samo Penic

2010-12-05 f74313919463dd08d93faeefd40900e2917c0157

Some rewritting done.

 src/Makefile.am

@@ -1,5 +1,5 @@
trisurfdir=../
trisurf_PROGRAMS=trisurf
trisurf_SOURCES=general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c main.c
trisurf_SOURCES=general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c main.c
trisurf_LDFLAGS = -lm -lconfuse
trisurf_CFLAGS = -Wall -Werror

 src/energy.c

@@ -6,14 +6,13 @@
#include<stdio.h>
ts_bool mean_curvature_and_energy(ts_vesicle *vesicle){

    ts_uint i, jj, jjp;
    ts_uint i;
    
    ts_vertex_list *vlist=&vesicle->vlist;
    ts_vertex *vtx=vlist->vertex;
    ts_vertex_list *vlist=vesicle->vlist;
    ts_vertex **vtx=vlist->vtx;

    for(i=0;i<vlist->n;i++){
        //should call with zero index!!!
        energy_vertex(&vtx[i]);
        energy_vertex(vtx[i]);
    }

    return TS_SUCCESS;
@@ -23,6 +22,7 @@
inline ts_bool energy_vertex(ts_vertex *vtx){
//    ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value!
//    ts_triangle *tristar=vtx->tristar-1;
    ts_vertex_data *data=vtx->data;
    ts_uint jj;
    ts_uint jjp,jjm;
    ts_vertex *j,*jp, *jm;
@@ -30,20 +30,20 @@
    ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0;
    ts_double x1,x2,x3,ctp,ctm,tot,xlen;
    ts_double h,ht;
    for(jj=1; jj<=vtx->neigh_no;jj++){
    for(jj=1; jj<=data->neigh_no;jj++){
        jjp=jj+1;
        if(jjp>vtx->neigh_no) jjp=1;
        if(jjp>data->neigh_no) jjp=1;
        jjm=jj-1;
        if(jjm<1) jjm=vtx->neigh_no;
        j=vtx->neigh[jj-1];
        jp=vtx->neigh[jjp-1];
        jm=vtx->neigh[jjm-1];
        jt=vtx->tristar[jj-1];
        x1=vertex_distance_sq(vtx,jp); //shouldn't be zero!
        x2=vertex_distance_sq(j,jp); // shouldn't be zero!
        x3=(j->x-jp->x)*(vtx->x-jp->x)+
           (j->y-jp->y)*(vtx->y-jp->y)+
           (j->z-jp->z)*(vtx->z-jp->z);
        if(jjm<1) jjm=data->neigh_no;
        j=data->neigh[jj-1];
        jp=data->neigh[jjp-1];
        jm=data->neigh[jjm-1];
        jt=data->tristar[jj-1];
        x1=vtx_distance_sq(vtx,jp); //shouldn't be zero!
        x2=vtx_distance_sq(j,jp); // shouldn't be zero!
        x3=(j->data->x-jp->data->x)*(data->x-jp->data->x)+
           (j->data->y-jp->data->y)*(data->y-jp->data->y)+
           (j->data->z-jp->data->z)*(data->z-jp->data->z);
        
#ifdef TS_DOUBLE_DOUBLE
        ctp=x3/sqrt(x1*x2-x3*x3);
@@ -54,11 +54,11 @@
#ifdef TS_DOUBLE_LONGDOUBLE
        ctp=x3/sqrtl(x1*x2-x3*x3);
#endif
        x1=vertex_distance_sq(vtx,jm);
        x2=vertex_distance_sq(j,jm);
        x3=(j->x-jm->x)*(vtx->x-jm->x)+
           (j->y-jm->y)*(vtx->y-jm->y)+
           (j->z-jm->z)*(vtx->z-jm->z);
        x1=vtx_distance_sq(vtx,jm);
        x2=vtx_distance_sq(j,jm);
        x3=(j->data->x-jm->data->x)*(data->x-jm->data->x)+
           (j->data->y-jm->data->y)*(data->y-jm->data->y)+
           (j->data->z-jm->data->z)*(data->z-jm->data->z);
#ifdef TS_DOUBLE_DOUBLE
        ctm=x3/sqrt(x1*x2-x3*x3);
#endif
@@ -70,26 +70,26 @@
#endif
        tot=ctp+ctm;
        tot=0.5*tot;
        xlen=vertex_distance_sq(j,vtx);
        xlen=vtx_distance_sq(j,vtx);
#ifdef  TS_DOUBLE_DOUBLE 
        vtx->bond_length[jj-1]=sqrt(xlen); 
        data->bond[jj-1]->data->bond_length=sqrt(xlen); 
#endif
#ifdef  TS_DOUBLE_FLOAT
        vtx->bond_length[jj-1]=sqrtf(xlen); 
        data->bond[jj-1]->data->bond_length=sqrtf(xlen); 
#endif
#ifdef  TS_DOUBLE_LONGDOUBLE 
        vtx->bond_length[jj-1]=sqrtl(xlen); 
        data->bond[jj-1]->data->bond_length=sqrtl(xlen); 
#endif

        vtx->bond_length_dual[jj-1]=tot*vtx->bond_length[jj-1];
        data->bond[jj-1]->data->bond_length_dual=tot*data->bond[jj-1]->data->bond_length;

        s+=tot*xlen;
        xh+=tot*(j->x - vtx->x);
        yh+=tot*(j->y - vtx->y);
        zh+=tot*(j->z - vtx->z);
        txn+=jt->xnorm;
        tyn+=jt->ynorm;
        tzn+=jt->znorm;
        xh+=tot*(j->data->x - data->x);
        yh+=tot*(j->data->y - data->y);
        zh+=tot*(j->data->z - data->z);
        txn+=jt->data->xnorm;
        tyn+=jt->data->ynorm;
        tzn+=jt->data->znorm;
    }
    
    h=xh*xh+yh*yh+zh*zh;
@@ -97,26 +97,27 @@
    s=s/4.0; 
#ifdef TS_DOUBLE_DOUBLE
    if(ht>=0.0) {
        vtx->curvature=sqrt(h);
        data->curvature=sqrt(h);
    } else {
        vtx->curvature=-sqrt(h);
        data->curvature=-sqrt(h);
    }
#endif
#ifdef TS_DOUBLE_FLOAT
    if(ht>=0.0) {
        vtx->curvature=sqrtf(h);
        data->curvature=sqrtf(h);
    } else {
        vtx->curvature=-sqrtf(h);
        data->curvature=-sqrtf(h);
    }
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
    if(ht>=0.0) {
        vtx->curvature=sqrtl(h);
        data->curvature=sqrtl(h);
    } else {
        vtx->curvature=-sqrtl(h);
        data->curvature=-sqrtl(h);
    }
#endif
    vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c);
//TODO: MAJOR!!!! What is vtx->data->c?????????????? Here it is 0!
    data->energy=0.5*s*(data->curvature/s-data->c)*(data->curvature/s-data->c);

    return TS_SUCCESS;
}

 src/frame.c

@@ -1,38 +1,38 @@
#include<stdlib.h>
#include "general.h"
#include "cell.h"

#include "frame.h"
ts_bool centermass(ts_vesicle *vesicle){
    ts_uint i;
    ts_uint i, n=vesicle->vlist->n;
    ts_vertex **vtx=vesicle->vlist->vtx;
    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;
    for(i=0;i<vesicle->vlist.n;i++){
        vesicle->cm[0]+=vesicle->vlist.vertex[i].x;
        vesicle->cm[1]+=vesicle->vlist.vertex[i].y;
        vesicle->cm[2]+=vesicle->vlist.vertex[i].z; 
    for(i=0;i<n;i++){
        vesicle->cm[0]+=vtx[i]->data->x;
        vesicle->cm[1]+=vtx[i]->data->y;
        vesicle->cm[2]+=vtx[i]->data->z; 
    } 
    vesicle->cm[0]/=(float)vesicle->vlist.n;
    vesicle->cm[1]/=(float)vesicle->vlist.n;
    vesicle->cm[2]/=(float)vesicle->vlist.n;
    vesicle->cm[0]/=(float)n;
    vesicle->cm[1]/=(float)n;
    vesicle->cm[2]/=(float)n;
    return TS_SUCCESS;
}

ts_bool cell_ocupation(ts_vesicle *vesicle){
    ts_uint i,j,cellidx;
ts_bool cell_occupation(ts_vesicle *vesicle){
    ts_uint i,cellidx, n=vesicle->vlist->n;
    ts_double shift;
    ts_double dcell;
    shift=(ts_double) vesicle->clist.ncmax[0]/2;
    shift=(ts_double) vesicle->clist->ncmax[0]/2;
    dcell=1.0/(1.0 + vesicle->stepsize);
    ts_uint ncx, ncy,ncz;

    cell_list_cell_ocupation_clear(&vesicle->clist);
    for(i=0;i<vesicle->vlist.n;i++){
    cell_list_cell_occupation_clear(vesicle->clist);
    for(i=0;i<n;i++){
  
    cellidx=vertex_self_avoidance(vesicle, &vesicle->vlist.vertex[i]);
    vesicle->vlist.vertex[i].cell=&(vesicle->clist.cell[cellidx]);
    cellidx=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]);
    vesicle->vlist->vtx[i]->data->cell=vesicle->clist->cell[cellidx];

    cell_add_vertex(&vesicle->clist.cell[cellidx],&vesicle->vlist.vertex[i]);
    cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->vlist->vtx[i]);

  //  if(ncx > vesicle->clist.ncmax[0]) vesicle->clist.ncmax[0]=ncx;
  //  if(ncy > vesicle->clist.ncmax[1]) vesicle->clist.ncmax[1]=ncy;

 src/frame.h

@@ -0,0 +1,6 @@
#ifndef _H_FRAME
#define _H_FRAME
ts_bool centermass(ts_vesicle *vesicle);
ts_bool cell_occupation(ts_vesicle *vesicle);

#endif

 src/initial_distribution.c

@@ -8,6 +8,7 @@
#include "vertex.h"
#include "triangle.h"
#include "initial_distribution.h"
#include "energy.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);
@@ -24,6 +25,7 @@
    retval = init_triangles(vesicle);
    retval = init_triangle_neighbours(vesicle);
    retval = init_common_vertex_triangle_neighbours(vesicle);
    retval = mean_curvature_and_energy(vesicle);
 ts_fprintf(stderr,"initial_distribution finished!\n");
    return vesicle;
} 

 src/io.c

@@ -1,7 +1,9 @@
#include<general.h>
#include "general.h"
#include<stdio.h>
#include "io.h"
#include "confuse.h"
#include <confuse.h>
#include "vertex.h"
#include "bond.h"
#include<string.h>
#include<stdlib.h>
#include <sys/types.h>
@@ -194,6 +196,7 @@
            }  
        }
    }
    free(dir);
    fprintf(fh,"</Collection>\n</VTKFile>\n");
    fclose(fh);
    return TS_SUCCESS;
@@ -330,3 +333,49 @@
    return TS_SUCCESS;

}


ts_bool read_geometry_file(char *fname, ts_vesicle *vesicle){
    FILE *fh;
    ts_uint i, nvtx,nedges,ntria;
    ts_uint vtxi1,vtxi2;
    float x,y,z;
    ts_vertex_list *vlist;
    fh=fopen(fname, "r");
        if(fh==NULL){
                err("Cannot open file for reading... Nonexistant file?");
                return TS_FAIL;
        }
    ts_uint retval;
    retval=fscanf(fh,"%u %u %u",&nvtx, &nedges, &ntria);
    vesicle->vlist=init_vertex_list(nvtx);
    vlist=vesicle->vlist;
    for(i=0;i<nvtx;i++){
   //     fscanf(fh,"%F %F %F",&vlist->vtx[i]->data->x,&vlist->vtx[i]->data->y,&vlist->vtx[i]->data->z);
       retval=fscanf(fh,"%F %F %F",&x,&y,&z);
        vlist->vtx[i]->data->x=x;
        vlist->vtx[i]->data->y=y;
        vlist->vtx[i]->data->z=z;
    }
    for(i=0;i<nedges;i++){
        retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
        bond_add(vesicle->blist,vesicle->vlist->vtx[vtxi1-1],vesicle->vlist->vtx[vtxi2-1]);
    }
    //TODO: neighbours from bonds,
    //TODO: triangles from neigbours

//    Don't need to read triangles. Already have enough data
    /*
    for(i=0;i<ntria;i++){
        retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
        vtxi1=vesicle->blist->data->vertex1->idx;
        vtxi2=vesicle->blist->data->vertex1->idx;
        
    }
    */
    fclose(fh);	



    return TS_SUCCESS;
}

 src/tape

@@ -21,3 +21,14 @@

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

#what type of multiprocessing? (*none, smp, cluster, distributed, cuda, auto)
#currently only none makes sense.
multiprocessing=none
#how many cores are allowed to process in SMP?
smp_cores=2
#how many nodes in cluster?
cluster_nodes=50
#max number of processes in distributed (voluntary) environment
distributed_processes=50
#cuda options???

 src/vertex.c

@@ -214,3 +214,37 @@
    return TS_SUCCESS;
}


/* ****************************************************************** */
/* ***** New vertex copy operations. Inherently they are slow.  ***** */
/* ****************************************************************** */

ts_vertex *vtx_copy(ts_vertex *ovtx){
    ts_vertex *cvtx=(ts_vertex *)malloc(sizeof(ts_vertex));
    memcpy((void *)cvtx,(void *)ovtx,sizeof(ts_vertex));
    cvtx->data=(ts_vertex_data *)malloc(sizeof(ts_vertex_data));
    memcpy((void *)cvtx->data,(void *)ovtx->data,sizeof(ts_vertex_data));
    cvtx->data->neigh=NULL;
    cvtx->data->neigh_no=0;
    cvtx->data->tristar_no=0;
    cvtx->data->bond_no=0;
    cvtx->data->tristar=NULL;
    cvtx->data->bond=NULL;
    cvtx->data->cell=NULL; 
    return cvtx;
}
//TODO: needs to be done
ts_vertex **vtx_neigh_copy(ts_vertex_list *vlist,ts_vertex *ovtx){
        return NULL;
}

ts_vertex_list *vertex_list_copy(ts_vertex_list *ovlist){
    ts_uint i;
    ts_vertex_list *vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
    memcpy((void *)vlist, (void *)ovlist, sizeof(ts_vertex_list));
    ts_vertex **vtx=(ts_vertex **)malloc(vlist->n*sizeof(ts_vertex *));
    for(i=0;i<vlist->n;i++){
        vtx[i]=vtx_copy(ovlist->vtx[i]);
    }
    return vlist;
}

 src/vertex.h

@@ -26,4 +26,9 @@
ts_bool vtx_set_global_values(ts_vesicle *vesicle);
inline ts_double vtx_direct(ts_vertex *vtx1, ts_vertex *vtx2, ts_vertex *vtx3);
inline ts_bool vertex_add_tristar(ts_vertex *vtx, ts_triangle *tristarmem);

ts_vertex *vtx_copy(ts_vertex *ovtx);
ts_vertex **vtx_neigh_copy(ts_vertex_list *vlist,ts_vertex *ovtx);
ts_vertex_list *vertex_list_copy(ts_vertex_list *ovlist);

#endif