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