Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2019-10-18 842378276ab4c4a8dabd02324bb24b8da39c4384
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
842378 176 char *ts_compress(char *data, ts_uint data_len){
692fb5 177     size_t nbase1, nbase2;
842378 178     unsigned char *compr=(unsigned char *)malloc(data_len);
SP 179     size_t number_of_compressed_bytes=data_len;
180
181     compress(compr,&number_of_compressed_bytes, (unsigned char *)data, data_len);
182 //    printf("Compressdion error code=%d, Z_OK=%d, Z_BUF_ERROR=%d", errcode, Z_OK, Z_BUF_ERROR);
183
692fb5 184     char *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase1);
842378 185
aecb11 186     free(compr);
842378 187
SP 188     ts_uint header[4]={1, data_len, data_len, number_of_compressed_bytes};
189
692fb5 190     char *encoded_header=(char *)base64_encode((unsigned char *)header, 4*sizeof(ts_uint), &nbase2);
a69203 191     char *return_value=malloc((nbase1+nbase2+1)*sizeof(char));
SP 192     strncpy(return_value,encoded_header,nbase2);
193     strncpy(return_value+nbase2,encoded_compressed,nbase1);
194     *(return_value+nbase1+nbase2)=0;
842378 195
SP 196 //    printf("Compressed size in bytes= %ld, size of encoded header= %ld, size of encoded compressed= %ld.\n",number_of_compressed_bytes, nbase2, nbase1);
aecb11 197     free(encoded_compressed);
a69203 198     free(encoded_header);
SP 199     return return_value;
aecb11 200 }
SP 201
8c1bb1 202 ts_uint ts_decompress_string64(char *b64, ts_uint data_len, char **decompressed){
SP 203 return TS_SUCCESS;
204
205 }
f0bcea 206
SP 207 /* base64 encoding, taken from http://stackoverflow.com/questions/342409/how-do-i-base64-encode-decode-in-c */
208 static char encoding_table[] = {'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H',
209                                 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P',
210                                 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X',
211                                 'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f',
212                                 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n',
213                                 'o', 'p', 'q', 'r', 's', 't', 'u', 'v',
214                                 'w', 'x', 'y', 'z', '0', '1', '2', '3',
215                                 '4', '5', '6', '7', '8', '9', '+', '/'};
216 static char *decoding_table = NULL;
217 static int mod_table[] = {0, 2, 1};
218
219
7f6076 220 char *base64_encode(const unsigned char *data, size_t input_length, size_t *output_length) {
f0bcea 221     *output_length = 4 * ((input_length + 2) / 3);
SP 222     int i,j;
223     char *encoded_data = malloc(*output_length);
224     if (encoded_data == NULL) return NULL;
225     
226     for (i = 0, j = 0; i < input_length;) {
227
7f6076 228            uint32_t octet_a = i < input_length ? (unsigned char)data[i++] : 0;
SP 229         uint32_t octet_b = i < input_length ? (unsigned char)data[i++] : 0;
230            uint32_t octet_c = i < input_length ? (unsigned char)data[i++] : 0;
231         uint32_t triple = (octet_a << 0x10) + (octet_b << 0x08) + octet_c;
f0bcea 232
7f6076 233            encoded_data[j++] = encoding_table[(triple >> 3 * 6) & 0x3F];
SP 234            encoded_data[j++] = encoding_table[(triple >> 2 * 6) & 0x3F];
235            encoded_data[j++] = encoding_table[(triple >> 1 * 6) & 0x3F];
236            encoded_data[j++] = encoding_table[(triple >> 0 * 6) & 0x3F];
f0bcea 237     }
SP 238
239     for (i = 0; i < mod_table[input_length % 3]; i++)
7f6076 240            encoded_data[*output_length - 1 - i] = '=';
f0bcea 241
SP 242     return encoded_data;
243 }
244
245
7f6076 246 unsigned char *base64_decode(const char *data, size_t input_length, size_t *output_length) {
f0bcea 247     int i,j;
SP 248     if (decoding_table == NULL) build_decoding_table();
249
250     if (input_length % 4 != 0) return NULL;
251
252     *output_length = input_length / 4 * 3;
253     if (data[input_length - 1] == '=') (*output_length)--;
254     if (data[input_length - 2] == '=') (*output_length)--;
255
256     unsigned char *decoded_data = malloc(*output_length);
257     if (decoded_data == NULL) return NULL;
258
259     for (i = 0, j = 0; i < input_length;) {
7f6076 260            uint32_t sextet_a = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
SP 261            uint32_t sextet_b = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
262            uint32_t sextet_c = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
263            uint32_t sextet_d = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
f0bcea 264
7f6076 265            uint32_t triple = (sextet_a << 3 * 6)
f0bcea 266             + (sextet_b << 2 * 6)
SP 267             + (sextet_c << 1 * 6)
268             + (sextet_d << 0 * 6);
269
7f6076 270            if (j < *output_length) decoded_data[j++] = (triple >> 2 * 8) & 0xFF;
SP 271            if (j < *output_length) decoded_data[j++] = (triple >> 1 * 8) & 0xFF;
272            if (j < *output_length) decoded_data[j++] = (triple >> 0 * 8) & 0xFF;
f0bcea 273     }
8c1bb1 274     if(decoding_table !=NULL) free(decoding_table);
f0bcea 275     return decoded_data;
SP 276 }
277
278
279 void build_decoding_table() {
280
281     decoding_table = malloc(256);
282     int i;
283     for (i = 0; i < 64; i++)
7f6076 284         decoding_table[(unsigned char) encoding_table[i]] = i;
f0bcea 285 }
SP 286
287
288 void base64_cleanup() {
289     free(decoding_table);
290 }
291
292
293
294
295
296
297