Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2019-10-17 aecb11873dde5425c1acdf8791736461f3628b14
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;    
156     char *compr=(char *)malloc(data_len*sizeof(char *));
157     defstream.avail_out = data_len+1;
158     defstream.next_out = (unsigned char *)compr;
159     deflateInit(&defstream, Z_BEST_COMPRESSION);
160         deflate(&defstream, Z_FINISH);
161         deflateEnd(&defstream);
aecb11 162     *compressed=compr;
SP 163 //    *compressed=base64_encode((unsigned char *)compr,(size_t)defstream.total_out,&nbase);
f0bcea 164     //fwrite(base64, sizeof(unsigned char), nbase, fh);
aecb11 165 //    free(compr);
SP 166     return defstream.total_out;
167 }
168
169 ts_uint ts_compress_string64(char *data, ts_uint data_len, char **encoded_compressed){
170     size_t nbase;
171     char *compr;
172     size_t number_of_compressed_bytes=ts_compress_data(data, data_len, &compr);
173     *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase);
f0bcea 174     free(compr);
SP 175     return nbase;
176 }
177
aecb11 178 char *ts_compress_intlist(int *data, ts_uint data_len){
SP 179     size_t nbase;
180     char *compr;
181     size_t number_of_compressed_bytes=ts_compress_data((char *)data, data_len*sizeof(int), &compr);
182     char *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase);
183     free(compr);
184     ts_uint header[4]={1, data_len, data_len, nbase};
185     char *encoded_header=(char *)base64_encode((unsigned char *)header, 4*sizeof(ts_uint), &nbase);
186     encoded_header=realloc(encoded_header, 4*sizeof(ts_uint)+strlen(encoded_compressed));
187     encoded_header=strcat(encoded_header,encoded_compressed);
188     free(encoded_compressed);
189     return encoded_header;
190 }
191
8c1bb1 192 ts_uint ts_decompress_string64(char *b64, ts_uint data_len, char **decompressed){
SP 193 return TS_SUCCESS;
194
195 }
f0bcea 196
SP 197 /* base64 encoding, taken from http://stackoverflow.com/questions/342409/how-do-i-base64-encode-decode-in-c */
198 static char encoding_table[] = {'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H',
199                                 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P',
200                                 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X',
201                                 'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f',
202                                 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n',
203                                 'o', 'p', 'q', 'r', 's', 't', 'u', 'v',
204                                 'w', 'x', 'y', 'z', '0', '1', '2', '3',
205                                 '4', '5', '6', '7', '8', '9', '+', '/'};
206 static char *decoding_table = NULL;
207 static int mod_table[] = {0, 2, 1};
208
209
7f6076 210 char *base64_encode(const unsigned char *data, size_t input_length, size_t *output_length) {
f0bcea 211     *output_length = 4 * ((input_length + 2) / 3);
SP 212     int i,j;
213     char *encoded_data = malloc(*output_length);
214     if (encoded_data == NULL) return NULL;
215     
216     for (i = 0, j = 0; i < input_length;) {
217
7f6076 218            uint32_t octet_a = i < input_length ? (unsigned char)data[i++] : 0;
SP 219         uint32_t octet_b = i < input_length ? (unsigned char)data[i++] : 0;
220            uint32_t octet_c = i < input_length ? (unsigned char)data[i++] : 0;
221         uint32_t triple = (octet_a << 0x10) + (octet_b << 0x08) + octet_c;
f0bcea 222
7f6076 223            encoded_data[j++] = encoding_table[(triple >> 3 * 6) & 0x3F];
SP 224            encoded_data[j++] = encoding_table[(triple >> 2 * 6) & 0x3F];
225            encoded_data[j++] = encoding_table[(triple >> 1 * 6) & 0x3F];
226            encoded_data[j++] = encoding_table[(triple >> 0 * 6) & 0x3F];
f0bcea 227     }
SP 228
229     for (i = 0; i < mod_table[input_length % 3]; i++)
7f6076 230            encoded_data[*output_length - 1 - i] = '=';
f0bcea 231
SP 232     return encoded_data;
233 }
234
235
7f6076 236 unsigned char *base64_decode(const char *data, size_t input_length, size_t *output_length) {
f0bcea 237     int i,j;
SP 238     if (decoding_table == NULL) build_decoding_table();
239
240     if (input_length % 4 != 0) return NULL;
241
242     *output_length = input_length / 4 * 3;
243     if (data[input_length - 1] == '=') (*output_length)--;
244     if (data[input_length - 2] == '=') (*output_length)--;
245
246     unsigned char *decoded_data = malloc(*output_length);
247     if (decoded_data == NULL) return NULL;
248
249     for (i = 0, j = 0; i < input_length;) {
7f6076 250            uint32_t sextet_a = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
SP 251            uint32_t sextet_b = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
252            uint32_t sextet_c = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
253            uint32_t sextet_d = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
f0bcea 254
7f6076 255            uint32_t triple = (sextet_a << 3 * 6)
f0bcea 256             + (sextet_b << 2 * 6)
SP 257             + (sextet_c << 1 * 6)
258             + (sextet_d << 0 * 6);
259
7f6076 260            if (j < *output_length) decoded_data[j++] = (triple >> 2 * 8) & 0xFF;
SP 261            if (j < *output_length) decoded_data[j++] = (triple >> 1 * 8) & 0xFF;
262            if (j < *output_length) decoded_data[j++] = (triple >> 0 * 8) & 0xFF;
f0bcea 263     }
8c1bb1 264     if(decoding_table !=NULL) free(decoding_table);
f0bcea 265     return decoded_data;
SP 266 }
267
268
269 void build_decoding_table() {
270
271     decoding_table = malloc(256);
272     int i;
273     for (i = 0; i < 64; i++)
7f6076 274         decoding_table[(unsigned char) encoding_table[i]] = i;
f0bcea 275 }
SP 276
277
278 void base64_cleanup() {
279     free(decoding_table);
280 }
281
282
283
284
285
286
287