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