trisurf-ng.git - Gitblit

Samo Penic

2014-04-30 fe5069cc3872b513f4715b0cfd303175ae80c468

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