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