trisurf-ng.git - Gitblit

Samo Penic

2014-03-05 719c9febac2eaff9483fda487b57684afbb59bb2

commit \| author \| age
f74313	1	#include "general.h"
d7639a	2	#include<stdio.h>
SP	3	#include "io.h"
f74313	4	#include <confuse.h>
SP	5	#include "vertex.h"
	6	#include "bond.h"
d7639a	7	#include<string.h>
a6b1b5	8	#include<stdlib.h>
d7639a	9	#include <sys/types.h>
SP	10	#include <dirent.h>
b14a8d	11	#include "initial_distribution.h"
a2db52	12	#include "poly.h"
d7639a	13
e9eab4	14
SP	15
	16	/ DUMP STATE TO DISK DRIVE /
	17
	18	ts_bool dump_state(ts_vesicle *vesicle){
	19
	20	/* save current state with wrong pointers. Will fix that later */
	21	ts_uint i,j,k;
	22	ts_ulong off;
	23	FILE *fh=fopen("dump.bin","wb");
	24
	25	/* dump vesicle */
	26	fwrite(vesicle, sizeof(vesicle),1,fh);
	27	/* dump vertex list */
	28	fwrite(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
	29	/* dump bond list */
	30	fwrite(vesicle->blist, sizeof(ts_bond_list),1,fh);
	31	/* dump triangle list */
	32	fwrite(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
	33	/* dump cell list */
	34	fwrite(vesicle->clist, sizeof(ts_cell_list),1,fh);
	35	/* dump poly list */
	36	fwrite(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
	37	/* level 1 complete */
	38
	39	/dump vertices/
	40	for(i=0;i<vesicle->vlist->n;i++){
	41	fwrite(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
	42	/* dump pointer offsets for:
	43	neigh
	44	bond
	45	tria
	46	cell is ignored
	47	*/
	48	for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
	49	off=(ts_ulong)(vesicle->vlist->vtx[i]->neigh[j]-vesicle->vlist->vtx[0]);
	50	fwrite(&off,sizeof(ts_ulong),1,fh);
	51	}
	52	for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
	53	off=(ts_ulong)(vesicle->vlist->vtx[i]->bond[j]-vesicle->blist->bond[0]);
	54	fwrite(&off,sizeof(ts_ulong),1,fh);
	55	}
	56	for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
	57	off=(ts_ulong)(vesicle->vlist->vtx[i]->tristar[j]-vesicle->tlist->tria[0]);
	58	fwrite(&off,sizeof(ts_ulong),1,fh);
	59	}
	60	// off=(ts_ulong)(vesicle->vlist->vtx[i]->cell-vesicle->clist->cell[0]);
	61	// fwrite(&off,sizeof(ts_ulong),1,fh);
	62
	63	}
	64
	65	/dump bonds/
	66	for(i=0;i<vesicle->blist->n;i++){
	67	fwrite(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
	68	/* dump pointer offsets for vtx1 and vtx2 */
	69	off=(ts_ulong)(vesicle->blist->bond[i]->vtx1-vesicle->vlist->vtx[0]);
	70	fwrite(&off,sizeof(ts_ulong),1,fh);
	71	off=(ts_ulong)(vesicle->blist->bond[i]->vtx2-vesicle->vlist->vtx[0]);
	72	fwrite(&off,sizeof(ts_ulong),1,fh);
	73	}
	74
	75	/dump triangles/
	76	for(i=0;i<vesicle->tlist->n;i++){
	77	fwrite(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
	78	/* dump pointer offsets for vertex */
	79	off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[0]-vesicle->vlist->vtx[0]);
	80	fwrite(&off,sizeof(ts_ulong),1,fh);
	81	off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[1]-vesicle->vlist->vtx[0]);
	82	fwrite(&off,sizeof(ts_ulong),1,fh);
	83	off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[2]-vesicle->vlist->vtx[0]);
	84	fwrite(&off,sizeof(ts_ulong),1,fh);
	85	/* dump pointer offsets for neigh */
	86	for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
	87	off=(ts_ulong)(vesicle->tlist->tria[i]->neigh[j]-vesicle->tlist->tria[0]);
	88	fwrite(&off,sizeof(ts_ulong),1,fh);
	89	}
	90	}
	91
	92	/dump cells and its contents is not dumped but recalculated on restore/
	93
	94	/* for(i=0;i<vesicle->clist->cellno;i++){
	95	fwrite(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
	96	for(j=0;j<vesicle->clist->cell[i]->nvertex;j++){
	97	off=(ts_ulong)(vesicle->clist->cell[i]->vertex[j]-vesicle->vlist->vtx[0]);
	98	fwrite(&off,sizeof(ts_ulong),1,fh);
	99	}
	100	}
	101	*/
	102
	103
	104	/dump polymeres /
	105	for(i=0;i<vesicle->poly_list->n;i++){
	106	fwrite(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
	107	}
	108
	109	/* dump poly vertex(monomer) list*/
	110	for(i=0;i<vesicle->poly_list->n;i++){
	111	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
	112	fwrite(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
	113	/* dump offset for neigh and bond */
	114	for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
	115	off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]-vesicle->poly_list->poly[i]->vlist->vtx[0]);
	116	fwrite(&off,sizeof(ts_ulong),1,fh);
	117	}
	118	for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
	119	off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]-vesicle->poly_list->poly[i]->blist->bond[0]);
	120	fwrite(&off,sizeof(ts_ulong),1,fh);
	121	}
	122	}
	123	}
	124	/* dump poly bonds between monomers list*/
	125	for(i=0;i<vesicle->poly_list->n;i++){
	126	for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
	127	fwrite(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
	128	/* dump vtx1 and vtx2 offsets */
	129	off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx1-vesicle->poly_list->poly[i]->vlist->vtx[0]);
	130	fwrite(&off,sizeof(ts_ulong),1,fh);
	131	off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]);
	132	fwrite(&off,sizeof(ts_ulong),1,fh);
	133	}
	134	}
	135
	136	/* pointer offsets for fixing the restored pointers */
	137	/* need pointers for
	138	vlist->vtx
	139	blist->bond
	140	tlist->tria
	141	clist->cell
	142	poly_list->poly
	143	and for each poly:
	144	poly_list->poly->vtx
	145	poly_list->poly->bond
	146	*/
	147
	148	fwrite(vesicle->vlist->vtx, sizeof(ts_vertex *),1,fh);
	149	fwrite(vesicle->blist->bond, sizeof(ts_bond *),1,fh);
	150	fwrite(vesicle->tlist->tria, sizeof(ts_triangle *),1,fh);
	151	fwrite(vesicle->clist->cell, sizeof(ts_cell *),1,fh);
	152	fwrite(vesicle->poly_list->poly, sizeof(ts_poly *),1,fh);
	153	for(i=0;i<vesicle->poly_list->n;i++){
	154	fwrite(vesicle->poly_list->poly[i]->vlist->vtx, sizeof(ts_vertex *),1,fh);
	155	fwrite(vesicle->poly_list->poly[i]->blist->bond, sizeof(ts_bond *),1,fh);
	156	}
	157	fclose(fh);
	158	return TS_SUCCESS;
	159	}
	160
	161
	162	/ RESTORE DUMP FROM DISK /
	163	ts_vesicle *restore_state(){
	164	ts_uint i,j,k;
	165	FILE *fh=fopen("dump.bin","rb");
	166	ts_uint retval;
	167	ts_ulong off;
	168
	169	/* we restore all the data from the dump */
	170	/* restore vesicle */
	171	ts_vesicle vesicle=(ts_vesicle )calloc(1,sizeof(ts_vesicle));
	172	retval=fread(vesicle, sizeof(vesicle),1,fh);
	173	/* restore vertex list */
	174	vesicle->vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
	175	retval=fread(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
	176	/* restore bond list */
	177	vesicle->blist=(ts_bond_list *)malloc(sizeof(ts_bond_list));
	178	retval=fread(vesicle->blist, sizeof(ts_bond_list),1,fh);
	179	/* restore triangle list */
	180	vesicle->tlist=(ts_triangle_list *)malloc(sizeof(ts_triangle_list));
	181	retval=fread(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
	182	/* restore cell list */
	183	vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
	184	retval=fread(vesicle->clist, sizeof(ts_cell_list),1,fh);
	185	/* restore poly list */
	186	vesicle->poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
	187	retval=fread(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
	188	/* level 1 complete */
	189
	190	/restore vertices/
	191	vesicle->vlist->vtx=(ts_vertex *)calloc(vesicle->vlist->n,sizeof(ts_vertex ));
	192	for(i=0;i<vesicle->vlist->n;i++){
	193	vesicle->vlist->vtx[i]=(ts_vertex *)malloc(sizeof(ts_vertex));
	194	retval=fread(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
	195	/restore neigh, bond, tristar. Ignoring cell /
	196	vesicle->vlist->vtx[i]->neigh=(ts_vertex *)calloc(vesicle->vlist->vtx[i]->neigh_no, sizeof(ts_vertex ));
	197	for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
	198	retval=fread(&off,sizeof(ts_ulong),1,fh);
	199	vesicle->vlist->vtx[i]->neigh[j]=vesicle->vlist->vtx[0]+off;
	200	}
	201	vesicle->vlist->vtx[i]->bond=(ts_bond *)calloc(vesicle->vlist->vtx[i]->bond_no, sizeof(ts_bond ));
	202	for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
	203	retval=fread(&off,sizeof(ts_ulong),1,fh);
	204	vesicle->vlist->vtx[i]->bond[j]=vesicle->blist->bond[0]+off;
	205	}
	206
	207	vesicle->vlist->vtx[i]->tristar=(ts_triangle *)calloc(vesicle->vlist->vtx[i]->tristar_no, sizeof(ts_triangle ));
	208	for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
	209	retval=fread(&off,sizeof(ts_ulong),1,fh);
	210	vesicle->vlist->vtx[i]->tristar[j]=vesicle->tlist->tria[0]+off;
	211	}
	212
	213	}
	214
	215	/restore bonds/
	216	vesicle->blist->bond=(ts_bond *)calloc(vesicle->blist->n,sizeof(ts_bond ));
	217	for(i=0;i<vesicle->blist->n;i++){
	218	vesicle->blist->bond[i]=(ts_bond *)malloc(sizeof(ts_bond));
	219	retval=fread(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
	220	/* restore vtx1 and vtx2 */
	221	retval=fread(&off,sizeof(ts_ulong),1,fh);
	222	vesicle->blist->bond[i]->vtx1=vesicle->vlist->vtx[0]+off;
	223	retval=fread(&off,sizeof(ts_ulong),1,fh);
	224	vesicle->blist->bond[i]->vtx2=vesicle->vlist->vtx[0]+off;
	225	}
	226
	227	/restore triangles/
	228	vesicle->tlist->tria=(ts_triangle *)calloc(vesicle->tlist->n,sizeof(ts_triangle ));
	229	for(i=0;i<vesicle->tlist->n;i++){
	230	vesicle->tlist->tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle));
	231	retval=fread(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
	232	/* restore pointers for vertices */
	233	retval=fread(&off,sizeof(ts_ulong),1,fh);
	234	vesicle->tlist->tria[i]->vertex[0]=vesicle->vlist->vtx[0]+off;
	235	retval=fread(&off,sizeof(ts_ulong),1,fh);
	236	vesicle->tlist->tria[i]->vertex[1]=vesicle->vlist->vtx[0]+off;
	237	retval=fread(&off,sizeof(ts_ulong),1,fh);
	238	vesicle->tlist->tria[i]->vertex[2]=vesicle->vlist->vtx[0]+off;
	239	/* restore pointers for neigh */
	240	for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
	241	retval=fread(&off,sizeof(ts_ulong),1,fh);
	242	vesicle->tlist->tria[i]->neigh[j]=vesicle->tlist->tria[0]+off;
	243	}
	244
	245	}
	246
	247	/restore cells /
	248	/TODO: do we need to recalculate cells here? /
	249	/* vesicle->clist->cell=(ts_cell *)malloc(vesicle->clist->cellnosizeof(ts_cell *));
	250	for(i=0;i<vesicle->clist->cellno;i++){
	251	vesicle->clist->cell[i]=(ts_cell *)malloc(sizeof(ts_cell));
	252	retval=fread(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
	253	}
	254	*/
	255	/restore polymeres /
	256	vesicle->poly_list->poly = (ts_poly *)calloc(vesicle->poly_list->n,sizeof(ts_poly ));
	257	for(i=0;i<vesicle->poly_list->n;i++){
	258	vesicle->poly_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly));
	259	retval=fread(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
	260	}
	261
	262	/* restore poly vertex(monomer) list*/
	263	for(i=0;i<vesicle->poly_list->n;i++){
	264	vesicle->poly_list->poly[i]->vlist->vtx=(ts_vertex *)calloc(vesicle->poly_list->poly[i]->vlist->n,sizeof(ts_vertex ));
	265	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
	266	vesicle->poly_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex));
	267	retval=fread(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
	268
	269	/* restore neigh and bonds */
	270	vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh=(ts_vertex *)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no, sizeof(ts_vertex ));
	271	for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
	272	retval=fread(&off,sizeof(ts_ulong),1,fh);
	273	vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
	274	}
	275	vesicle->poly_list->poly[i]->vlist->vtx[j]->bond=(ts_bond *)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no, sizeof(ts_bond ));
	276	for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
	277	retval=fread(&off,sizeof(ts_ulong),1,fh);
	278	vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->poly_list->poly[i]->blist->bond[0]+off;
	279	}
	280
	281	}
	282	}
	283	/* restore poly bonds between monomers list*/
	284	for(i=0;i<vesicle->poly_list->n;i++){
	285	vesicle->poly_list->poly[i]->blist->bond=(ts_bond *)calloc(vesicle->poly_list->poly[i]->blist->n,sizeof(ts_bond ));
	286	for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
	287	vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
	288	retval=fread(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
	289	/* restore vtx1 and vtx2 */
	290	retval=fread(&off,sizeof(ts_ulong),1,fh);
	291	vesicle->poly_list->poly[i]->blist->bond[j]->vtx1=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
	292	retval=fread(&off,sizeof(ts_ulong),1,fh);
	293	vesicle->poly_list->poly[i]->blist->bond[j]->vtx2=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
	294	}
	295	}
	296
	297
	298	/* recalculate pointers using saved information on pointer addresses */
	299	/* following pointers need to be recalculated:
	300	all vtx->neigh
	301	all vtx->tristar
	302	all vtx->bond
	303	vtx->cell
	304	vtx->grafted_poly
	305	bond->vtx1
	306	bond->vtx2
	307	tria->vertex[*]
	308	tria->neigh
	309	cell->vertex
	310	poly->grafted_vtx
	311	poly->vlist->vtx->neigh
	312	poly->vlist->vtx->bond
	313	poly->vlist->vtx->cell
	314	*/
	315
	316	/* read pointers from disk */
	317	/* fwrite(vesicle->vlist->vtx, sizeof(ts_vertex *),1,fh);
	318	fwrite(vesicle->blist->bond, sizeof(ts_bond *),1,fh);
	319	fwrite(vesicle->tlist->tria, sizeof(ts_triangle *),1,fh);
	320	fwrite(vesicle->clist->cell, sizeof(ts_cell *),1,fh);
	321	fwrite(vesicle->poly_list->poly, sizeof(ts_poly *),1,fh);
	322	for(i=0;i<vesicle->poly_list->n;i++){
	323	fwrite(vesicle->poly_list->poly[i]->vlist->vtx, sizeof(ts_vertex *),1,fh);
	324	fwrite(vesicle->poly_list->poly[i]->blist->bond, sizeof(ts_bond *),1,fh);
	325	}
	326	*/
	327
	328
	329	if(retval);
	330	fclose(fh);
	331	return vesicle;
	332	}
	333
	334
d7639a	335
SP	336	ts_bool print_vertex_list(ts_vertex_list *vlist){
	337	ts_uint i;
	338	printf("Number of vertices: %u\n",vlist->n);
	339	for(i=0;i<vlist->n;i++){
a6b1b5	340	printf("%u: %f %f %f\n",
8f6a69	341	vlist->vtx[i]->idx,vlist->vtx[i]->x, vlist->vtx[i]->y, vlist->vtx[i]->z);
d7639a	342	}
SP	343	return TS_SUCCESS;
	344	}
	345
	346	ts_bool print_vertex_neighbours(ts_vertex_list *vlist){
	347	ts_uint i,j;
a6b1b5	348	ts_vertex **vtx=vlist->vtx;
d7639a	349	printf("Vertex id(neigh no): (neighvertex coord) (neighvertex coord) ...\n");
SP	350	for(i=0;i<vlist->n;i++){
b5ba7b	351	printf("%u(%u): ",vtx[i]->idx,vtx[i]->neigh->n);
SP	352	for(j=0;j<vtx[i]->neigh->n;j++){
	353	printf("(%f,%f,%f)",vtx[i]->neigh->vtx[j]->x,
	354	vtx[i]->neigh->vtx[j]->y,vtx[i]->neigh->vtx[j]->z);
d7639a	355	}
SP	356	printf("\n");
	357	}
	358
	359	return TS_SUCCESS;
	360	}
	361
	362	ts_bool write_vertex_fcompat_file(ts_vertex_list vlist,ts_char filename){
a6b1b5	363	ts_vertex **vtx=vlist->vtx;
d7639a	364	ts_uint i;
SP	365	FILE *fh;
	366
	367	fh=fopen(filename, "w");
	368	if(fh==NULL){
	369	err("Cannot open file %s for writing");
	370	return TS_FAIL;
	371	}
	372	for(i=0;i<vlist->n;i++)
8f6a69	373	fprintf(fh," %E\t%E\t%E\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a	374
SP	375	fclose(fh);
	376	return TS_SUCCESS;
	377	}
	378
	379
	380	ts_bool fprint_vertex_list(FILE fh,ts_vertex_list vlist){
	381	ts_uint i,j;
	382	for(i=0;i<vlist->n;i++){
8f6a69	383	fprintf(fh," %.17E\t%.17E\t%.17E\t%u\n",vlist->vtx[i]->x,
SP	384	vlist->vtx[i]->y, vlist->vtx[i]->z,
20ae83	385	vlist->vtx[i]->neigh->n);
SP	386	for(j=0;j<vlist->vtx[i]->neigh->n;j++){
	387	fprintf(fh,"\t%u",(ts_uint)(vlist->vtx[i]->neigh->vtx[j]->idx));
a6b1b5	388	//-vlist->vtx+1));
d7639a	389	}
SP	390	fprintf(fh,"\n");
	391	}
	392	return TS_SUCCESS;
	393	}
	394
	395	ts_bool fprint_tristar(FILE fh, ts_vesicle vesicle){
	396	ts_uint i,j;
a6b1b5	397	for(i=0;i<vesicle->vlist->n;i++){
8f6a69	398	fprintf(fh,"\t%u",vesicle->vlist->vtx[i]->tristar_no);
SP	399	for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
	400	fprintf(fh,"\t%u",(ts_uint)(vesicle->vlist->vtx[i]->tristar[j]->idx));//-vesicle->tlist->tria+1));
d7639a	401	}
SP	402	fprintf(fh,"\n");
	403	}
	404	return TS_SUCCESS;
	405	}
	406
	407	ts_bool fprint_triangle_list(FILE fh, ts_vesicle vesicle){
a6b1b5	408	ts_triangle_list *tlist=vesicle->tlist;
d7639a	409	ts_uint i,j;
SP	410	for(i=0;i<tlist->n;i++){
41a035	411	fprintf(fh,"\t%u",tlist->tria[i]->neigh_no);
SP	412	for(j=0;j<tlist->tria[i]->neigh_no;j++){
	413	fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->neigh[j]->idx));//-tlist->tria+1));
d7639a	414	}
SP	415	fprintf(fh,"\n");
	416	for(j=0;j<3;j++){
41a035	417	fprintf(fh,"\t%u",(ts_uint)(tlist->tria[i]->vertex[j]->idx));//-vesicle->vlist->vtx+1));
d7639a	418	}
SP	419	fprintf(fh,"\n");
41a035	420	fprintf(fh,"%.17E\t%.17E\t%.17E\n",tlist->tria[i]->xnorm,
SP	421	tlist->tria[i]->ynorm,tlist->tria[i]->znorm);
d7639a	422	fprintf(fh,"0.00000000000000000\n0.00000000000000000\n");
SP	423	}
	424	return TS_SUCCESS;
	425	}
	426
	427	ts_bool fprint_vertex_data(FILE fh,ts_vertex_list vlist){
	428	ts_uint i,j;
	429	for(i=0;i<vlist->n;i++){
	430	fprintf(fh," %.17E\t%.17E\t%.17E\t%.17E\t%.17E\t%u\n",
8f6a69	431	vlist->vtx[i]->xk,vlist->vtx[i]->c,vlist->vtx[i]->energy,
SP	432	vlist->vtx[i]->energy_h, vlist->vtx[i]->curvature, 0);
20ae83	433	for(j=0;j<vlist->vtx[i]->neigh->n;j++){
8f6a69	434	fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length_dual);
d7639a	435	}
SP	436	fprintf(fh,"\n");
20ae83	437	for(j=0;j<vlist->vtx[i]->neigh->n;j++){
8f6a69	438	fprintf(fh," %.17E", vlist->vtx[i]->bond[j]->bond_length);
d7639a	439	}
SP	440	fprintf(fh,"\n");
	441	}
	442	return TS_SUCCESS;
	443	}
	444
	445	ts_bool fprint_bonds(FILE fh,ts_vesicle vesicle){
	446	ts_uint i;
a6b1b5	447	for(i=0;i<vesicle->blist->n;i++){
e016c4	448	fprintf(fh,"\t%u\t%u\n",(ts_uint)(vesicle->blist->bond[i]->vtx1->idx),
a6b1b5	449	//-vesicle->vlist->vtx+1),
e016c4	450	(ts_uint)(vesicle->blist->bond[i]->vtx2->idx));
a6b1b5	451	//-vesicle->vlist.vtx+1));
d7639a	452	}
SP	453	return TS_SUCCESS;
	454	}
	455
	456
	457	ts_bool write_dout_fcompat_file(ts_vesicle vesicle, ts_char filename){
	458	FILE *fh;
	459	fh=fopen(filename, "w");
	460	if(fh==NULL){
	461	err("Cannot open file %s for writing");
	462	return TS_FAIL;
	463	}
	464	fprintf(fh,"%.17E\n%.17E\n",vesicle->stepsize,vesicle->dmax);
a6b1b5	465	fprint_vertex_list(fh,vesicle->vlist);
d7639a	466	fprint_tristar(fh,vesicle);
SP	467	fprint_triangle_list(fh,vesicle);
a6b1b5	468	fprint_vertex_data(fh,vesicle->vlist);
d7639a	469	fprint_bonds(fh,vesicle);
SP	470	fclose(fh);
	471	return TS_SUCCESS;
	472	}
	473
	474	ts_bool read_tape_fcompat_file(ts_vesicle vesicle, ts_char filename){
	475	FILE *fh;
	476	char line[255];
	477	fh=fopen(filename, "r");
	478	if(fh==NULL){
	479	err("Cannot open file for reading... Nonexistant file?");
	480	return TS_FAIL;
	481	}
a6b1b5	482	ts_uint retval=1;
d7639a	483	while(retval!=EOF){
SP	484	retval=fscanf(fh,"%s",line);
	485
	486	fprintf(stderr,"%s",line);
	487	}
	488	fclose(fh);
	489	return TS_SUCCESS;
	490	}
	491
	492	ts_bool write_master_xml_file(ts_char *filename){
	493	FILE *fh;
	494	ts_char i,j;
	495	ts_uint tstep;
	496	ts_char *number;
	497	fh=fopen(filename, "w");
	498	if(fh==NULL){
	499	err("Cannot open file %s for writing");
	500	return TS_FAIL;
	501	}
	502
	503	fprintf(fh,"<?xml version=\"1.0\"?>\n<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n<Collection>");
	504	DIR *dir = opendir(".");
	505	if(dir){
	506	struct dirent *ent;
	507	tstep=0;
	508	while((ent = readdir(dir)) != NULL)
	509	{
	510	i=rindex(ent->d_name,'.');
	511	if(i==NULL) continue;
	512	if(strcmp(i+1,"vtu")==0){
	513	j=rindex(ent->d_name,'_');
	514	if(j==NULL) continue;
	515	number=strndup(j+1,j-i);
a6b1b5	516	fprintf(fh,"<DataSet timestep=\"%u\" group=\"\" part=\"0\" file=\"%s\"/>\n",atoi(number),ent->d_name);
d7639a	517	tstep++;
SP	518	free(number);
	519	}
	520	}
	521	}
f74313	522	free(dir);
d7639a	523	fprintf(fh,"</Collection>\n</VTKFile>\n");
SP	524	fclose(fh);
	525	return TS_SUCCESS;
	526	}
	527
	528	ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno){
a6b1b5	529	ts_vertex_list *vlist=vesicle->vlist;
SP	530	ts_bond_list *blist=vesicle->blist;
	531	ts_vertex **vtx=vlist->vtx;
40aa5b	532	ts_uint i,j;
d7639a	533	char filename[255];
SP	534	FILE *fh;
	535
	536	sprintf(filename,"timestep_%.6u.vtu",timestepno);
	537	fh=fopen(filename, "w");
	538	if(fh==NULL){
	539	err("Cannot open file %s for writing");
	540	return TS_FAIL;
	541	}
	542	/* Here comes header of the file */
40aa5b	543
SP	544	//find number of extra vtxs and bonds of polymeres
	545	ts_uint monono=0, polyno=0;
	546	ts_bool poly=0;
	547	if(vesicle->poly_list!=NULL){
	548	if(vesicle->poly_list->poly[0]!=NULL){
	549	polyno=vesicle->poly_list->n;
	550	monono=vesicle->poly_list->poly[0]->vlist->n;
	551	poly=1;
	552	}
	553	}
d7639a	554	fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n <UnstructuredGrid>\n");
40aa5b	555	fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+mononopolyno, blist->n+mononopolyno);
d7639a	556	fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
SP	557	for(i=0;i<vlist->n;i++){
a6b1b5	558	fprintf(fh,"%u ",vtx[i]->idx);
d7639a	559	}
40aa5b	560	//polymeres
SP	561	if(poly){
	562	for(i=0;i<vesicle->poly_list->n;i++){
	563	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
	564	fprintf(fh,"%u ", vesicle->poly_list->poly[i]->vlist->vtx[j]->idx);
	565	}
	566	}
	567	}
d7639a	568
SP	569	fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
	570	for(i=0;i<vlist->n;i++){
8f6a69	571	fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
40aa5b	572	}
SP	573	//polymeres
	574	if(poly){
	575	for(i=0;i<vesicle->poly_list->n;i++){
	576	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
	577	fprintf(fh,"%e %e %e\n", vesicle->poly_list->poly[i]->vlist->vtx[j]->x,vesicle->poly_list->poly[i]->vlist->vtx[j]->y, vesicle->poly_list->poly[i]->vlist->vtx[j]->z );
	578	}
	579	}
d7639a	580	}
SP	581
	582	fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
	583	for(i=0;i<blist->n;i++){
e016c4	584	fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a	585	}
40aa5b	586	//polymeres
SP	587	if(poly){
	588	for(i=0;i<vesicle->poly_list->n;i++){
	589	for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
	590	fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx);
	591	}
	592	//grafted bonds
	593	fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->grafted_vtx->idx, vesicle->poly_list->poly[i]->vlist->vtx[0]->idx);
	594	}
	595
	596	}
	597
	598
d7639a	599	fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
40aa5b	600	for (i=2;i<(blist->n+mononopolyno)2+1;i+=2){
d7639a	601	fprintf(fh,"%u ",i);
SP	602	}
	603	fprintf(fh,"\n");
	604	fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
40aa5b	605	for (i=0;i<blist->n+monono*polyno;i++){
d7639a	606	fprintf(fh,"3 ");
SP	607	}
	608
	609	fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
	610	fclose(fh);
	611	return TS_SUCCESS;
	612
	613	}
	614
	615	ts_bool write_vertex_vtk_file(ts_vesicle vesicle,ts_char filename, ts_char *text){
a6b1b5	616	ts_vertex_list *vlist=vesicle->vlist;
SP	617	ts_bond_list *blist=vesicle->blist;
	618	ts_vertex **vtx=vlist->vtx;
	619	ts_uint i;
d7639a	620	FILE *fh;
SP	621
	622	fh=fopen(filename, "w");
	623	if(fh==NULL){
	624	err("Cannot open file %s for writing");
	625	return TS_FAIL;
	626	}
	627	/* Here comes header of the file */
	628	// fprintf(stderr,"NSHELL=%u\n",nshell);
	629	fprintf(fh, "# vtk DataFile Version 2.0\n");
	630	/* TODO: Do a sanity check on text. Max 255 char, must not me \n terminated */
	631	fprintf(fh, "%s\n", text);
	632	fprintf(fh,"ASCII\n");
	633	fprintf(fh,"DATASET UNSTRUCTURED_GRID\n");
	634	fprintf(fh,"POINTS %u double\n", vlist->n);
	635	for(i=0;i<vlist->n;i++){
8f6a69	636	fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
d7639a	637	}
SP	638
	639	fprintf(fh,"CELLS %u %u\n",blist->n,3*blist->n);
	640	for(i=0;i<blist->n;i++){
e016c4	641	fprintf(fh,"2 %u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
d7639a	642	}
SP	643	fprintf(fh,"CELL_TYPES %u\n",blist->n);
	644	for(i=0;i<blist->n;i++)
	645	fprintf(fh,"3\n");
	646
	647	fprintf(fh,"POINT_DATA %u\n", vlist->n);
	648	fprintf(fh,"SCALARS scalars long 1\n");
	649	fprintf(fh,"LOOKUP_TABLE default\n");
	650
	651	for(i=0;i<vlist->n;i++)
8f6a69	652	fprintf(fh,"%u\n",vtx[i]->idx);
d7639a	653
SP	654	fclose(fh);
	655	return TS_SUCCESS;
	656	}
	657
	658
	659
d7a113	660	ts_vesicle parsetape(ts_uint mcsweeps, ts_uint inititer, ts_uint iterations){
a2db52	661	long int nshell=17,ncxmax=60, ncymax=60, nczmax=60, npoly=10, nmono=20; // THIS IS DUE TO CONFUSE BUG!
b14a8d	662	char buf=malloc(255sizeof(char));
SP	663	long int brezveze0=1;
	664	long int brezveze1=1;
	665	long int brezveze2=1;
48bb92	666	ts_double xk0=25.0, dmax=1.67,stepsize=0.15,kspring=800.0;
d7a113	667	long int iter=1000, init=1000, mcsw=1000;
40aa5b	668
SP	669
d7639a	670	cfg_opt_t opts[] = {
SP	671	CFG_SIMPLE_INT("nshell", &nshell),
a2db52	672	CFG_SIMPLE_INT("npoly", &npoly),
SP	673	CFG_SIMPLE_INT("nmono", &nmono),
d7639a	674	CFG_SIMPLE_FLOAT("dmax", &dmax),
SP	675	CFG_SIMPLE_FLOAT("xk0",&xk0),
48bb92	676	CFG_SIMPLE_FLOAT("k_spring",&kspring),
d7639a	677	CFG_SIMPLE_FLOAT("stepsize",&stepsize),
SP	678	CFG_SIMPLE_INT("nxmax", &ncxmax),
	679	CFG_SIMPLE_INT("nymax", &ncymax),
	680	CFG_SIMPLE_INT("nzmax", &nczmax),
d7a113	681	CFG_SIMPLE_INT("iterations",&iter),
SP	682	CFG_SIMPLE_INT("mcsweeps",&mcsw),
	683	CFG_SIMPLE_INT("inititer", &init),
d7639a	684	CFG_SIMPLE_BOOL("quiet",&quiet),
314f2d	685	CFG_SIMPLE_STR("multiprocessing",buf),
b14a8d	686	CFG_SIMPLE_INT("smp_cores",&brezveze0),
SP	687	CFG_SIMPLE_INT("cluster_nodes",&brezveze1),
	688	CFG_SIMPLE_INT("distributed_processes",&brezveze2),
d7639a	689	CFG_END()
SP	690	};
	691	cfg_t *cfg;
	692	ts_uint retval;
	693	cfg = cfg_init(opts, 0);
a6b1b5	694	retval=cfg_parse(cfg, "tape");
d7639a	695	if(retval==CFG_FILE_ERROR){
a6b1b5	696	fatal("No tape file.",100);
d7639a	697	}
SP	698	else if(retval==CFG_PARSE_ERROR){
	699	fatal("Invalid tape!",100);
	700	}
b14a8d	701	ts_vesicle *vesicle;
d7a113	702	*iterations=iter;
SP	703	*inititer=init;
	704	*mcsweeps=mcsw;
b14a8d	705	vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize);
a2db52	706	vesicle->poly_list=init_poly_list(npoly,nmono, vesicle->vlist);
48bb92	707	vesicle->spring_constant=kspring;
M	708	poly_assign_spring_const(vesicle);
	709
a2db52	710
d7639a	711	vesicle->nshell=nshell;
SP	712	vesicle->dmax=dmax*dmax;
	713	vesicle->bending_rigidity=xk0;
	714	vesicle->stepsize=stepsize;
a6b1b5	715	vesicle->clist->ncmax[0]=ncxmax;
SP	716	vesicle->clist->ncmax[1]=ncymax;
	717	vesicle->clist->ncmax[2]=nczmax;
	718	vesicle->clist->max_occupancy=8;
b14a8d	719
d7639a	720	cfg_free(cfg);
b14a8d	721	free(buf);
SP	722	// fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
40aa5b	723
SP	724
	725
b14a8d	726	return vesicle;
d7639a	727
SP	728	}
f74313	729
SP	730
	731	ts_bool read_geometry_file(char fname, ts_vesicle vesicle){
	732	FILE *fh;
	733	ts_uint i, nvtx,nedges,ntria;
	734	ts_uint vtxi1,vtxi2;
	735	float x,y,z;
	736	ts_vertex_list *vlist;
	737	fh=fopen(fname, "r");
	738	if(fh==NULL){
	739	err("Cannot open file for reading... Nonexistant file?");
	740	return TS_FAIL;
	741	}
	742	ts_uint retval;
	743	retval=fscanf(fh,"%u %u %u",&nvtx, &nedges, &ntria);
	744	vesicle->vlist=init_vertex_list(nvtx);
	745	vlist=vesicle->vlist;
	746	for(i=0;i<nvtx;i++){
8f6a69	747	// fscanf(fh,"%F %F %F",&vlist->vtx[i]->x,&vlist->vtx[i]->y,&vlist->vtx[i]->z);
f74313	748	retval=fscanf(fh,"%F %F %F",&x,&y,&z);
8f6a69	749	vlist->vtx[i]->x=x;
SP	750	vlist->vtx[i]->y=y;
	751	vlist->vtx[i]->z=z;
f74313	752	}
SP	753	for(i=0;i<nedges;i++){
	754	retval=fscanf(fh,"%u %u",&vtxi1,&vtxi2);
	755	bond_add(vesicle->blist,vesicle->vlist->vtx[vtxi1-1],vesicle->vlist->vtx[vtxi2-1]);
	756	}
	757	//TODO: neighbours from bonds,
	758	//TODO: triangles from neigbours
	759
	760	// Don't need to read triangles. Already have enough data
	761	/*
	762	for(i=0;i<ntria;i++){
	763	retval=fscanf(fh,"%u %u %u", &bi1, &bi2, &bi3);
8f6a69	764	vtxi1=vesicle->blist->vertex1->idx;
SP	765	vtxi2=vesicle->blist->vertex1->idx;
f74313	766
SP	767	}
	768	*/
41a035	769	if(retval);
f74313	770	fclose(fh);
SP	771
	772
	773
	774	return TS_SUCCESS;
	775	}