trisurf-ng.git - Gitblit

Samo Penic

2014-04-28 460c2a326b332e1b3e03621970a1b48dfb872050

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
	12	ts_bool constvolume(ts_vesicle vesicle, ts_vertex vtx_avoid, ts_double Vol, ts_double retEnergy, ts_vertex vtx_moved, ts_vertex *vtx_backup){
	13
	14	ts_uint vtxind,i,j;
	15	ts_uint Ntries=20;
	16	ts_vertex *backupvtx;
	17	ts_double Rv, dh, dvol, voldiff, oenergy,delta_energy;
	18
	19	backupvtx=(ts_vertex *)calloc(sizeof(ts_vertex),10);
	20
	21	for(i=0;i<Ntries;i++){
	22	vtxind=rand() % vesicle->vlist->n;
	23	vtx_moved=vesicle->vlist->vtx[vtxind];
	24	if(vtx_moved==vtx_avoid) continue;
	25
	26	for(j=0;j<vtx_moved->neigh_no;j++){
	27	if(vtx_moved->neigh[j]==vtx_avoid) continue;
	28	}
	29
	30	memcpy((void )&backupvtx[0],(void )vtx_moved,sizeof(ts_vertex));
	31	//move vertex in specified direction. first try, test move!
	32
	33	Rv=sqrt(pow(vtx_moved->x,2)+pow(vtx_moved->y,2)+pow(vtx_moved->z,2));
	34	dh=2Rvvesicle->dmax/sqrt(3);
	35	vtx_moved->x=vtx_moved->x*(1-dh/Rv);
	36	vtx_moved->y=vtx_moved->y*(1-dh/Rv);
	37	vtx_moved->z=vtx_moved->z*(1-dh/Rv);
	38
	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
	45	// All checks OK!
	46
	47	// doing second and final move.
	48	for(j=0;j<vtx_moved->neigh_no;j++){
	49	memcpy((void )&backupvtx[j+1],(void )vtx_moved->neigh[j],sizeof(ts_vertex));
	50	}
	51	dvol=0.0;
	52	for(j=0;j<vtx_moved->tristar_no;j++){
	53	dvol-=vtx_moved->tristar[j]->volume;
	54	triangle_normal_vector(vtx_moved->tristar[j]);
	55	dvol+=vtx_moved->tristar[j]->volume;
	56	}
	57
	58	voldiff=dvol-Vol;
	59
	60	if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
	61	//calculate energy, return change in energy...
	62	oenergy=vtx_moved->energy;
	63	energy_vertex(vtx_moved);
	64	delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
	65	//the same is done for neighbouring vertices
	66	for(i=0;i<vtx_moved->neigh_no;i++){
	67	oenergy=vtx_moved->neigh[i]->energy;
	68	energy_vertex(vtx_moved->neigh[i]);
	69	delta_energy+=vtx_moved->neigh[i]->xk*(vtx_moved->neigh[i]->energy-oenergy);
	70	}
	71	*retEnergy=delta_energy;
	72	vtx_backup=backupvtx;
	73	return TS_SUCCESS;
	74	}
	75	//do it again ;)
	76	dh=Vol*dh/dvol;
	77	vtx_moved=memcpy((void )vtx_moved,(void )&backupvtx[0],sizeof(ts_vertex));
	78	vtx_moved->x=vtx_moved->x*(1-dh/Rv);
	79	vtx_moved->y=vtx_moved->y*(1-dh/Rv);
	80	vtx_moved->z=vtx_moved->z*(1-dh/Rv);
	81	//check for constraints
	82	if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
	83	for(j=0;j<vtx_moved->neigh_no;j++){
	84	memcpy((void )vtx_moved->neigh[j],(void )&backupvtx[j+1],sizeof(ts_vertex));
	85	}
	86	vtx_moved=memcpy((void )vtx_moved,(void )&backupvtx[0],sizeof(ts_vertex));
	87	continue;
	88	}
	89
	90	voldiff=dvol-Vol;
	91	if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
	92	//calculate energy, return change in energy...
	93	oenergy=vtx_moved->energy;
	94	energy_vertex(vtx_moved);
	95	delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
	96	//the same is done for neighbouring vertices
	97	for(i=0;i<vtx_moved->neigh_no;i++){
	98	oenergy=vtx_moved->neigh[i]->energy;
	99	energy_vertex(vtx_moved->neigh[i]);
	100	delta_energy+=vtx_moved->neigh[i]->xk*(vtx_moved->neigh[i]->energy-oenergy);
	101	}
	102	*retEnergy=delta_energy;
	103	vtx_backup=backupvtx;
	104	return TS_SUCCESS;
	105	}
	106
	107
	108	}
	109	free(backupvtx);
	110	return TS_FAIL;
	111	}
	112
	113
	114	ts_bool constvolConstraintCheck(ts_vesicle vesicle, ts_vertex vtx){
	115	ts_uint i;
	116	ts_double dist;
	117	ts_uint cellidx;
	118	//distance with neighbours check
	119	for(i=0;i<vtx->neigh_no;i++){
	120	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
	121	if(dist<1.0 \|\| dist>vesicle->dmax) {
	122	return TS_FAIL;
	123	}
	124	}
	125	// Distance with grafted poly-vertex check:
	126	if(vtx->grafted_poly!=NULL){
	127	dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
	128	if(dist<1.0 \|\| dist>vesicle->dmax) {
	129	return TS_FAIL;
	130	}
	131	}
	132
	133	// Nucleus penetration check:
	134	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	135	return TS_FAIL;
	136	}
	137
	138	//self avoidance check with distant vertices
	139	cellidx=vertex_self_avoidance(vesicle, vtx);
	140	//check occupation number
	141	return cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	142	}
	143
	144
	145
	146	ts_bool constvolumerestore(ts_vertex vtx_moved,ts_vertex vtx_backup){
	147	ts_uint j;
	148	for(j=0;j<vtx_moved->neigh_no;j++){
	149	memcpy((void )vtx_moved->neigh[j],(void )&vtx_backup[j+1],sizeof(ts_vertex));
	150	}
	151	vtx_moved=memcpy((void )vtx_moved,(void )&vtx_backup[0],sizeof(ts_vertex));
	152
	153	return TS_SUCCESS;
	154	}
	155
	156	ts_bool constvolumeaccept(ts_vertex vtx_moved, ts_vertex vtx_backup){
	157
	158
	159	return TS_SUCCESS;
	160	}