trisurf-ng.git - Gitblit

To be sure a fix in comparing the floating point numbers

Samo Penic

2014-06-15 7d62fc765adc11bbcb7db0716ac5dc19145dca80

commit \| author \| age
460c2a	1	#include<stdlib.h>
SP	2	#include<stdio.h>
	3	#include<string.h>
	4	#include<math.h>
	5	#include "general.h"
	6	#include "constvol.h"
	7	#include "triangle.h"
	8	#include "energy.h"
	9	#include "vertex.h"
	10	#include "cell.h"
	11
fbcbdf	12	ts_bool constvolume(ts_vesicle vesicle, ts_vertex vtx_avoid, ts_double Vol, ts_double retEnergy, ts_vertex vtx_moved_retval, ts_vertex *vtx_backup){
SP	13	ts_vertex *vtx_moved;
460c2a	14	ts_uint vtxind,i,j;
41637a	15	ts_uint Ntries=100;
460c2a	16	ts_vertex *backupvtx;
fe5069	17	ts_double Rv, dh, dvol, volFirst, voldiff, oenergy,delta_energy;
460c2a	18	backupvtx=(ts_vertex *)calloc(sizeof(ts_vertex),10);
fbcbdf	19	ts_double l0 = (1.0 + sqrt(vesicle->dmax))/2.0; //make this a global constant if necessary
460c2a	20	for(i=0;i<Ntries;i++){
SP	21	vtxind=rand() % vesicle->vlist->n;
	22	vtx_moved=vesicle->vlist->vtx[vtxind];
41637a	23	/* chosen vertex must not be a nearest neighbour. WASTODO: probably must.
SP	24	* DONE: solved below, when using different algorithm.
fe5069	25	* extend search in case of bondflip */
41637a	26	/* if(vtx_moved==vtx_avoid) continue;
460c2a	27	for(j=0;j<vtx_moved->neigh_no;j++){
SP	28	if(vtx_moved->neigh[j]==vtx_avoid) continue;
	29	}
41637a	30	*/
SP	31	/* different check of nearest neighbour (and second nearest neighbour) check.
	32	* Checking the distance between vertex and vertex to be moved to assure
	33	* constant volume. Solves upper todo problem. */
	34
	35	if(vtx_distance_sq(vtx_moved,vtx_avoid)<16.0*vesicle->dmax){
	36	continue;
	37	}
	38
460c2a	39	memcpy((void )&backupvtx[0],(void )vtx_moved,sizeof(ts_vertex));
SP	40
fe5069	41	//move vertex in specified direction. first try, test move!
460c2a	42	Rv=sqrt(pow(vtx_moved->x,2)+pow(vtx_moved->y,2)+pow(vtx_moved->z,2));
fbcbdf	43	dh=2.0Vol/(sqrt(3.0)l0*l0);
SP	44	vtx_moved->x=vtx_moved->x*(1.0-dh/Rv);
	45	vtx_moved->y=vtx_moved->y*(1.0-dh/Rv);
	46	vtx_moved->z=vtx_moved->z*(1.0-dh/Rv);
460c2a	47
41637a	48
SP	49	//SKIPPING FIRST CHECK of CONSTRAINTS. This is not a final move.
460c2a	50	//check for constraints
41637a	51	/* if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
460c2a	52	vtx_moved=memcpy((void )vtx_moved,(void )&backupvtx[0],sizeof(ts_vertex));
41637a	53	// continue;
SP	54	break;
	55	} */
460c2a	56	// All checks OK!
SP	57
	58	dvol=0.0;
	59	for(j=0;j<vtx_moved->tristar_no;j++){
	60	dvol-=vtx_moved->tristar[j]->volume;
fe5069	61	}
SP	62	volFirst=dvol;
	63	for(j=0;j<vtx_moved->tristar_no;j++){
460c2a	64	triangle_normal_vector(vtx_moved->tristar[j]);
SP	65	dvol+=vtx_moved->tristar[j]->volume;
	66	}
fe5069	67	//TODO: here there is a bug. Don't know where, but preliminary success can
SP	68	//happen sometimes. And when I do this checks, constant value is not achieved
	69	//anymore
	70	/* voldiff=dvol-Vol;
460c2a	71	if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
SP	72	//calculate energy, return change in energy...
	73	oenergy=vtx_moved->energy;
	74	energy_vertex(vtx_moved);
	75	delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
	76	//the same is done for neighbouring vertices
fe5069	77	for(j=0;j<vtx_moved->neigh_no;j++){
SP	78	oenergy=vtx_moved->neigh[j]->energy;
	79	energy_vertex(vtx_moved->neigh[j]);
	80	delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy);
460c2a	81	}
SP	82	*retEnergy=delta_energy;
fbcbdf	83	*vtx_backup=backupvtx;
SP	84	*vtx_moved_retval=vtx_moved;
b2fa8c	85	fprintf(stderr, "Preliminary success\n");
460c2a	86	return TS_SUCCESS;
fe5069	87	} */
SP	88
	89	// fprintf(stderr, "Step 2 success\n");
460c2a	90	//do it again ;)
SP	91	dh=Vol*dh/dvol;
	92	vtx_moved=memcpy((void )vtx_moved,(void )&backupvtx[0],sizeof(ts_vertex));
	93	vtx_moved->x=vtx_moved->x*(1-dh/Rv);
	94	vtx_moved->y=vtx_moved->y*(1-dh/Rv);
	95	vtx_moved->z=vtx_moved->z*(1-dh/Rv);
	96	//check for constraints
	97	if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
	98	vtx_moved=memcpy((void )vtx_moved,(void )&backupvtx[0],sizeof(ts_vertex));
fe5069	99	//also, restore normals
SP	100	for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]);
41637a	101	// continue;
SP	102	break;
460c2a	103	}
SP	104
fe5069	105	dvol=volFirst;
b2fa8c	106	for(j=0;j<vtx_moved->tristar_no;j++){
SP	107	triangle_normal_vector(vtx_moved->tristar[j]);
	108	dvol+=vtx_moved->tristar[j]->volume;
	109	}
460c2a	110	voldiff=dvol-Vol;
SP	111	if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
	112	//calculate energy, return change in energy...
fe5069	113	// fprintf(stderr, "Constvol success! %e\n",voldiff);
3de289	114	for(j=0;j<vtx_moved->neigh_no;j++){
SP	115	memcpy((void )&backupvtx[j+1],(void )vtx_moved->neigh[j],sizeof(ts_vertex));
	116	}
	117
460c2a	118	oenergy=vtx_moved->energy;
SP	119	energy_vertex(vtx_moved);
	120	delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
	121	//the same is done for neighbouring vertices
fe5069	122	for(j=0;j<vtx_moved->neigh_no;j++){
SP	123	oenergy=vtx_moved->neigh[j]->energy;
	124	energy_vertex(vtx_moved->neigh[j]);
	125	delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy);
460c2a	126	}
SP	127	*retEnergy=delta_energy;
fbcbdf	128	*vtx_backup=backupvtx;
SP	129	*vtx_moved_retval=vtx_moved;
460c2a	130	return TS_SUCCESS;
SP	131	}
	132	}
	133	free(backupvtx);
	134	return TS_FAIL;
	135	}
	136
	137
	138	ts_bool constvolConstraintCheck(ts_vesicle vesicle, ts_vertex vtx){
	139	ts_uint i;
	140	ts_double dist;
	141	ts_uint cellidx;
	142	//distance with neighbours check
	143	for(i=0;i<vtx->neigh_no;i++){
	144	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
	145	if(dist<1.0 \|\| dist>vesicle->dmax) {
	146	return TS_FAIL;
	147	}
	148	}
	149	// Distance with grafted poly-vertex check:
	150	if(vtx->grafted_poly!=NULL){
	151	dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
	152	if(dist<1.0 \|\| dist>vesicle->dmax) {
	153	return TS_FAIL;
	154	}
	155	}
	156
	157	// Nucleus penetration check:
	158	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	159	return TS_FAIL;
	160	}
	161
	162	//self avoidance check with distant vertices
	163	cellidx=vertex_self_avoidance(vesicle, vtx);
	164	//check occupation number
	165	return cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	166	}
	167
	168
	169
	170	ts_bool constvolumerestore(ts_vertex vtx_moved,ts_vertex vtx_backup){
	171	ts_uint j;
fbcbdf	172	memcpy((void )vtx_moved,(void )&vtx_backup[0],sizeof(ts_vertex));
fe5069	173	for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]);
3de289	174	for(j=0;j<vtx_moved->neigh_no;j++){
SP	175	// memcpy((void )vtx_moved->neigh[j],(void )&vtx_backup[j+1],sizeof(ts_vertex));
	176	energy_vertex(vtx_moved->neigh[j]);
	177	}
fe5069	178
fbcbdf	179	free(vtx_backup);
460c2a	180	return TS_SUCCESS;
SP	181	}
	182
fbcbdf	183	ts_bool constvolumeaccept(ts_vesicle vesicle,ts_vertex vtx_moved, ts_vertex *vtx_backup){
SP	184	ts_bool retval;
	185	ts_uint cellidx=vertex_self_avoidance(vesicle, vtx_moved);
	186	if(vtx_moved->cell!=vesicle->clist->cell[cellidx]){
	187	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx_moved);
	188	if(retval==TS_SUCCESS) cell_remove_vertex(vtx_backup[0].cell,vtx_moved);
	189
	190	}
	191	free(vtx_backup);
460c2a	192
SP	193	return TS_SUCCESS;
	194	}