trisurf-ng.git - Gitblit

Moved compression/decompression routines to separate file

Samo Penic

2019-10-19 dba66b4f8a25cb7eb18de5f7ffca97e179d7a5d9

commit \| author \| age
7f6076	1	/* vim: set ts=4 sts=4 sw=4 noet : */
f74313	2	#include "general.h"
d7639a	3	#include<stdio.h>
SP	4	#include "io.h"
f74313	5	#include "vertex.h"
SP	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"
055dd7	13	#include "cell.h"
083e03	14	#include <getopt.h>
SP	15	#include <sys/stat.h>
	16	#include <sys/types.h>
	17	#include <dirent.h>
	18	#include <errno.h>
854cb6	19	#include <snapshot.h>
dba66b	20	#include <b64zlib_compression.h>
e9eab4	21
083e03	22	ts_bool parse_args(int argc, char **argv){
3f5c83	23	int c, retval;
SP	24	struct stat sb;
	25	sprintf(command_line_args.path, "./"); //clear string;
7b0c07	26	sprintf(command_line_args.output_fullfilename,"output.pvd");
SP	27	sprintf(command_line_args.dump_fullfilename,"dump.bin");
	28	sprintf(command_line_args.tape_fullfilename,"tape");
1bf3c3	29	sprintf(command_line_args.tape_templatefull,"./tape");
3f5c83	30	FILE *file;
SP	31
083e03	32	while (1)
SP	33	{
	34	static struct option long_options[] =
	35	{
8a6614	36	{"force-from-tape", no_argument, &(command_line_args.force_from_tape), 1},
SP	37	{"reset-iteration-count", no_argument, &(command_line_args.reset_iteration_count), 1},
083e03	38	{"tape", no_argument, 0, 't'},
fd8126	39	{"version", no_argument, 0, 'v'},
083e03	40	{"output-file", required_argument, 0, 'o'},
SP	41	{"directory", required_argument, 0, 'd'},
3f5c83	42	{"dump-filename", required_argument,0, 'f'},
7b0c07	43	{"tape-options",required_argument,0,'c'},
1bf3c3	44	{"tape-template", required_argument,0,0},
ab798b	45	{"restore-from-vtk",required_argument,0,0},
083e03	46	{0, 0, 0, 0}
SP	47	};
	48	/* getopt_long stores the option index here. */
	49	int option_index = 0;
	50
fd8126	51	c = getopt_long (argc, argv, "d:f:o:t:c:v",
083e03	52	long_options, &option_index);
SP	53
	54	/* Detect the end of the options. */
	55	if (c == -1)
	56	break;
	57
	58	switch (c)
	59	{
	60	case 0:
	61	/* If this option set a flag, do nothing else now. */
	62	if (long_options[option_index].flag != 0)
	63	break;
1bf3c3	64	/* printf ("option %s", long_options[option_index].name);
083e03	65	if (optarg)
1bf3c3	66	printf (" with arg %s", optarg);
SP	67	printf ("\n"); */
	68	//TODO: find a better way.
	69	if(strcmp(long_options[option_index].name,"tape-template")==0){
	70	strcpy(command_line_args.tape_templatefull,optarg);
	71	}
ab798b	72	if(strcmp(long_options[option_index].name,"restore-from-vtk")==0){
ee84bd	73	strcpy(command_line_args.dump_from_vtk,optarg);
SP	74	}
083e03	75	break;
fd8126	76	case 'v':
SP	77	fprintf(stdout,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);
	78	fprintf(stdout,"Programming done by: Samo Penic and Miha Fosnaric\n");
	79	fprintf(stdout,"Released under terms of GPLv3\n");
	80	exit(0);
083e03	81
7b0c07	82	case 'c':
SP	83	strcpy(command_line_args.tape_opts,optarg);
	84	break;
	85	case 't': //tape
3f5c83	86	strcpy(command_line_args.tape_fullfilename,optarg);
083e03	87	break;
SP	88
7b0c07	89	case 'o': //set filename of master pvd output file
3f5c83	90	strcpy(command_line_args.output_fullfilename, optarg);
SP	91	break;
083e03	92
SP	93	case 'd':
	94	//check if directory exists. If not create one. If creation is
	95	//successful, set directory for output files.
	96	//printf ("option -d with value `%s'\n", optarg);
3f5c83	97	if (stat(optarg, &sb) == -1) {
SP	98	//directory does not exist
083e03	99	retval=mkdir(optarg, 0700);
SP	100	if(retval){
3f5c83	101	fatal("Could not create requested directory. Check if you have permissions",1);
083e03	102	}
SP	103	}
3f5c83	104	//check if is a proper directory
SP	105	else if((sb.st_mode & S_IFMT) != S_IFDIR) {
	106	//it is not a directory. fatal error.
	107	ts_fprintf(stderr,"%s is not a directory!\n",optarg);
	108	fatal("Cannot continue",1);
083e03	109	}
3f5c83	110	strcpy(command_line_args.path, optarg);
083e03	111	break;
SP	112
	113	case 'f':
3f5c83	114	strcpy(command_line_args.dump_fullfilename, optarg);
083e03	115	break;
SP	116
	117	case '?':
	118	/* getopt_long already printed an error message. */
ba73ab	119	print_help(stdout);
SP	120	//ts_fprintf(stderr,"\n\nhere comes the help.\n\n");
083e03	121	fatal("Ooops, read help first",1);
SP	122	break;
	123
	124	default:
	125	exit (1);
	126	}
	127	}
	128
7b0c07	129	//Here we set correct values for full filenames!
SP	130	char buffer=(char )malloc(10000*sizeof(char));
	131	//correct the path and add trailing /
	132	if(command_line_args.path[strlen(command_line_args.path)-1]!='/') strcat(command_line_args.path,"/");
	133
	134	/* master pvd output file */
	135	strcpy(buffer,command_line_args.path);
	136	strcat(buffer,command_line_args.output_fullfilename);
	137	if ((file = fopen(buffer, "w")) == NULL) {
	138	fprintf(stderr,"Could not create output file %s!\n", buffer);
	139	fatal("Please specify correct output file or check permissions of the file",1);
1bf3c3	140	//there is a tape template. make a copy into desired directory
7b0c07	141
SP	142	} else {
	143	fclose(file);
	144	strcpy(command_line_args.output_fullfilename,buffer);
	145	}
	146
	147	/* tape file */
	148	strcpy(buffer,command_line_args.path);
	149	strcat(buffer,command_line_args.tape_fullfilename);
	150	if (stat(buffer, &sb) == -1) {
1bf3c3	151
SP	152	//tape does not exist. does tape template exist?
	153	if(stat(command_line_args.tape_templatefull, &sb)==-1){
	154	ts_fprintf(stderr,"Tape '%s' does not exist and no tape template was specified (or does not exist)!\n",buffer);
	155	fatal("Please select correct tape or check permissions of the file",1);
	156	} else {
	157	//tape template found
	158	fatal("Samo did not program template copy yet",1);
	159	}
7b0c07	160	} else {
SP	161	strcpy(command_line_args.tape_fullfilename,buffer);
	162	}
	163
	164
	165	/* dump file */
	166	strcpy(buffer,command_line_args.path);
	167	strcat(buffer,command_line_args.dump_fullfilename);
	168	//check if dump file exist first.
	169	if (stat(buffer, &sb) == -1) {
	170	//no dump file. check if we can create one.
	171	if ((file = fopen(buffer, "w")) == NULL) {
	172	fprintf(stderr,"Could not create dump file '%s'!\n",buffer);
	173	fatal("Please specify correct dump file or check permissions of the file",1);
	174	} else {
	175	fclose(file);
	176	//good, file is writeable. delete it for now.
	177	remove(buffer);
	178	}
	179	}
	180	strcpy(command_line_args.dump_fullfilename, buffer);
	181
	182
	183	free(buffer);
083e03	184	return TS_SUCCESS;
SP	185
	186	}
	187
	188
ba73ab	189	ts_bool print_help(FILE *fd){
SP	190	fprintf(fd,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);
	191	fprintf(fd,"Programming done by: Samo Penic and Miha Fosnaric\n");
	192	fprintf(fd,"Released under terms of GPLv3\n\n");
	193
	194	fprintf(fd, "Invoking trisurf-ng without any flags results in default operation. Program reads 'tape' file and 'dump.bin' binary representation of the vesicle from disk and continues the simulation where it was aborted (as specified in 'dump.bin').\n\n");
	195	fprintf(fd,"If 'tape' has different values than binary dump, those are used (if possible -- some values in tape, such as nshell cannot be modified).\n\n\n");
	196	fprintf(fd,"However, if dump.bin is not present, user is notified to specify --force-from-tape flag. The vesicle will be created from specifications in tape.\n\n\n");
	197	fprintf(fd,"Flags:\n\n");
	198	fprintf(fd,"--force-from-tape\t\t makes initial shape of the vesicle from tape. Ignores already existing binary dump and possible simulation results.\n");
	199	fprintf(fd,"--restore-from-vtk\t\t VTK's file ending with '.vtu' are preferred way to make state snapshots for restoration. With this flag the restoration of the vesicle from vtk is possible. The simulation will continue if hidden '.status' file with last iteration done is present. Otherwise it will start simulation from timestep 0.\n");
c4178c	200	fprintf(fd,"--reset-iteration-count\t\t starts simulation from the beginning (using binary dump).\n");
ba73ab	201	fprintf(fd,"--tape (or -t)\t\t specifies tape filename. For --force-from-tape and restoring from binary dump. Defaults to 'tape'.\n");
SP	202	fprintf(fd,"--version (or -v)\t\t Prints version information.\n");
	203	fprintf(fd,"--output-file (or -o)\t\t Specifies filename of .PVD file. Defaults to 'output.pvd'\n");
c4178c	204	fprintf(fd,"--dump-filename (or -f)\t\t specifies filename for binary dump&restore. Defaults to 'dump.bin'\n\n\n");
SP	205	fprintf(fd,"Examples:\n\n");
	206	fprintf(fd,"trisurf --force-from-tape\n");
	207	fprintf(fd,"trisurf --reset-iteration-count\n");
	208	fprintf(fd,"trisurf --restore-from-vtk filename.vtu\n");
	209	fprintf(fd,"\n\n");
	210
d7639a	211	return TS_SUCCESS;
SP	212	}
	213
0a2c81	214	ts_bool write_vertex_xml_file(ts_vesicle vesicle, ts_uint timestepno, ts_cluster_list cstlist){
a6b1b5	215	ts_vertex_list *vlist=vesicle->vlist;
SP	216	ts_bond_list *blist=vesicle->blist;
	217	ts_vertex **vtx=vlist->vtx;
40aa5b	218	ts_uint i,j;
13d445	219	//ts_double senergy=0.0;
267db5	220	char filename[10000];
SP	221	char just_name[255];
d7639a	222	FILE *fh;
267db5	223	strcpy(filename,command_line_args.path);
SP	224	sprintf(just_name,"timestep_%.6u.vtu",timestepno);
	225	strcat(filename,just_name);
d7639a	226
SP	227	fh=fopen(filename, "w");
	228	if(fh==NULL){
	229	err("Cannot open file %s for writing");
	230	return TS_FAIL;
	231	}
	232	/* Here comes header of the file */
40aa5b	233
SP	234	//find number of extra vtxs and bonds of polymeres
58230a	235	ts_uint monono=0, polyno=0, poly_idx=0, filno=0, fonono=0;
M	236	ts_bool poly=0, fil=0;
40aa5b	237	if(vesicle->poly_list!=NULL){
SP	238	if(vesicle->poly_list->poly[0]!=NULL){
	239	polyno=vesicle->poly_list->n;
	240	monono=vesicle->poly_list->poly[0]->vlist->n;
	241	poly=1;
	242	}
	243	}
58230a	244
M	245	if(vesicle->filament_list!=NULL){
	246	if(vesicle->filament_list->poly[0]!=NULL){
	247	filno=vesicle->filament_list->n;
	248	fonono=vesicle->filament_list->poly[0]->vlist->n;
	249	fil=1;
	250	}
	251	}
	252
854cb6	253	fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n");
SP	254	xml_trisurf_data(fh,vesicle);
	255	fprintf(fh, " <UnstructuredGrid>\n");
8db203	256	fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+mononopolyno+fononofilno, blist->n+mononopolyno+filno(fonono-1)+vesicle->tlist->n);
80ebbe	257	fprintf(fh,"<PointData Scalars=\"vertices_idx\">\n<DataArray type=\"Int64\" Name=\"vertices_idx\" format=\"binary\">");
a69203	258	long int_vector=(long )malloc((vlist->n+mononopolyno+fononofilno)*sizeof(long));
842378	259	ts_double float_vector=(double )malloc((vlist->n+mononopolyno+fononofilno)*sizeof(double));
SP	260	ts_double float_cell_vector=(double )malloc((blist->n+mononopolyno+filno(fonono-1)+vesicle->tlist->n)*sizeof(double));
	261	ts_double coordinate_vector=(double )malloc(3(vlist->n+mononopolyno+fononofilno)sizeof(double));
a69203	262	// sizes are as expected 4 bytes for int and 1 byte for char
SP	263	// printf("size of ts_uint %ld, of int %ld, of char %ld",sizeof(ts_uint), sizeof(int), sizeof(char));
80ebbe	264	int offset=0;
d7639a	265	for(i=0;i<vlist->n;i++){
80ebbe	266	// fprintf(fh,"%u ",vtx[i]->idx);
SP	267	int_vector[i+offset]=vtx[i]->idx;
d7639a	268	}
a69203	269	offset=offset+i;
40aa5b	270	//polymeres
SP	271	if(poly){
3c1ac1	272	poly_idx=vlist->n;
40aa5b	273	for(i=0;i<vesicle->poly_list->n;i++){
3c1ac1	274	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
80ebbe	275	//fprintf(fh,"%u ", poly_idx);
SP	276	int_vector[j+offset]=poly_idx;
58230a	277	}
80ebbe	278	offset=offset+j;
58230a	279	}
M	280	}
	281	//filaments
	282	if(fil){
	283	poly_idx=vlist->n+monono*polyno;
	284	for(i=0;i<vesicle->filament_list->n;i++){
	285	for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
	286	// fprintf(stderr,"was here\n");
80ebbe	287	//fprintf(fh,"%u ", poly_idx);
SP	288	int_vector[j+offset]=poly_idx;
40aa5b	289	}
80ebbe	290	offset=offset+j;
40aa5b	291	}
SP	292	}
a69203	293	printf("Offset in the end is %d, should be %d",offset,(vlist->n+mononopolyno+fononofilno) );
842378	294	char printout=ts_compress((char )int_vector,offset*sizeof(long));
80ebbe	295	fprintf(fh,"%s",printout);
SP	296	free(printout);
d7639a	297
45c708	298	fprintf(fh,"</DataArray>\n");
0a2c81	299	if(cstlist!=NULL){
842378	300	fprintf(fh,"<DataArray type=\"Int64\" Name=\"vertices_in_cluster\" format=\"binary\">");
SP	301	offset=0;
	302	for(i=0;i<vlist->n;i++,offset++){
0a2c81	303	if(vtx[i]->cluster!=NULL){
842378	304	// fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
SP	305	int_vector[offset]=vtx[i]->cluster->nvtx;
0a2c81	306	} else {
842378	307	// fprintf(fh,"-1 ");
SP	308	int_vector[offset]=-1;
0a2c81	309	}
SP	310	}
	311	//polymeres
	312	if(poly){
	313	poly_idx=vlist->n;
	314	for(i=0;i<vesicle->poly_list->n;i++){
842378	315	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
SP	316	//fprintf(fh,"-1 ");
	317	int_vector[offset]=-1;
0a2c81	318	}
SP	319	}
	320	}
	321	//filaments
	322	if(fil){
	323	poly_idx=vlist->n+monono*polyno;
	324	for(i=0;i<vesicle->filament_list->n;i++){
842378	325	for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
SP	326	// fprintf(fh,"-1 ");
	327	int_vector[offset]=-1;
0a2c81	328	}
SP	329	}
	330	}
842378	331	char printout=ts_compress((char )int_vector,offset*sizeof(long));
SP	332	fprintf(fh,"%s",printout);
	333	free(printout);
0a2c81	334
SP	335	fprintf(fh,"</DataArray>\n");
	336
	337
	338	}
	339
45c708	340	//here comes additional data as needed. Currently only spontaneous curvature
842378	341	fprintf(fh,"<DataArray type=\"Float64\" Name=\"spontaneous_curvature\" format=\"binary\">");
SP	342	offset=0;
	343	for(i=0;i<vlist->n;i++, offset++){
	344	// fprintf(fh,"%.17e ",vtx[i]->c);
	345	float_vector[offset]=vtx[i]->c;
45c708	346	}
SP	347	//polymeres
	348	if(poly){
	349	poly_idx=vlist->n;
	350	for(i=0;i<vesicle->poly_list->n;i++){
842378	351	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
SP	352	//fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->c);
	353	float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->c;
45c708	354	}
SP	355	}
	356	}
	357	//filaments
	358	if(fil){
	359	poly_idx=vlist->n+monono*polyno;
	360	for(i=0;i<vesicle->filament_list->n;i++){
842378	361	for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
45c708	362	// fprintf(stderr,"was here\n");
842378	363	//fprintf(fh,"%.17e ", vesicle->filament_list->poly[i]->vlist->vtx[j]->c);
SP	364	float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->c;
45c708	365	}
SP	366	}
	367	}
842378	368	printout=ts_compress((char )float_vector,offsetsizeof(double));
SP	369	fprintf(fh,"%s",printout);
	370	free(printout);
	371
	372
45c708	373	fprintf(fh,"</DataArray>\n");
SP	374
5c64e2	375	//here comes additional data. Energy!
842378	376	fprintf(fh,"<DataArray type=\"Float64\" Name=\"bending_energy\" format=\"binary\">");
SP	377	offset=0;
	378	for(i=0;i<vlist->n;i++,offset++){
	379	// fprintf(fh,"%.17e ",vtx[i]->energy*vtx[i]->xk);
	380	float_vector[offset]=vtx[i]->energy*vtx[i]->xk;
5c64e2	381	}
SP	382	//polymeres
	383	if(poly){
	384	poly_idx=vlist->n;
	385	for(i=0;i<vesicle->poly_list->n;i++){
842378	386	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
SP	387	//fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k);
	388	float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k;
5c64e2	389	}
SP	390	}
	391	}
	392	//filaments
	393	if(fil){
	394	poly_idx=vlist->n+monono*polyno;
	395	for(i=0;i<vesicle->filament_list->n;i++){
842378	396	for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
5c64e2	397	// fprintf(stderr,"was here\n");
842378	398	// fprintf(fh,"%.17e ", vesicle->filament_list->poly[i]->vlist->vtx[j]->energy* vesicle->filament_list->poly[i]->k);
SP	399	float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->energy* vesicle->filament_list->poly[i]->k;
5c64e2	400	}
SP	401	}
	402	}
842378	403
SP	404	printout=ts_compress((char )float_vector,offsetsizeof(double));
	405	fprintf(fh,"%s",printout);
	406	free(printout);
	407
5c64e2	408	fprintf(fh,"</DataArray>\n");
SP	409
45c708	410
13d445	411	fprintf(fh,"</PointData>\n<CellData>\n");
SP	412
7ec6fb	413	if(vesicle->tape->stretchswitch==1){
842378	414	fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"binary\">");
SP	415	offset=0;
	416	for(i=0;i<blist->n;i++,offset++){
	417	//fprintf(fh, "0.0 ");
	418	float_cell_vector[offset]=0.0;
cf7aba	419	}
842378	420	for(i=0;i<mononopolyno+filno(fonono-1);i++,offset++){
SP	421	float_cell_vector[offset]=0.0;
	422	//fprintf(fh,"0.0 ");
7ec6fb	423	}
842378	424	for(i=0;i<vesicle->tlist->n;i++,offset++){
SP	425	float_cell_vector[offset]=vesicle->tlist->tria[i]->energy;
	426	//fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
13d445	427	}
842378	428
SP	429	printout=ts_compress((char )float_cell_vector,offsetsizeof(double));
	430	fprintf(fh,"%s",printout);
	431	free(printout);
	432
7ec6fb	433	fprintf(fh,"</DataArray>\n");
SP	434	}
45c708	435
13d445	436
SP	437
842378	438	fprintf(fh,"</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"binary\">\n");
SP	439	offset=0;
	440	for(i=0;i<vlist->n;i++,offset+=3){
	441	//fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
	442	coordinate_vector[offset]=vtx[i]->x;
	443	coordinate_vector[offset+1]=vtx[i]->y;
	444	coordinate_vector[offset+2]=vtx[i]->z;
40aa5b	445	}
SP	446	//polymeres
	447	if(poly){
	448	for(i=0;i<vesicle->poly_list->n;i++){
842378	449	for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,offset+=3){
SP	450	//fprintf(fh,"%.17e %.17e %.17e\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 );
	451	coordinate_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->x;
	452	coordinate_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[j]->y;
	453	coordinate_vector[offset+2]=vesicle->poly_list->poly[i]->vlist->vtx[j]->z;
58230a	454	}
M	455	}
	456	}
	457	//filaments
	458	if(fil){
842378	459	for(i=0;i<vesicle->filament_list->n;i++,offset+=3){
58230a	460	for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
842378	461	// fprintf(fh,"%.17e %.17e %.17e\n", vesicle->filament_list->poly[i]->vlist->vtx[j]->x,vesicle->filament_list->poly[i]->vlist->vtx[j]->y, vesicle->filament_list->poly[i]->vlist->vtx[j]->z );
SP	462	coordinate_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->x;
	463	coordinate_vector[offset+1]=vesicle->filament_list->poly[i]->vlist->vtx[j]->y;
	464	coordinate_vector[offset+2]=vesicle->filament_list->poly[i]->vlist->vtx[j]->z;
40aa5b	465	}
SP	466	}
d7639a	467	}
842378	468	printout=ts_compress((char )coordinate_vector,offsetsizeof(double));
SP	469	fprintf(fh,"%s",printout);
	470	free(printout);
d7639a	471
842378	472
SP	473
	474	fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"binary\">");
	475	free(int_vector);
	476	int_vector=(long )malloc(((blist->n+mononopolyno+filno(fonono-1))2+vesicle->tlist->n3)sizeof(long));
	477	offset=0;
	478	for(i=0;i<blist->n;i++,offset+=2){
	479	// fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
	480	int_vector[offset]=blist->bond[i]->vtx1->idx;
	481	int_vector[offset+1]=blist->bond[i]->vtx2->idx;
d7639a	482	}
40aa5b	483	//polymeres
SP	484	if(poly){
3c1ac1	485	poly_idx=vlist->n;
40aa5b	486	for(i=0;i<vesicle->poly_list->n;i++){
842378	487	for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++,offset+=2){
SP	488	//fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+imonono,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+imonono);
	489	int_vector[offset]= vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono;
	490	int_vector[offset+1]=vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono;
	491
40aa5b	492	}
SP	493	//grafted bonds
842378	494	//fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->grafted_vtx->idx, vesicle->poly_list->poly[i]->vlist->vtx[0]->idx+vlist->n+i*monono);
SP	495	int_vector[offset]=vesicle->poly_list->poly[i]->grafted_vtx->idx;
	496	int_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[0]->idx+vlist->n+i*monono;
	497	offset+=2;
40aa5b	498	}
SP	499
	500	}
	501
58230a	502	//filaments
M	503	if(fil){
	504	poly_idx=vlist->n+monono*polyno;
	505	for(i=0;i<vesicle->filament_list->n;i++){
842378	506	for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++,offset+=2){
SP	507	// fprintf(fh,"%u %u\n", vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+mononopolyno+ifonono,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+mononopolyno+ifonono);
	508	int_vector[offset]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+mononopolyno+ifonono;
	509	int_vector[offset+1]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+mononopolyno+ifonono;
58230a	510	// fprintf(stderr,"was here\n");
M	511
	512	}
	513	}
	514
	515	}
842378	516	for(i=0;i<vesicle->tlist->n;i++,offset+=3){
SP	517	// fprintf(fh,"%u %u %u\n", vesicle->tlist->tria[i]->vertex[0]->idx, vesicle->tlist->tria[i]->vertex[1]->idx, vesicle->tlist->tria[i]->vertex[2]->idx);
	518	int_vector[offset]=vesicle->tlist->tria[i]->vertex[0]->idx;
	519	int_vector[offset+1]=vesicle->tlist->tria[i]->vertex[1]->idx;
	520	int_vector[offset+2]=vesicle->tlist->tria[i]->vertex[2]->idx;
	521
8db203	522	}
842378	523	printout=ts_compress((char )int_vector,offsetsizeof(long));
SP	524	fprintf(fh,"%s",printout);
	525	free(printout);
	526
	527
	528	fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"binary\">");
	529	offset=0;
	530	for (i=2;i<(blist->n+mononopolyno+(fonono-1)filno)*2+1;i+=2,offset++){
	531	// fprintf(fh,"%u ",i);
	532	int_vector[offset]=i;
d7639a	533	}
842378	534	for(j=i+1;j<i+3*(vesicle->tlist->n);j+=3,offset++){ //let's continue counting from where we left of
SP	535	// fprintf(fh,"%u ", j);
	536	int_vector[offset]=j;
8db203	537	}
842378	538	// fprintf(fh,"\n");
SP	539	printout=ts_compress((char )int_vector,offsetsizeof(long));
	540	fprintf(fh,"%s",printout);
	541	free(printout);
	542	fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"binary\">\n");
	543	free(int_vector);
	544	free(float_vector);
	545	free(float_cell_vector);
	546	free(coordinate_vector);
	547	unsigned char char_vector=(unsigned char )malloc((blist->n+mononopolyno+filno(fonono-1)+vesicle->tlist->n)*sizeof(char));
	548	offset=0;
	549	for (i=0;i<blist->n+mononopolyno+(fonono-1)filno;i++,offset++){
	550	//fprintf(fh,"3 ");
	551	char_vector[offset]=3;
d7639a	552	}
842378	553	for(i=0;i<vesicle->tlist->n;i++,offset++){
SP	554	// fprintf(fh,"5 ");
	555	char_vector[offset]=5;
8db203	556	}
842378	557	printout=ts_compress((char )char_vector,offsetsizeof(char));
SP	558	fprintf(fh,"%s",printout);
	559	free(printout);
	560	free(char_vector);
	561
d7639a	562	fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
SP	563	fclose(fh);
	564	return TS_SUCCESS;
	565
	566	}
	567
144784	568	ts_bool write_pov_file(ts_vesicle vesicle, char filename){
SP	569	FILE *fh;
	570	ts_uint i;
	571
	572	fh=fopen(filename, "w");
	573	if(fh==NULL){
	574	err("Cannot open file %s for writing");
	575	return TS_FAIL;
	576	}
	577
	578	for(i=0;i<vesicle->tlist->n;i++){
	579
	580	fprintf(fh,"\ttriangle {");
	581	fprintf(fh,"\t<%e,%e,%e> <%e,%e,%e> <%e,%e,%e> }\n",
	582	vesicle->tlist->tria[i]->vertex[0]->x,
	583	vesicle->tlist->tria[i]->vertex[0]->y,
	584	vesicle->tlist->tria[i]->vertex[0]->z,
	585
	586	vesicle->tlist->tria[i]->vertex[1]->x,
	587	vesicle->tlist->tria[i]->vertex[1]->y,
	588	vesicle->tlist->tria[i]->vertex[1]->z,
	589
	590	vesicle->tlist->tria[i]->vertex[2]->x,
	591	vesicle->tlist->tria[i]->vertex[2]->y,
	592	vesicle->tlist->tria[i]->vertex[2]->z
	593	);
	594	}
	595
	596	fclose(fh);
	597	return TS_SUCCESS;
	598	}
	599
	600
1ab449	601	ts_tape parsetape(char filename){
698ae1	602	FILE *fd = fopen (filename, "r");
SP	603	long length;
	604	size_t size;
	605	fseek (fd, 0, SEEK_END);
	606	length = ftell (fd);
	607	fseek (fd, 0, SEEK_SET);
	608	size=fread (tapetxt, 1, length, fd);
	609	fclose(fd);
	610	if(size);
	611	ts_tape *tape=parsetapebuffer(tapetxt);
	612	return tape;
	613	}
	614
	615	ts_tape parsetapebuffer(char buffer){
1ab449	616	ts_tape tape=(ts_tape )calloc(1,sizeof(ts_tape));
SP	617	tape->multiprocessing=calloc(255,sizeof(char));
40aa5b	618
d7639a	619	cfg_opt_t opts[] = {
1ab449	620	CFG_SIMPLE_INT("nshell", &tape->nshell),
SP	621	CFG_SIMPLE_INT("npoly", &tape->npoly),
	622	CFG_SIMPLE_INT("nmono", &tape->nmono),
58230a	623	CFG_SIMPLE_INT("nfil",&tape->nfil),
M	624	CFG_SIMPLE_INT("nfono",&tape->nfono),
e98482	625	CFG_SIMPLE_INT("internal_poly",&tape->internal_poly),
58230a	626	CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
37791b	627	CFG_SIMPLE_FLOAT("R_nucleusX",&tape->R_nucleusX),
SP	628	CFG_SIMPLE_FLOAT("R_nucleusY",&tape->R_nucleusY),
	629	CFG_SIMPLE_FLOAT("R_nucleusZ",&tape->R_nucleusZ),
58230a	630	CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
ea1cce	631	CFG_SIMPLE_FLOAT("dmin_interspecies", &tape->dmin_interspecies),
1ab449	632	CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
SP	633	CFG_SIMPLE_INT("pswitch",&tape->pswitch),
9166cb	634	CFG_SIMPLE_INT("constvolswitch",&tape->constvolswitch),
c0ae90	635	CFG_SIMPLE_INT("constareaswitch",&tape->constareaswitch),
43c042	636	CFG_SIMPLE_FLOAT("constvolprecision",&tape->constvolprecision),
2ae815	637	CFG_SIMPLE_INT("stretchswitch",&tape->stretchswitch),
SP	638	CFG_SIMPLE_FLOAT("xkA0",&tape->xkA0),
1ab449	639	CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
SP	640	CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
b30f45	641	CFG_SIMPLE_FLOAT("xi",&tape->xi),
1ab449	642	CFG_SIMPLE_FLOAT("stepsize",&tape->stepsize),
SP	643	CFG_SIMPLE_INT("nxmax", &tape->ncxmax),
	644	CFG_SIMPLE_INT("nymax", &tape->ncymax),
	645	CFG_SIMPLE_INT("nzmax", &tape->nczmax),
	646	CFG_SIMPLE_INT("iterations",&tape->iterations),
	647	CFG_SIMPLE_INT("mcsweeps",&tape->mcsweeps),
	648	CFG_SIMPLE_INT("inititer", &tape->inititer),
819a09	649	CFG_SIMPLE_BOOL("quiet",(cfg_bool_t *)&tape->quiet),
S	650	CFG_SIMPLE_STR("multiprocessing",&tape->multiprocessing),
1ab449	651	CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
SP	652	CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
	653	CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
dc77e8	654	CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
e5858f	655	CFG_SIMPLE_INT("number_of_vertices_with_c0", &tape->number_of_vertices_with_c0),
SP	656	CFG_SIMPLE_FLOAT("c0",&tape->c0),
	657	CFG_SIMPLE_FLOAT("w",&tape->w),
250de4	658	CFG_SIMPLE_FLOAT("F",&tape->F),
6c274d	659	/* Variables related to plane confinement */
514ebc	660	CFG_INT("plane_confinement_switch", 0, CFGF_NONE),
SP	661	CFG_FLOAT("plane_d", 15, CFGF_NONE),
	662	CFG_FLOAT("plane_F", 1000, CFGF_NONE),
6c274d	663	/* Variables related to stretching */
36bc6d	664	// CFG_FLOAT("stretchswitch", 0, CFGF_NONE),
SP	665	// CFG_FLOAT("xkA0",0,CFGF_NONE),
d7639a	666	CFG_END()
SP	667	};
	668	cfg_t *cfg;
	669	ts_uint retval;
	670	cfg = cfg_init(opts, 0);
698ae1	671	retval=cfg_parse_buf(cfg, buffer);
514ebc	672	tape->plane_confinement_switch=cfg_getint(cfg,"plane_confinement_switch");
88bdd7	673	tape->plane_d=cfg_getfloat(cfg,"plane_d");
SP	674	tape->plane_F=cfg_getfloat(cfg,"plane_F");
514ebc	675
d7639a	676	if(retval==CFG_FILE_ERROR){
a6b1b5	677	fatal("No tape file.",100);
d7639a	678	}
SP	679	else if(retval==CFG_PARSE_ERROR){
	680	fatal("Invalid tape!",100);
	681	}
a2db52	682
7b0c07	683	/* here we override all values read from tape with values from commandline*/
SP	684	getcmdline_tape(cfg,command_line_args.tape_opts);
d7639a	685	cfg_free(cfg);
40aa5b	686
SP	687
1ab449	688	/* global variables are set automatically */
SP	689	quiet=tape->quiet;
	690	return tape;
	691	}
d7639a	692
1ab449	693	ts_bool tape_free(ts_tape *tape){
SP	694	free(tape->multiprocessing);
	695	free(tape);
	696	return TS_SUCCESS;
d7639a	697	}
f74313	698
SP	699
7b0c07	700
SP	701	ts_bool getcmdline_tape(cfg_t cfg, char opts){
	702
	703	char commands, backup, saveptr, saveopptr, command, operator[2];
07e3de	704	operator[0]=0;
SP	705	operator[1]=0;
7b0c07	706	ts_uint i,j;
SP	707	commands=(char )malloc(10000sizeof(char));
	708	backup=commands; //since the pointer to commands will be lost, we acquire a pointer that will serve as backup.
	709	strcpy(commands,opts);
	710	for(i=0; ;i++, commands=NULL){
	711	//breaks comma separated list of commands into specific commands.
	712	command=strtok_r(commands,",",&saveptr);
	713	if(command==NULL) break;
	714	// fprintf(stdout,"Command %d: %s\n",i,command);
	715	//extracts name of command and value of command into operator[2] array.
	716	for(j=0; j<2;j++,command=NULL){
	717	operator[j]=strtok_r(command,"=",&saveopptr);
	718	if(operator[j]==NULL) break;
	719	// fprintf(stdout," ---> Operator %d: %s\n",j,operator[j]);
	720	}
	721	//1. check: must have 2 operators.
	722	if(j!=2) fatal("Error. Command line tape options are not formatted properly",1);
	723
	724	// cfg_setstr(cfg,operator[0],operator[1]);
	725	cmdline_to_tape(cfg,operator[0],operator[1]);
	726	//2. check: must be named properly.
	727	//3. check: must be of right format (integer, double, string, ...)
	728
	729	}
	730	free(backup);
	731	return TS_SUCCESS;
	732	}
	733
	734
	735	ts_bool cmdline_to_tape(cfg_t cfg, char key, char *val){
	736
	737	cfg_opt_t *cfg_opt=cfg_getopt(cfg,key);
	738	if(cfg_opt==NULL) fatal("Commandline tape option not recognised",1); //return TS_FAIL;
	739	switch (cfg_opt->type){
	740	case CFGT_INT:
	741	cfg_setint(cfg,key,atol(val));
	742	break;
	743	case CFGT_FLOAT:
	744	cfg_setfloat(cfg,key,atof(val));
	745	break;
	746	/* case CFGT_BOOL:
	747	cfg_setbool(cfg,operator[0],operator[1]);
	748	break; */
	749	case CFGT_STR:
	750	cfg_setstr(cfg,key,val);
	751	break;
	752	default:
	753	break;
	754
	755	}
f74313	756	return TS_SUCCESS;
SP	757	}