From 2afc2f4f1dd89518995f1b5a539aea932aecab65 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Fri, 08 Mar 2019 19:50:18 +0000
Subject: [PATCH] Full write of constant volume plugin

---
 src/io.c |  127 ++++++++++++++++++++++++++++++++++++++---
 1 files changed, 116 insertions(+), 11 deletions(-)

diff --git a/src/io.c b/src/io.c
index 39dc1e0..30bd4d6 100644
--- a/src/io.c
+++ b/src/io.c
@@ -1,3 +1,4 @@
+
 /* vim: set ts=4 sts=4 sw=4 noet : */
 #include "general.h"
 #include<stdio.h>
@@ -423,7 +424,7 @@
     vesicle->tape=parsetape(command_line_args.tape_fullfilename);
 // recreating space for cells // 
     vesicle->clist=init_cell_list(vesicle->tape->ncxmax, vesicle->tape->ncymax, vesicle->tape->nczmax, vesicle->tape->stepsize);
-	vesicle->clist->max_occupancy=8;
+	vesicle->clist->max_occupancy=16;
 //    vesicle->tape=(ts_tape *)malloc(sizeof(ts_tape));
 //    retval=fread(vesicle->tape, sizeof(ts_tape),1,fh);
 	retval=fread(iteration,sizeof(ts_uint),1,fh);
@@ -820,11 +821,12 @@
 	return TS_SUCCESS;
 }
 
-ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno){
+ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno, ts_cluster_list *cstlist){
 	ts_vertex_list *vlist=vesicle->vlist;
 	ts_bond_list *blist=vesicle->blist;
 	ts_vertex **vtx=vlist->vtx;
     ts_uint i,j;
+	//ts_double senergy=0.0;
     	char filename[10000];
         char just_name[255];
 	FILE *fh;
@@ -861,8 +863,8 @@
 	fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n");
 	xml_trisurf_data(fh,vesicle);
 	fprintf(fh, " <UnstructuredGrid>\n");
-    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1));
-    fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
+    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n);
+    fprintf(fh,"<PointData Scalars=\"vertices_idx\">\n<DataArray type=\"Int64\" Name=\"vertices_idx\" format=\"ascii\">");
    	for(i=0;i<vlist->n;i++){
 		fprintf(fh,"%u ",vtx[i]->idx);
     }
@@ -887,7 +889,40 @@
 	}
 
     	fprintf(fh,"</DataArray>\n");
-	
+	if(cstlist!=NULL){
+		fprintf(fh,"<DataArray type=\"Int64\" Name=\"vertices_in_cluster\" format=\"ascii\">");
+		for(i=0;i<vlist->n;i++){
+			if(vtx[i]->cluster!=NULL){
+				fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
+			} else {
+				fprintf(fh,"-1 ");
+			}
+	    	}
+		//polymeres
+		if(poly){
+			poly_idx=vlist->n;
+			for(i=0;i<vesicle->poly_list->n;i++){
+				for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
+					fprintf(fh,"-1 ");
+				}
+			}
+		}
+		//filaments
+		if(fil){
+			poly_idx=vlist->n+monono*polyno;
+			for(i=0;i<vesicle->filament_list->n;i++){
+				for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
+		//	fprintf(stderr,"was here\n");
+					fprintf(fh,"-1 ");
+				}
+			}
+		}
+
+		fprintf(fh,"</DataArray>\n");
+
+
+	}
+
 	//here comes additional data as needed. Currently only spontaneous curvature
 	fprintf(fh,"<DataArray type=\"Float64\" Name=\"spontaneous_curvature\" format=\"ascii\">");
 	for(i=0;i<vlist->n;i++){
@@ -914,9 +949,52 @@
 		}
     fprintf(fh,"</DataArray>\n");
 
+	//here comes additional data. Energy!
+	fprintf(fh,"<DataArray type=\"Float64\" Name=\"bending_energy\" format=\"ascii\">");
+	for(i=0;i<vlist->n;i++){
+		fprintf(fh,"%.17e ",vtx[i]->energy*vtx[i]->xk);
+	}
+		//polymeres
+		if(poly){
+			poly_idx=vlist->n;
+			for(i=0;i<vesicle->poly_list->n;i++){
+				for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
+					fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k);
+				}
+			}
+		}
+		//filaments
+		if(fil){
+			poly_idx=vlist->n+monono*polyno;
+			for(i=0;i<vesicle->filament_list->n;i++){
+				for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
+		//	fprintf(stderr,"was here\n");
+					fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k);
+				}
+			}
+		}
+    fprintf(fh,"</DataArray>\n");
 
 	
-	fprintf(fh,"</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
+	fprintf(fh,"</PointData>\n<CellData>\n");
+
+	if(vesicle->tape->stretchswitch==1){
+		fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"ascii\">");
+		for(i=0;i<blist->n;i++){
+			fprintf(fh, "0.0 ");
+		}
+		for(i=0;i<monono*polyno+filno*(fonono-1);i++){
+			fprintf(fh,"0.0 ");
+		}
+		for(i=0;i<vesicle->tlist->n;i++){
+			fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
+		}
+	    fprintf(fh,"</DataArray>\n");
+	}
+
+
+
+	fprintf(fh,"</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
 	for(i=0;i<vlist->n;i++){
 		fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
 	}
@@ -967,17 +1045,24 @@
 		}
 
 	}
-
+	for(i=0;i<vesicle->tlist->n;i++){
+		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);
+	}
     fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
     for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2){
     fprintf(fh,"%u ",i);
     }
+	for(j=i+1;j<i+3*(vesicle->tlist->n);j+=3){ //let's continue counting from where we left of
+		fprintf(fh,"%u ", j);
+	}
     fprintf(fh,"\n");
     fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
-     for (i=0;i<blist->n+monono*polyno+fonono*filno;i++){
+     for (i=0;i<blist->n+monono*polyno+(fonono-1)*filno;i++){
         fprintf(fh,"3 ");
     }
-
+	for(i=0;i<vesicle->tlist->n;i++){
+		fprintf(fh,"5 ");
+	}
     fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
     fclose(fh);
     return TS_SUCCESS;
@@ -1086,6 +1171,7 @@
         CFG_SIMPLE_INT("nmono", &tape->nmono),
 	CFG_SIMPLE_INT("nfil",&tape->nfil),
 	CFG_SIMPLE_INT("nfono",&tape->nfono),
+	CFG_SIMPLE_INT("internal_poly",&tape->internal_poly),
 	CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
 	CFG_SIMPLE_FLOAT("R_nucleusX",&tape->R_nucleusX),
 	CFG_SIMPLE_FLOAT("R_nucleusY",&tape->R_nucleusY),
@@ -1097,6 +1183,8 @@
 	CFG_SIMPLE_INT("constvolswitch",&tape->constvolswitch),
 	CFG_SIMPLE_INT("constareaswitch",&tape->constareaswitch),
 	CFG_SIMPLE_FLOAT("constvolprecision",&tape->constvolprecision),
+	CFG_SIMPLE_INT("stretchswitch",&tape->stretchswitch),
+	CFG_SIMPLE_FLOAT("xkA0",&tape->xkA0),	
 	CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
 	CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
 	CFG_SIMPLE_FLOAT("xi",&tape->xi),
@@ -1107,18 +1195,33 @@
         CFG_SIMPLE_INT("iterations",&tape->iterations),
 	CFG_SIMPLE_INT("mcsweeps",&tape->mcsweeps),
 	CFG_SIMPLE_INT("inititer", &tape->inititer),
-        CFG_SIMPLE_BOOL("quiet",&tape->quiet),
-        CFG_SIMPLE_STR("multiprocessing",tape->multiprocessing),
+				CFG_SIMPLE_BOOL("quiet",(cfg_bool_t *)&tape->quiet),
+        CFG_SIMPLE_STR("multiprocessing",&tape->multiprocessing),
         CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
         CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
         CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
 	CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
+	CFG_SIMPLE_INT("number_of_vertices_with_c0", &tape->number_of_vertices_with_c0),
+	CFG_SIMPLE_FLOAT("c0",&tape->c0),
+	CFG_SIMPLE_FLOAT("w",&tape->w),
+	CFG_SIMPLE_FLOAT("F",&tape->F),
+/* Variables related to plane confinement */
+	CFG_INT("plane_confinement_switch", 0, CFGF_NONE),
+	CFG_FLOAT("plane_d", 15, CFGF_NONE),
+	CFG_FLOAT("plane_F", 1000, CFGF_NONE),
+/* Variables related to stretching */
+//	CFG_FLOAT("stretchswitch", 0, CFGF_NONE),
+//	CFG_FLOAT("xkA0",0,CFGF_NONE),
         CFG_END()
     };
     cfg_t *cfg;    
     ts_uint retval;
     cfg = cfg_init(opts, 0);
     retval=cfg_parse_buf(cfg, buffer);
+	tape->plane_confinement_switch=cfg_getint(cfg,"plane_confinement_switch");
+	tape->plane_d=cfg_getfloat(cfg,"plane_d");
+	tape->plane_F=cfg_getfloat(cfg,"plane_F");
+
     if(retval==CFG_FILE_ERROR){
 	fatal("No tape file.",100);
 	}
@@ -1147,6 +1250,8 @@
 ts_bool getcmdline_tape(cfg_t *cfg, char *opts){
 
 	char *commands, *backup, *saveptr, *saveopptr, *command, *operator[2];
+	operator[0]=0;
+	operator[1]=0;
 	ts_uint i,j;
 	commands=(char *)malloc(10000*sizeof(char));
     backup=commands; //since the pointer to commands will be lost, we acquire a pointer that will serve as backup.

--
Gitblit v1.9.3