From 523bf18206f550a315c6c17e5a0a253381b0f8bf Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@fe.uni-lj.si> Date: Thu, 07 Jun 2012 11:16:16 +0000 Subject: [PATCH] Spherical harmonics. Almost everyhing is done. Missing triangle area calculation when vertex is moved or bond is flipped. Also missing volume calculation on vertex move or bondflip. Calculation of co coefficient is not done completely yet. Problems are in numbering the coefficients. Newly added data structure ts_spharm is referenced from ts_vesicle. Missing function for initialization and freeing the memory of that datastructure -- but that memory is already used by some functions. --- src/vertexmove.c | 70 +++++++++++++++++++---------------- 1 files changed, 38 insertions(+), 32 deletions(-) diff --git a/src/vertexmove.c b/src/vertexmove.c index 243463f..c62aeda 100644 --- a/src/vertexmove.c +++ b/src/vertexmove.c @@ -16,55 +16,59 @@ *rn){ ts_uint i; ts_double dist; - ts_vertex *tvtx=(ts_vertex *)malloc(sizeof(ts_vertex)); - tvtx->data=init_vertex_data(); ts_bool retval; ts_uint cellidx; ts_double xold,yold,zold; ts_double delta_energy,oenergy; ts_vertex *ovtx; + ts_vertex *tvtx=(ts_vertex *)calloc(1,sizeof(ts_vertex)); //randomly we move the temporary vertex - tvtx->data->x=vtx->data->x+vesicle->stepsize*(2.0*rn[0]-1.0); - tvtx->data->y=vtx->data->y+vesicle->stepsize*(2.0*rn[1]-1.0); - tvtx->data->z=vtx->data->z+vesicle->stepsize*(2.0*rn[2]-1.0); + tvtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0); + tvtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0); + tvtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0); //check we if some length to neighbours are too much - for(i=0;i<vtx->data->neigh_no;i++){ - dist=vtx_distance_sq(tvtx,vtx->data->neigh[i]); - if(dist<1.0 || dist>vesicle->dmax) return TS_FAIL; + for(i=0;i<vtx->neigh_no;i++){ + dist=vtx_distance_sq(tvtx,vtx->neigh[i]); + if(dist<1.0 || dist>vesicle->dmax) { + vtx_free(tvtx); +// fprintf(stderr,"Fail 1, dist=%f, vesicle->dmax=%f\n", dist, vesicle->dmax); + return TS_FAIL; + } } -fprintf(stderr,"Was here!\n"); //self avoidance check with distant vertices cellidx=vertex_self_avoidance(vesicle, tvtx); //check occupation number retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx,tvtx); if(retval==TS_FAIL){ + vtx_free(tvtx); +// fprintf(stderr,"Fail 2\n"); return TS_FAIL; } //if all the tests are successful, then we update the vertex position - xold=vtx->data->x; - yold=vtx->data->y; - zold=vtx->data->z; + xold=vtx->x; + yold=vtx->y; + zold=vtx->z; ovtx=malloc(sizeof(ts_vertex)); vtx_copy(ovtx,vtx); - vtx->data->x=tvtx->data->x; - vtx->data->y=tvtx->data->y; - vtx->data->z=tvtx->data->z; + vtx->x=tvtx->x; + vtx->y=tvtx->y; + vtx->z=tvtx->z; delta_energy=0; //update the normals of triangles that share bead i. - for(i=0;i<vtx->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]); + for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); //energy and curvature energy_vertex(vtx); - delta_energy=vtx->data->xk*(vtx->data->energy - ovtx->data->energy); + delta_energy=vtx->xk*(vtx->energy - ovtx->energy); //the same is done for neighbouring vertices - for(i=0;i<vtx->data->neigh_no;i++){ - oenergy=vtx->data->neigh[i]->data->energy; - energy_vertex(vtx->data->neigh[i]); - delta_energy+=vtx->data->neigh[i]->data->xk*(vtx->data->neigh[i]->data->energy-oenergy); + for(i=0;i<vtx->neigh_no;i++){ + oenergy=vtx->neigh[i]->energy; + energy_vertex(vtx->neigh[i]); + delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy); } - fprintf(stderr, "DE=%f\n",delta_energy); +// fprintf(stderr, "DE=%f\n",delta_energy); //MONTE CARLOOOOOOOO if(delta_energy>=0){ #ifdef TS_DOUBLE_DOUBLE @@ -78,28 +82,30 @@ #endif { //not accepted, reverting changes - vtx->data->x=xold; - vtx->data->y=yold; - vtx->data->z=zold; + vtx->x=xold; + vtx->y=yold; + vtx->z=zold; //update the normals of triangles that share bead i. - for(i=0;i<vtx->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]); + for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); //energy and curvature energy_vertex(vtx); //the same is done for neighbouring vertices - for(i=0;i<vtx->data->neigh_no;i++) energy_vertex(vtx->data->neigh[i]); - free(ovtx->data->bond_length); - free(ovtx->data->bond_length_dual); + for(i=0;i<vtx->neigh_no;i++) energy_vertex(vtx->neigh[i]); + free(ovtx->bond_length); + free(ovtx->bond_length_dual); free(ovtx); + vtx_free(tvtx); return TS_FAIL; } } //END MONTE CARLOOOOOOO //TODO: change cell occupation if necessary! - - free(ovtx->data->bond_length); - free(ovtx->data->bond_length_dual); +// fprintf(stderr,"Success!!\n"); + free(ovtx->bond_length); + free(ovtx->bond_length_dual); free(ovtx); + vtx_free(tvtx); return TS_SUCCESS; } -- Gitblit v1.9.3