From ea1cce4c80ddb68c0b46a6873ef3b64cc4a958a1 Mon Sep 17 00:00:00 2001 From: mihaf <miha.fosnaric@gmail.com> Date: Mon, 24 Mar 2014 16:37:44 +0000 Subject: [PATCH] Added dmin_interspecies. --- src/io.c | 365 +++++++++++++++++++++++++++++++++++++++------------ 1 files changed, 277 insertions(+), 88 deletions(-) diff --git a/src/io.c b/src/io.c index 93fda4c..a25b16c 100644 --- a/src/io.c +++ b/src/io.c @@ -18,7 +18,7 @@ #include <errno.h> /** DUMP STATE TO DISK DRIVE **/ -ts_bool dump_state(ts_vesicle *vesicle){ +ts_bool dump_state(ts_vesicle *vesicle, ts_uint iteration){ /* save current state with wrong pointers. Will fix that later */ ts_uint i,j,k; @@ -36,6 +36,8 @@ fwrite(vesicle->clist, sizeof(ts_cell_list),1,fh); /* dump poly list */ fwrite(vesicle->poly_list, sizeof(ts_poly_list),1,fh); + /* dump filament list */ + fwrite(vesicle->filament_list, sizeof(ts_poly_list),1,fh); /* level 1 complete */ /*dump vertices*/ @@ -118,6 +120,43 @@ } } + + /*dump filamentes grandes svinjas */ + for(i=0;i<vesicle->filament_list->n;i++){ + fwrite(vesicle->filament_list->poly[i],sizeof(ts_poly),1,fh); + fwrite(vesicle->filament_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh); + fwrite(vesicle->filament_list->poly[i]->blist,sizeof(ts_bond_list),1,fh); + } + + /* dump filamentes vertex(monomer) list*/ + for(i=0;i<vesicle->filament_list->n;i++){ + for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){ + fwrite(vesicle->filament_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh); + /* dump offset for neigh and bond */ + for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no;k++){ + // off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]-vesicle->poly_list->poly[i]->vlist->vtx[0]); + fwrite(&vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh[k]->idx,sizeof(ts_uint),1,fh); + } + for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no;k++){ + //off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]-vesicle->poly_list->poly[i]->blist->bond[0]); + fwrite(&vesicle->filament_list->poly[i]->vlist->vtx[j]->bond[k]->idx,sizeof(ts_uint),1,fh); + } + } + } + /* dump poly bonds between monomers list*/ + for(i=0;i<vesicle->filament_list->n;i++){ + for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){ + fwrite(vesicle->filament_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh); + /* dump vtx1 and vtx2 offsets */ + //off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx1-vesicle->poly_list->poly[i]->vlist->vtx[0]); + fwrite(&vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx,sizeof(ts_uint),1,fh); +// off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]); + fwrite(&vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx,sizeof(ts_uint),1,fh); + } + } + + + /* pointer offsets for fixing the restored pointers */ /* need pointers for vlist->vtx @@ -131,16 +170,31 @@ */ fwrite(vesicle->clist, sizeof(ts_cell_list),1, fh); - + + fwrite(&iteration, sizeof(ts_uint),1,fh); fclose(fh); return TS_SUCCESS; } /** RESTORE DUMP FROM DISK **/ -ts_vesicle *restore_state(){ +ts_vesicle *restore_state(ts_uint *iteration){ ts_uint i,j,k; FILE *fh=fopen("dump.bin","rb"); + + struct stat sb; + if (stat("dump.bin", &sb) == -1) { + //dump file does not exist. + return NULL; + } + + //check if it is regular file + if((sb.st_mode & S_IFMT) != S_IFREG) { + //dump file is not a regular file. + ts_fprintf(stderr,"Dump file is not a regular file!\n"); + return NULL; + } + ts_uint retval; ts_uint idx; @@ -165,6 +219,9 @@ /* restore poly list */ vesicle->poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list)); retval=fread(vesicle->poly_list, sizeof(ts_poly_list),1,fh); + /* restore filament list */ + vesicle->filament_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list)); + retval=fread(vesicle->filament_list, sizeof(ts_poly_list),1,fh); /* level 1 complete */ /* prerequisity. Bonds must be malloced before vertexes are recreated */ @@ -307,6 +364,61 @@ } } + /*restore filaments */ + vesicle->filament_list->poly = (ts_poly **)calloc(vesicle->filament_list->n,sizeof(ts_poly *)); + for(i=0;i<vesicle->filament_list->n;i++){ + vesicle->filament_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly)); + retval=fread(vesicle->filament_list->poly[i],sizeof(ts_poly),1,fh); + vesicle->filament_list->poly[i]->vlist=(ts_vertex_list *)calloc(1,sizeof(ts_vertex_list)); + retval=fread(vesicle->filament_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh); + vesicle->filament_list->poly[i]->blist=(ts_bond_list *)calloc(1,sizeof(ts_bond_list)); + retval=fread(vesicle->filament_list->poly[i]->blist,sizeof(ts_bond_list),1,fh); + /* initialize adress space for pointers that will hold specific vertices (monomers) and bonds */ + vesicle->filament_list->poly[i]->vlist->vtx=(ts_vertex **)calloc(vesicle->filament_list->poly[i]->vlist->n,sizeof(ts_vertex *)); + vesicle->filament_list->poly[i]->blist->bond=(ts_bond **)calloc(vesicle->filament_list->poly[i]->blist->n,sizeof(ts_bond *)); + for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){ + vesicle->filament_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex)); + } + for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){ + vesicle->filament_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond)); + } + + } + + + /* restore poly vertex(monomer) list*/ + for(i=0;i<vesicle->filament_list->n;i++){ + for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){ + retval=fread(vesicle->filament_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh); + + /* restore neigh and bonds */ + vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh=(ts_vertex **)calloc(vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no, sizeof(ts_vertex *)); + for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no;k++){ + retval=fread(&idx,sizeof(ts_uint),1,fh); + vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->filament_list->poly[i]->vlist->vtx[idx]; + } + vesicle->filament_list->poly[i]->vlist->vtx[j]->bond=(ts_bond **)calloc(vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no, sizeof(ts_bond *)); + for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no;k++){ + retval=fread(&idx,sizeof(ts_uint),1,fh); + vesicle->filament_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->filament_list->poly[i]->blist->bond[idx]; + } + + } + } + + /* restore poly bonds between monomers list*/ + for(i=0;i<vesicle->filament_list->n;i++){ + for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){ + // vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond)); + retval=fread(vesicle->filament_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh); + /* restore vtx1 and vtx2 */ + retval=fread(&idx,sizeof(ts_uint),1,fh); + vesicle->filament_list->poly[i]->blist->bond[j]->vtx1=vesicle->filament_list->poly[i]->vlist->vtx[idx]; + retval=fread(&idx,sizeof(ts_uint),1,fh); + vesicle->filament_list->poly[i]->blist->bond[j]->vtx2=vesicle->filament_list->poly[i]->vlist->vtx[idx]; + } + } + // recreating space for cells // vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list)); retval=fread(vesicle->clist, sizeof(ts_cell_list), 1,fh); @@ -316,6 +428,7 @@ vesicle->clist->cell[i]->idx=i+1; // We enumerate cells! Probably never required! } + retval=fread(iteration,sizeof(ts_uint),1,fh); if(retval); fclose(fh); return vesicle; @@ -324,24 +437,30 @@ ts_bool parse_args(int argc, char **argv){ - int c, retval; - DIR *dir; - path[0]=0; + int c, retval; + struct stat sb; + sprintf(command_line_args.path, "./"); //clear string; + sprintf(command_line_args.output_fullfilename,"./output.pvd"); + sprintf(command_line_args.dump_fullfilename,"./dump.bin"); + sprintf(command_line_args.tape_fullfilename,"./tape"); + FILE *file; + while (1) { static struct option long_options[] = { - {"force-from-tape", no_argument, &force_from_tape, 1}, + {"force-from-tape", no_argument, &(command_line_args.force_from_tape), 1}, + {"reset-iteration-count", no_argument, &(command_line_args.reset_iteration_count), 1}, {"tape", no_argument, 0, 't'}, {"output-file", required_argument, 0, 'o'}, {"directory", required_argument, 0, 'd'}, - {"dump-file", required_argument,0, 'f'}, + {"dump-filename", required_argument,0, 'f'}, {0, 0, 0, 0} }; /* getopt_long stores the option index here. */ int option_index = 0; - c = getopt_long (argc, argv, "d:fot", + c = getopt_long (argc, argv, "d:f:o:t:", long_options, &option_index); /* Detect the end of the options. */ @@ -362,46 +481,56 @@ case 't': //check if tape exists. If not, fail immediately. - puts ("option -t\n"); + if (stat(optarg, &sb) == -1) { + ts_fprintf(stderr,"Tape '%s' does not exist!\n",optarg); + fatal("Please select correct tape",1); + } else { + strcpy(command_line_args.tape_fullfilename,optarg); + } break; case 'o': //set filename of master output file - printf ("option -o with value `%s'\n", optarg); - break; + if ((file = fopen(optarg, "w")) == NULL) { + fprintf(stderr,"Could not create output file!\n"); + fatal("Please specify correct output file",1); + + } else { + fclose(file); + } + strcpy(command_line_args.output_fullfilename, optarg); + break; case 'd': //check if directory exists. If not create one. If creation is //successful, set directory for output files. //printf ("option -d with value `%s'\n", optarg); - dir = opendir(optarg); - if (dir) - { - /* Directory exists. */ - closedir(dir); - } - else if (ENOENT == errno) - { - /* Directory does not exist. */ + if (stat(optarg, &sb) == -1) { + //directory does not exist retval=mkdir(optarg, 0700); if(retval){ - fatal("Could not create requested directory. Check if you have permissions",1); + fatal("Could not create requested directory. Check if you have permissions",1); } } - else - { - /* opendir() failed for some other reason. */ - fatal("Could not check if directory exists. Reason unknown",1); + //check if is a proper directory + else if((sb.st_mode & S_IFMT) != S_IFDIR) { + //it is not a directory. fatal error. + ts_fprintf(stderr,"%s is not a directory!\n",optarg); + fatal("Cannot continue",1); } - ts_fprintf(stdout,"\n*** Using output directory: %s\n\n", optarg); -// sprintf(path,"%s", optarg); - strcpy(path, optarg); - // ts_fprintf(stdout,"ok!\n"); - + strcpy(command_line_args.path, optarg); break; case 'f': //check if dump file specified exists. Defaults to dump.bin + if ((file = fopen(optarg, "w")) == NULL) { + fprintf(stderr,"Could not create dump file!\n"); + fatal("Please specify correct dump file",1); + + } else { + fclose(file); + } + strcpy(command_line_args.dump_fullfilename, optarg); break; case '?': @@ -632,8 +761,8 @@ /* Here comes header of the file */ //find number of extra vtxs and bonds of polymeres - ts_uint monono=0, polyno=0, poly_idx=0; - ts_bool poly=0; + ts_uint monono=0, polyno=0, poly_idx=0, filno=0, fonono=0; + ts_bool poly=0, fil=0; if(vesicle->poly_list!=NULL){ if(vesicle->poly_list->poly[0]!=NULL){ polyno=vesicle->poly_list->n; @@ -641,8 +770,17 @@ poly=1; } } + + if(vesicle->filament_list!=NULL){ + if(vesicle->filament_list->poly[0]!=NULL){ + filno=vesicle->filament_list->n; + fonono=vesicle->filament_list->poly[0]->vlist->n; + fil=1; + } + } + fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n <UnstructuredGrid>\n"); - fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno, blist->n+monono*polyno); + 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\">"); for(i=0;i<vlist->n;i++){ fprintf(fh,"%u ",vtx[i]->idx); @@ -652,6 +790,16 @@ 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,"%u ", poly_idx); + } + } + } + //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,"%u ", poly_idx); } } @@ -666,6 +814,14 @@ for(i=0;i<vesicle->poly_list->n;i++){ for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){ 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 ); + } + } + } + //filaments + if(fil){ + for(i=0;i<vesicle->filament_list->n;i++){ + for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){ + fprintf(fh,"%e %e %e\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 ); } } } @@ -688,14 +844,26 @@ } + //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]->blist->n;j++){ + fprintf(fh,"%u %u\n", vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono); +// fprintf(stderr,"was here\n"); + + } + } + + } fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">"); - for (i=2;i<(blist->n+monono*polyno)*2+1;i+=2){ + for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2){ fprintf(fh,"%u ",i); } fprintf(fh,"\n"); fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n"); - for (i=0;i<blist->n+monono*polyno;i++){ + for (i=0;i<blist->n+monono*polyno+fonono*filno;i++){ fprintf(fh,"3 "); } @@ -750,80 +918,101 @@ -ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations){ - long int nshell=17,ncxmax=60, ncymax=60, nczmax=60, npoly=10, nmono=20, pswitch=0; // THIS IS DUE TO CONFUSE BUG! - char *buf=malloc(255*sizeof(char)); - long int brezveze0=1; +ts_bool write_pov_file(ts_vesicle *vesicle, char *filename){ + FILE *fh; + ts_uint i; + + fh=fopen(filename, "w"); + if(fh==NULL){ + err("Cannot open file %s for writing"); + return TS_FAIL; + } + + for(i=0;i<vesicle->tlist->n;i++){ + + fprintf(fh,"\ttriangle {"); + fprintf(fh,"\t<%e,%e,%e> <%e,%e,%e> <%e,%e,%e> }\n", + vesicle->tlist->tria[i]->vertex[0]->x, + vesicle->tlist->tria[i]->vertex[0]->y, + vesicle->tlist->tria[i]->vertex[0]->z, + + vesicle->tlist->tria[i]->vertex[1]->x, + vesicle->tlist->tria[i]->vertex[1]->y, + vesicle->tlist->tria[i]->vertex[1]->z, + + vesicle->tlist->tria[i]->vertex[2]->x, + vesicle->tlist->tria[i]->vertex[2]->y, + vesicle->tlist->tria[i]->vertex[2]->z + ); + } + + fclose(fh); + return TS_SUCCESS; +} + + +ts_tape *parsetape(char *filename){ + // long int nshell=17,ncxmax=60, ncymax=60, nczmax=60, npoly=10, nmono=20, pswitch=0; // THIS IS DUE TO CONFUSE BUG! + ts_tape *tape=(ts_tape *)calloc(1,sizeof(ts_tape)); + tape->multiprocessing=calloc(255,sizeof(char)); + /* long int brezveze0=1; long int brezveze1=1; long int brezveze2=1; ts_double xk0=25.0, dmax=1.67,stepsize=0.15,kspring=800.0,pressure=0.0; long int iter=1000, init=1000, mcsw=1000; - +*/ cfg_opt_t opts[] = { - CFG_SIMPLE_INT("nshell", &nshell), - CFG_SIMPLE_INT("npoly", &npoly), - CFG_SIMPLE_INT("nmono", &nmono), - CFG_SIMPLE_FLOAT("dmax", &dmax), - CFG_SIMPLE_FLOAT("xk0",&xk0), - CFG_SIMPLE_INT("pswitch",&pswitch), - CFG_SIMPLE_FLOAT("pressure",&pressure), - CFG_SIMPLE_FLOAT("k_spring",&kspring), - CFG_SIMPLE_FLOAT("stepsize",&stepsize), - CFG_SIMPLE_INT("nxmax", &ncxmax), - CFG_SIMPLE_INT("nymax", &ncymax), - CFG_SIMPLE_INT("nzmax", &nczmax), - CFG_SIMPLE_INT("iterations",&iter), - CFG_SIMPLE_INT("mcsweeps",&mcsw), - CFG_SIMPLE_INT("inititer", &init), - CFG_SIMPLE_BOOL("quiet",&quiet), - CFG_SIMPLE_STR("multiprocessing",buf), - CFG_SIMPLE_INT("smp_cores",&brezveze0), - CFG_SIMPLE_INT("cluster_nodes",&brezveze1), - CFG_SIMPLE_INT("distributed_processes",&brezveze2), + CFG_SIMPLE_INT("nshell", &tape->nshell), + CFG_SIMPLE_INT("npoly", &tape->npoly), + CFG_SIMPLE_INT("nmono", &tape->nmono), + CFG_SIMPLE_INT("nfil",&tape->nfil), + CFG_SIMPLE_INT("nfono",&tape->nfono), + CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus), + CFG_SIMPLE_FLOAT("dmax", &tape->dmax), + CFG_SIMPLE_FLOAT("dmin_interspecies", &tape->dmin_interspecies), + CFG_SIMPLE_FLOAT("xk0",&tape->xk0), + CFG_SIMPLE_INT("pswitch",&tape->pswitch), + CFG_SIMPLE_FLOAT("pressure",&tape->pressure), + CFG_SIMPLE_FLOAT("k_spring",&tape->kspring), + CFG_SIMPLE_FLOAT("xi",&tape->xi), + CFG_SIMPLE_FLOAT("stepsize",&tape->stepsize), + CFG_SIMPLE_INT("nxmax", &tape->ncxmax), + CFG_SIMPLE_INT("nymax", &tape->ncymax), + CFG_SIMPLE_INT("nzmax", &tape->nczmax), + 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_INT("smp_cores",&tape->brezveze0), + CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1), + CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2), CFG_END() }; cfg_t *cfg; ts_uint retval; cfg = cfg_init(opts, 0); - retval=cfg_parse(cfg, "tape"); + retval=cfg_parse(cfg, filename); if(retval==CFG_FILE_ERROR){ fatal("No tape file.",100); } else if(retval==CFG_PARSE_ERROR){ fatal("Invalid tape!",100); } - ts_vesicle *vesicle; - *iterations=iter; - *inititer=init; - *mcsweeps=mcsw; - vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize); - vesicle->poly_list=init_poly_list(npoly,nmono, vesicle->vlist); - vesicle->spring_constant=kspring; - poly_assign_spring_const(vesicle); - - vesicle->nshell=nshell; - vesicle->dmax=dmax*dmax; - vesicle->bending_rigidity=xk0; - vtx_set_global_values(vesicle); //copies xk0 to every vertex - - - vesicle->stepsize=stepsize; - vesicle->clist->ncmax[0]=ncxmax; - vesicle->clist->ncmax[1]=ncymax; - vesicle->clist->ncmax[2]=nczmax; - vesicle->clist->max_occupancy=8; - vesicle->pressure=pressure/vesicle->bending_rigidity; //all energy contributions need to be divided by bending_rigidity! - vesicle->pswitch=pswitch; cfg_free(cfg); - free(buf); - // fprintf(stderr,"NSHELL=%u\n",vesicle->nshell); - - return vesicle; + /* global variables are set automatically */ + quiet=tape->quiet; + return tape; +} +ts_bool tape_free(ts_tape *tape){ + free(tape->multiprocessing); + free(tape); + return TS_SUCCESS; } -- Gitblit v1.9.3