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