Trisurf Monte Carlo simulator
Samo Penic
2019-10-18 842378276ab4c4a8dabd02324bb24b8da39c4384
Done and implemented compression. It doesn't segfaults at the moment
3 files modified
232 ■■■■ changed files
src/io.c 212 ●●●● patch | view | raw | blame | history
src/snapshot.c 18 ●●●● patch | view | raw | blame | history
src/snapshot.h 2 ●●● patch | view | raw | blame | history
src/io.c
@@ -256,6 +256,9 @@
    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=\"binary\">");
    long *int_vector=(long *)malloc((vlist->n+monono*polyno+fonono*filno)*sizeof(long));
    ts_double *float_vector=(double *)malloc((vlist->n+monono*polyno+fonono*filno)*sizeof(double));
    ts_double *float_cell_vector=(double *)malloc((blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n)*sizeof(double));
    ts_double *coordinate_vector=(double *)malloc(3*(vlist->n+monono*polyno+fonono*filno)*sizeof(double));
// sizes are as expected 4 bytes for int and 1 byte for char
//    printf("size of ts_uint %ld, of int %ld, of char %ld",sizeof(ts_uint), sizeof(int), sizeof(char));
    int offset=0;
@@ -288,26 +291,30 @@
        }
    }
    printf("Offset in the end is %d, should be %d",offset,(vlist->n+monono*polyno+fonono*filno) );
    char *printout=ts_compress((char *)int_vector,offset*sizeof(long), offset);
    char *printout=ts_compress((char *)int_vector,offset*sizeof(long));
    fprintf(fh,"%s",printout);
    free(printout);
        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++){
        fprintf(fh,"<DataArray type=\"Int64\" Name=\"vertices_in_cluster\" format=\"binary\">");
        offset=0;
        for(i=0;i<vlist->n;i++,offset++){
            if(vtx[i]->cluster!=NULL){
                fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
//                fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
                int_vector[offset]=vtx[i]->cluster->nvtx;
            } else {
                fprintf(fh,"-1 ");
//                fprintf(fh,"-1 ");
                int_vector[offset]=-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 ");
                for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
                    //fprintf(fh,"-1 ");
                    int_vector[offset]=-1;
                }
            }
        }
@@ -315,12 +322,15 @@
        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 ");
                for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
    //                    fprintf(fh,"-1 ");
                    int_vector[offset]=-1;
                }
            }
        }
        char *printout=ts_compress((char *)int_vector,offset*sizeof(long));
        fprintf(fh,"%s",printout);
        free(printout);
        fprintf(fh,"</DataArray>\n");
@@ -328,16 +338,19 @@
    }
    //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);
    fprintf(fh,"<DataArray type=\"Float64\" Name=\"spontaneous_curvature\" format=\"binary\">");
    offset=0;
    for(i=0;i<vlist->n;i++, offset++){
//        fprintf(fh,"%.17e ",vtx[i]->c);
        float_vector[offset]=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);
                for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++, offset++){
                    //fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->c);
                    float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->c;
                }
            }
        }
@@ -345,25 +358,34 @@
        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++){
                for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
        //    fprintf(stderr,"was here\n");
                    fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->c);
                    //fprintf(fh,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->c);
                    float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->c;
                }
            }
        }
        printout=ts_compress((char *)float_vector,offset*sizeof(double));
        fprintf(fh,"%s",printout);
        free(printout);
    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);
    fprintf(fh,"<DataArray type=\"Float64\" Name=\"bending_energy\" format=\"binary\">");
    offset=0;
    for(i=0;i<vlist->n;i++,offset++){
//        fprintf(fh,"%.17e ",vtx[i]->energy*vtx[i]->xk);
        float_vector[offset]=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);
                for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
                    //fprintf(fh,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k);
                    float_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k;
                }
            }
        }
@@ -371,68 +393,108 @@
        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++){
                for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++,offset++){
        //    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,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k);
                    float_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k;
                }
            }
        }
    printout=ts_compress((char *)float_vector,offset*sizeof(double));
        fprintf(fh,"%s",printout);
        free(printout);
    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 ");
        fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"binary\">");
        offset=0;
        for(i=0;i<blist->n;i++,offset++){
            //fprintf(fh, "0.0 ");
            float_cell_vector[offset]=0.0;
        }
        for(i=0;i<monono*polyno+filno*(fonono-1);i++){
            fprintf(fh,"0.0 ");
        for(i=0;i<monono*polyno+filno*(fonono-1);i++,offset++){
            float_cell_vector[offset]=0.0;
            //fprintf(fh,"0.0 ");
        }
        for(i=0;i<vesicle->tlist->n;i++){
            fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
        for(i=0;i<vesicle->tlist->n;i++,offset++){
            float_cell_vector[offset]=vesicle->tlist->tria[i]->energy;
            //fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
        }
        printout=ts_compress((char *)float_cell_vector,offset*sizeof(double));
        fprintf(fh,"%s",printout);
        free(printout);
        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);
    fprintf(fh,"</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"binary\">\n");
    offset=0;
    for(i=0;i<vlist->n;i++,offset+=3){
        //fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
        coordinate_vector[offset]=vtx[i]->x;
        coordinate_vector[offset+1]=vtx[i]->y;
        coordinate_vector[offset+2]=vtx[i]->z;
    }
    //polymeres
    if(poly){
        for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
                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 );
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,offset+=3){
                //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 );
        coordinate_vector[offset]=vesicle->poly_list->poly[i]->vlist->vtx[j]->x;
        coordinate_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[j]->y;
        coordinate_vector[offset+2]=vesicle->poly_list->poly[i]->vlist->vtx[j]->z;
            }
        }
    }
    //filaments
    if(fil){
        for(i=0;i<vesicle->filament_list->n;i++){
        for(i=0;i<vesicle->filament_list->n;i++,offset+=3){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
                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 );
    //        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 );
        coordinate_vector[offset]=vesicle->filament_list->poly[i]->vlist->vtx[j]->x;
        coordinate_vector[offset+1]=vesicle->filament_list->poly[i]->vlist->vtx[j]->y;
        coordinate_vector[offset+2]=vesicle->filament_list->poly[i]->vlist->vtx[j]->z;
            }
        }
    }
    printout=ts_compress((char *)coordinate_vector,offset*sizeof(double));
        fprintf(fh,"%s",printout);
        free(printout);
    fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
    for(i=0;i<blist->n;i++){
            fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
    fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"binary\">");
    free(int_vector);
    int_vector=(long *)malloc(((blist->n+monono*polyno+filno*(fonono-1))*2+vesicle->tlist->n*3)*sizeof(long));
    offset=0;
    for(i=0;i<blist->n;i++,offset+=2){
//            fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
            int_vector[offset]=blist->bond[i]->vtx1->idx;
            int_vector[offset+1]=blist->bond[i]->vtx2->idx;
    }
    //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]->blist->n;j++){
//                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);
                fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono);
            for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++,offset+=2){
                //fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono);
                int_vector[offset]= vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+i*monono;
                int_vector[offset+1]=vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+i*monono;
            }
    //grafted bonds
        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);
        //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);
        int_vector[offset]=vesicle->poly_list->poly[i]->grafted_vtx->idx;
        int_vector[offset+1]=vesicle->poly_list->poly[i]->vlist->vtx[0]->idx+vlist->n+i*monono;
        offset+=2;
        }
    }
@@ -441,32 +503,62 @@
    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);
            for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++,offset+=2){
//                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);
                int_vector[offset]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono;
                int_vector[offset+1]=vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono;
//        fprintf(stderr,"was here\n");
            
            }
        }
    }
    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);
    for(i=0;i<vesicle->tlist->n;i++,offset+=3){
//        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);
        int_vector[offset]=vesicle->tlist->tria[i]->vertex[0]->idx;
        int_vector[offset+1]=vesicle->tlist->tria[i]->vertex[1]->idx;
        int_vector[offset+2]=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);
        printout=ts_compress((char *)int_vector,offset*sizeof(long));
        fprintf(fh,"%s",printout);
        free(printout);
    fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"binary\">");
    offset=0;
    for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2,offset++){
//    fprintf(fh,"%u ",i);
    int_vector[offset]=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);
    for(j=i+1;j<i+3*(vesicle->tlist->n);j+=3,offset++){ //let's continue counting from where we left of
//        fprintf(fh,"%u ", j);
    int_vector[offset]=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-1)*filno;i++){
        fprintf(fh,"3 ");
//    fprintf(fh,"\n");
    printout=ts_compress((char *)int_vector,offset*sizeof(long));
    fprintf(fh,"%s",printout);
    free(printout);
    fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"binary\">\n");
    free(int_vector);
    free(float_vector);
    free(float_cell_vector);
    free(coordinate_vector);
    unsigned char *char_vector=(unsigned char *)malloc((blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n)*sizeof(char));
    offset=0;
     for (i=0;i<blist->n+monono*polyno+(fonono-1)*filno;i++,offset++){
        //fprintf(fh,"3 ");
    char_vector[offset]=3;
    }
    for(i=0;i<vesicle->tlist->n;i++){
        fprintf(fh,"5 ");
    for(i=0;i<vesicle->tlist->n;i++,offset++){
    //    fprintf(fh,"5 ");
        char_vector[offset]=5;
    }
    printout=ts_compress((char *)char_vector,offset*sizeof(char));
    fprintf(fh,"%s",printout);
    free(printout);
    free(char_vector);
    fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
    fclose(fh);
    return TS_SUCCESS;
src/snapshot.c
@@ -173,19 +173,27 @@
    return nbase;
}
char *ts_compress(char *data, ts_uint data_len, ts_uint original_len){
char *ts_compress(char *data, ts_uint data_len){
    size_t nbase1, nbase2;
    char *compr;
    size_t number_of_compressed_bytes=ts_compress_data(data, data_len, &compr);
    unsigned char *compr=(unsigned char *)malloc(data_len);
    size_t number_of_compressed_bytes=data_len;
    compress(compr,&number_of_compressed_bytes, (unsigned char *)data, data_len);
//    printf("Compressdion error code=%d, Z_OK=%d, Z_BUF_ERROR=%d", errcode, Z_OK, Z_BUF_ERROR);
    char *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase1);
    free(compr);
    ts_uint header[4]={1, original_len, original_len, number_of_compressed_bytes};
    ts_uint header[4]={1, data_len, data_len, number_of_compressed_bytes};
    char *encoded_header=(char *)base64_encode((unsigned char *)header, 4*sizeof(ts_uint), &nbase2);
    char *return_value=malloc((nbase1+nbase2+1)*sizeof(char));
    strncpy(return_value,encoded_header,nbase2);
    strncpy(return_value+nbase2,encoded_compressed,nbase1);
    *(return_value+nbase1+nbase2)=0;
    printf("Compressed size in bytes= %ld, size of encoded header= %ld, size of encoded compressed= %ld.\n",number_of_compressed_bytes, nbase2, nbase1);
//    printf("Compressed size in bytes= %ld, size of encoded header= %ld, size of encoded compressed= %ld.\n",number_of_compressed_bytes, nbase2, nbase1);
    free(encoded_compressed);
    free(encoded_header);
    return return_value;
src/snapshot.h
@@ -26,5 +26,5 @@
void build_decoding_table();
void base64_cleanup();
ts_uint ts_compress_string64(char *data, ts_uint data_len, char **compressed);
char *ts_compress(char *data, ts_uint data_len, ts_uint original_len);
char *ts_compress(char *data, ts_uint data_len);
#endif