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