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