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