Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2019-10-18 a69203a95d66595b80891aafc1ca8a59303290d0
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
SP 2 /* vim: set ts=4 sts=4 sw=4 noet : */    
854cb6 3 #include<stdio.h>
SP 4 #include<general.h>
5 #include<snapshot.h>
22721d 6 #include<stdlib.h>
SP 7 #include<stdarg.h>
8 #include <zlib.h>
f0bcea 9 #include<inttypes.h>
11fef8 10 #include<config.h>
2e8c7f 11 #include <time.h>
aecb11 12 #include <string.h>
698ae1 13 #include "io.h"
a3dbf0 14 /* a helper function that utilizes ts_string data structure and performs same as sprintf */
22721d 15 ts_uint ts_sprintf(ts_string *str, char *fmt, ...){
SP 16     va_list ap;
17     va_start(ap,fmt);
18     ts_uint n=vsprintf(&(str->string[str->beg]),fmt,ap);
19     va_end(ap);
20     str->beg+=n;
21     return n;
22 }
23
a3dbf0 24 /* outputs additional data into paraview xml file */
854cb6 25 ts_bool xml_trisurf_data(FILE *fh, ts_vesicle *vesicle){
22721d 26
SP 27     ts_string *data=(ts_string *)malloc(sizeof(ts_sprintf));
02d65c 28     data->string=(char *)malloc(5120000*sizeof(char)); /*TODO: warning, can break if the string is to long */
22721d 29     data->beg=0;
SP 30     
854cb6 31     xml_trisurf_header(fh, vesicle);
22721d 32     xml_trisurf_tria(data,vesicle->tlist);
SP 33     xml_trisurf_tria_neigh(data,vesicle->tlist);
34     xml_trisurf_vtx_neigh(data,vesicle->vlist);    
35     xml_trisurf_vtx_tristar(data,vesicle->vlist);
4891eb 36     xml_trisurf_nucleus(data,vesicle);
49981c 37     xml_trisurf_constvolarea(data,V0,A0);
aecb11 38 //#ifdef COMPRESSION
f0bcea 39     char *compressed;
588bbb 40     ts_uint nbytes=ts_compress_string64(data->string, data->beg-1, &compressed); //suppress null character at the end with by substracting 1
f0bcea 41     fwrite(compressed, sizeof(unsigned char), nbytes, fh);
SP 42     free (compressed);
aecb11 43 //#else
SP 44 //    fprintf(fh,"%s", data->string);
45 //#endif
a3dbf0 46     free(data->string);  /* TODO: valgrind is not ok with this! */
22721d 47     free(data);
854cb6 48     xml_trisurf_footer(fh);
SP 49     return TS_SUCCESS;
50 }
51
52 ts_bool xml_trisurf_header(FILE *fh, ts_vesicle *vesicle){
2e8c7f 53 /* format current time */
SP 54     time_t current_time;
55         char *c_time_string;
56     current_time = time(NULL);
57     c_time_string = ctime(&current_time);
a3dbf0 58     int npoly, nfono;    
SP 59
e29756 60     fprintf(fh, "<trisurfversion>Trisurf (commit %s), compiled on %s %s</trisurfversion>\n",TS_VERSION, __DATE__,  __TIME__);
2e8c7f 61     fprintf(fh, "<dumpdate>%s</dumpdate>\n", c_time_string);
SP 62
cc95bd 63     fprintf(fh, "<tape>");
698ae1 64         fprintf(fh,"%s",tapetxt);    
f0bcea 65     fprintf(fh, "</tape>\n");
a3dbf0 66     if(vesicle->poly_list!=NULL){
SP 67         npoly=vesicle->poly_list->n;
68         if(npoly!=0){
69             nfono=vesicle->poly_list->poly[0]->vlist->n;
70         } else {
71             nfono=0;
72         }
73     } else {
74         npoly=0;
75         nfono=0;
76     }
86b69b 77     fprintf(fh, "<trisurf nvtx=\"%u\" npoly=\"%u\" nmono=\"%u\" compressed=\"false\">\n", vesicle->vlist->n, npoly, nfono);
854cb6 78     return TS_SUCCESS;
SP 79 }
80
81 ts_bool xml_trisurf_footer(FILE *fh){
82     fprintf(fh, "</trisurf>\n");
83     return TS_SUCCESS;
84 }
85
22721d 86 ts_bool xml_trisurf_tria(ts_string *data, ts_triangle_list *tlist){
854cb6 87     ts_uint i;
588bbb 88     ts_sprintf(data,"<tria>");
854cb6 89     for(i=0; i<tlist->n;i++){
e87e4e 90         ts_sprintf(data,"%u %u %u ",tlist->tria[i]->vertex[0]->idx, tlist->tria[i]->vertex[1]->idx, tlist->tria[i]->vertex[2]->idx);
854cb6 91     }
588bbb 92     ts_sprintf(data,"</tria>");
854cb6 93     return TS_SUCCESS;
SP 94 }
95
22721d 96 ts_bool xml_trisurf_tria_neigh(ts_string *data, ts_triangle_list *tlist){
854cb6 97     ts_uint i;
22721d 98     ts_sprintf(data,"<trianeigh>\n");
854cb6 99     for(i=0; i<tlist->n;i++){
e87e4e 100         ts_sprintf(data,"%u %u %u ",tlist->tria[i]->neigh[0]->idx, tlist->tria[i]->neigh[1]->idx, tlist->tria[i]->neigh[2]->idx);
854cb6 101     }
22721d 102     ts_sprintf(data,"</trianeigh>\n");
854cb6 103     return TS_SUCCESS;
SP 104 }
105
22721d 106 ts_bool xml_trisurf_vtx_neigh(ts_string *data, ts_vertex_list *vlist){
854cb6 107     ts_uint i,j;
SP 108     for(i=0;i<vlist->n;i++){
22721d 109         ts_sprintf(data,"<vtxn idx=\"%u\">",vlist->vtx[i]->idx);
854cb6 110         for(j=0;j<vlist->vtx[i]->neigh_no;j++){
22721d 111             ts_sprintf(data,"%u ",vlist->vtx[i]->neigh[j]->idx);
854cb6 112         }
588bbb 113         ts_sprintf(data, "</vtxn>");
854cb6 114     }
SP 115     return TS_SUCCESS;
116 }
117
22721d 118 ts_bool xml_trisurf_vtx_tristar(ts_string *data, ts_vertex_list *vlist){
854cb6 119     ts_uint i,j;
SP 120     for(i=0;i<vlist->n;i++){
22721d 121         ts_sprintf(data,"<tristar idx=\"%u\">",vlist->vtx[i]->idx);
854cb6 122         for(j=0;j<vlist->vtx[i]->tristar_no;j++){
22721d 123             ts_sprintf(data,"%u ",vlist->vtx[i]->tristar[j]->idx);
854cb6 124         }
588bbb 125         ts_sprintf(data, "</tristar>");
854cb6 126     }
SP 127     return TS_SUCCESS;
128 }
129
4891eb 130 ts_bool xml_trisurf_nucleus(ts_string *data, ts_vesicle* vesicle){
SP 131     if(vesicle->R_nucleus>0.0 || (vesicle->R_nucleusX>0.0 && vesicle->R_nucleusY>0.0 && vesicle->R_nucleusZ>0.0)){
49981c 132         ts_sprintf(data,"<nucleus>%.17e %.17e %.17e</nucleus>",vesicle->nucleus_center[0], vesicle->nucleus_center[1], vesicle->nucleus_center[2]);
4891eb 133     }
SP 134     return TS_SUCCESS;
135 }
49981c 136 ts_bool xml_trisurf_constvolarea(ts_string *data, ts_double volume, ts_double area){
SP 137     ts_sprintf(data,"<constant_volume>%.17e</constant_volume>",volume);
138     ts_sprintf(data,"<constant_area>%.17e</constant_area>",area);
139
140     return TS_SUCCESS;
141 }
f0bcea 142
SP 143
144 /* UTILITIES */
145
146 /* zlib compression base64 encoded */
147 /* compressed must not be pre-malloced */
8c1bb1 148 /* taken from https://gist.github.com/arq5x/5315739 */
aecb11 149 ts_uint ts_compress_data(char *data, ts_uint data_len, char **compressed){
f0bcea 150     z_stream defstream;
SP 151     defstream.zalloc = Z_NULL;
152     defstream.zfree = Z_NULL;
153     defstream.opaque = Z_NULL;
154     defstream.avail_in = data_len+1;
155     defstream.next_in = (unsigned char *)data;    
7d1025 156     char *compr=(char *)malloc(data_len*sizeof(char));
f0bcea 157     defstream.avail_out = data_len+1;
7d1025 158     defstream.next_out = (unsigned char *)compr;
a69203 159     deflateInit(&defstream, 6);
SP 160 //    deflateInit(&defstream, Z_BEST_COMPRESSION);
f0bcea 161         deflate(&defstream, Z_FINISH);
SP 162         deflateEnd(&defstream);
7d1025 163     *compressed=compr;
aecb11 164     return defstream.total_out;
SP 165 }
166
167 ts_uint ts_compress_string64(char *data, ts_uint data_len, char **encoded_compressed){
168     size_t nbase;
169     char *compr;
170     size_t number_of_compressed_bytes=ts_compress_data(data, data_len, &compr);
171     *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase);
f0bcea 172     free(compr);
SP 173     return nbase;
174 }
175
a69203 176 char *ts_compress(char *data, ts_uint data_len, ts_uint original_len){
692fb5 177     size_t nbase1, nbase2;
aecb11 178     char *compr;
a69203 179     size_t number_of_compressed_bytes=ts_compress_data(data, data_len, &compr);
692fb5 180     char *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase1);
aecb11 181     free(compr);
a69203 182     ts_uint header[4]={1, original_len, original_len, number_of_compressed_bytes};
692fb5 183     char *encoded_header=(char *)base64_encode((unsigned char *)header, 4*sizeof(ts_uint), &nbase2);
a69203 184     char *return_value=malloc((nbase1+nbase2+1)*sizeof(char));
SP 185     strncpy(return_value,encoded_header,nbase2);
186     strncpy(return_value+nbase2,encoded_compressed,nbase1);
187     *(return_value+nbase1+nbase2)=0;
188     printf("Compressed size in bytes= %ld, size of encoded header= %ld, size of encoded compressed= %ld.\n",number_of_compressed_bytes, nbase2, nbase1);
aecb11 189     free(encoded_compressed);
a69203 190     free(encoded_header);
SP 191     return return_value;
aecb11 192 }
SP 193
8c1bb1 194 ts_uint ts_decompress_string64(char *b64, ts_uint data_len, char **decompressed){
SP 195 return TS_SUCCESS;
196
197 }
f0bcea 198
SP 199 /* base64 encoding, taken from http://stackoverflow.com/questions/342409/how-do-i-base64-encode-decode-in-c */
200 static char encoding_table[] = {'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H',
201                                 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P',
202                                 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X',
203                                 'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f',
204                                 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n',
205                                 'o', 'p', 'q', 'r', 's', 't', 'u', 'v',
206                                 'w', 'x', 'y', 'z', '0', '1', '2', '3',
207                                 '4', '5', '6', '7', '8', '9', '+', '/'};
208 static char *decoding_table = NULL;
209 static int mod_table[] = {0, 2, 1};
210
211
7f6076 212 char *base64_encode(const unsigned char *data, size_t input_length, size_t *output_length) {
f0bcea 213     *output_length = 4 * ((input_length + 2) / 3);
SP 214     int i,j;
215     char *encoded_data = malloc(*output_length);
216     if (encoded_data == NULL) return NULL;
217     
218     for (i = 0, j = 0; i < input_length;) {
219
7f6076 220            uint32_t octet_a = i < input_length ? (unsigned char)data[i++] : 0;
SP 221         uint32_t octet_b = i < input_length ? (unsigned char)data[i++] : 0;
222            uint32_t octet_c = i < input_length ? (unsigned char)data[i++] : 0;
223         uint32_t triple = (octet_a << 0x10) + (octet_b << 0x08) + octet_c;
f0bcea 224
7f6076 225            encoded_data[j++] = encoding_table[(triple >> 3 * 6) & 0x3F];
SP 226            encoded_data[j++] = encoding_table[(triple >> 2 * 6) & 0x3F];
227            encoded_data[j++] = encoding_table[(triple >> 1 * 6) & 0x3F];
228            encoded_data[j++] = encoding_table[(triple >> 0 * 6) & 0x3F];
f0bcea 229     }
SP 230
231     for (i = 0; i < mod_table[input_length % 3]; i++)
7f6076 232            encoded_data[*output_length - 1 - i] = '=';
f0bcea 233
SP 234     return encoded_data;
235 }
236
237
7f6076 238 unsigned char *base64_decode(const char *data, size_t input_length, size_t *output_length) {
f0bcea 239     int i,j;
SP 240     if (decoding_table == NULL) build_decoding_table();
241
242     if (input_length % 4 != 0) return NULL;
243
244     *output_length = input_length / 4 * 3;
245     if (data[input_length - 1] == '=') (*output_length)--;
246     if (data[input_length - 2] == '=') (*output_length)--;
247
248     unsigned char *decoded_data = malloc(*output_length);
249     if (decoded_data == NULL) return NULL;
250
251     for (i = 0, j = 0; i < input_length;) {
7f6076 252            uint32_t sextet_a = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
SP 253            uint32_t sextet_b = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
254            uint32_t sextet_c = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
255            uint32_t sextet_d = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
f0bcea 256
7f6076 257            uint32_t triple = (sextet_a << 3 * 6)
f0bcea 258             + (sextet_b << 2 * 6)
SP 259             + (sextet_c << 1 * 6)
260             + (sextet_d << 0 * 6);
261
7f6076 262            if (j < *output_length) decoded_data[j++] = (triple >> 2 * 8) & 0xFF;
SP 263            if (j < *output_length) decoded_data[j++] = (triple >> 1 * 8) & 0xFF;
264            if (j < *output_length) decoded_data[j++] = (triple >> 0 * 8) & 0xFF;
f0bcea 265     }
8c1bb1 266     if(decoding_table !=NULL) free(decoding_table);
f0bcea 267     return decoded_data;
SP 268 }
269
270
271 void build_decoding_table() {
272
273     decoding_table = malloc(256);
274     int i;
275     for (i = 0; i < 64; i++)
7f6076 276         decoding_table[(unsigned char) encoding_table[i]] = i;
f0bcea 277 }
SP 278
279
280 void base64_cleanup() {
281     free(decoding_table);
282 }
283
284
285
286
287
288
289