From f977e8fe32b6a1a2ae85662ab35fd7c9f647cc3e Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@gmail.com> Date: Mon, 02 Jul 2018 10:50:58 +0000 Subject: [PATCH] Added a funciton that calculates everything needed for spherical harmonics --- src/io.c | 154 +++++++++++++++++++++++++++++++++++++++++++++++---- 1 files changed, 142 insertions(+), 12 deletions(-) diff --git a/src/io.c b/src/io.c index 7cb0d62..3ea3662 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); @@ -612,11 +613,17 @@ fprintf(fd,"Flags:\n\n"); 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"); 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"); - fprintf(fd,"--reset-iteration-count\t\t starts simulation from the beginning (using binary dump or tape).\n"); + fprintf(fd,"--reset-iteration-count\t\t starts simulation from the beginning (using binary dump).\n"); fprintf(fd,"--tape (or -t)\t\t specifies tape filename. For --force-from-tape and restoring from binary dump. Defaults to 'tape'.\n"); fprintf(fd,"--version (or -v)\t\t Prints version information.\n"); fprintf(fd,"--output-file (or -o)\t\t Specifies filename of .PVD file. Defaults to 'output.pvd'\n"); - fprintf(fd,"--dump-filename (or -f)\t\t specifies filename for binary dump&restore. Defaults to 'dump.bin'\n\n"); + fprintf(fd,"--dump-filename (or -f)\t\t specifies filename for binary dump&restore. Defaults to 'dump.bin'\n\n\n"); + fprintf(fd,"Examples:\n\n"); + fprintf(fd,"trisurf --force-from-tape\n"); + fprintf(fd,"trisurf --reset-iteration-count\n"); + fprintf(fd,"trisurf --restore-from-vtk filename.vtu\n"); + fprintf(fd,"\n\n"); + return TS_SUCCESS; } @@ -814,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; @@ -855,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); } @@ -880,7 +888,113 @@ } } - fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n"); + 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++){ + fprintf(fh,"%.17e ",vtx[i]->c); + } + //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]->c); + } + } + } + //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]->c); + } + } + } + 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"); + + 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); } @@ -931,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; @@ -1050,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), @@ -1061,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), @@ -1071,12 +1195,16 @@ 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), CFG_END() }; cfg_t *cfg; @@ -1111,6 +1239,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