trisurf-ng.git - Gitblit

mihaf

2014-05-13 0f9d229fb830435d92cfc73c85aa7771a420286c

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