Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2016-02-29 ecdd71929b38a504e2b0ae1a9c105a5e1ac28e4a
commit | author | age
8c1bb1 1 #include <stdio.h>
SP 2 #include <string.h>
3 #include <stdlib.h>
4 #include <libxml/xmlmemory.h>
5 #include <libxml/parser.h>
6 #include <general.h>
7 #include <restore.h>
8 #include <snapshot.h>
9 #include <zlib.h>
10 #include "vesicle.h"
487968 11 #include "vertex.h"
a011d2 12 #include "triangle.h"
a3dbf0 13 #include "bond.h"
SP 14 #include "energy.h"
15 #include "initial_distribution.h"
a011d2 16
8c1bb1 17 ts_bool parseDump(char *dumpfname) {
SP 18     xmlDocPtr doc;
8fa1c0 19     xmlNodePtr cur, cur1,cur2;
SP 20     ts_vesicle *vesicle=NULL;
8c1bb1 21
SP 22     doc = xmlParseFile(dumpfname);
23     
24     if (doc == NULL ) {
25         fatal("Dump file could not be found or parsed. It is correct file?",1);
26     }
27     
28     cur = xmlDocGetRootElement(doc);
29     
30     if (cur == NULL) {
31         fatal("Dump file is empty.",1);
32     }
33     
34     if (xmlStrcmp(cur->name, (const xmlChar *) "VTKFile")) {
35         fatal("document of the wrong type, root node != story",1);
36     }
37     
38     cur = cur->xmlChildrenNode;
39     while (cur != NULL) {
40         if ((!xmlStrcmp(cur->name, (const xmlChar *)"trisurf"))){
3c772b 41             vesicle=parseTrisurfTag(doc, cur);
8c1bb1 42         }
8fa1c0 43         // START Point Position data &  Bonds
SP 44         if ((!xmlStrcmp(cur->name, (const xmlChar *)"UnstructuredGrid"))){
45             cur1 = cur->xmlChildrenNode;
46             while(cur1!=NULL){
47                 if ((!xmlStrcmp(cur1->name, (const xmlChar *)"Piece"))){
48                     cur2=cur1->xmlChildrenNode;
49                     while(cur2!=NULL){
50                         if ((!xmlStrcmp(cur2->name, (const xmlChar *)"Points"))){
51                             fprintf(stderr,"Found point data\n");
52                             if(vesicle!=NULL)
28efdb 53                                 parseXMLVertexPosition(vesicle, doc, cur2);
8fa1c0 54                         }
SP 55                         if ((!xmlStrcmp(cur2->name, (const xmlChar *)"Cells"))){
56                         fprintf(stderr,"Found cell(Bonds) data\n");
a3dbf0 57                             if(vesicle!=NULL)
28efdb 58                                 parseXMLBonds(vesicle, doc, cur2);
8fa1c0 59                         }
SP 60                         cur2=cur2->next;
61                     }    
62                 }
63                 cur1 = cur1->next;
64             }
65         }
66         // END Point Position data & Bonds
8c1bb1 67     cur = cur->next;
SP 68     }
69     
70     xmlFreeDoc(doc);
a3dbf0 71     init_normal_vectors(vesicle->tlist);
SP 72     mean_curvature_and_energy(vesicle);
73
74 /* TODO: cells, polymeres, filaments, core, tape */
75
8c1bb1 76     fprintf(stderr,"Restoration completed\n");
28efdb 77     write_vertex_xml_file(vesicle,999);
3c772b 78     vesicle_free(vesicle);
e87e4e 79     exit(0);
8c1bb1 80     return TS_SUCCESS;
SP 81 }
82
8fa1c0 83
SP 84
85 /* this is a parser of additional data in xml */
8c1bb1 86 ts_vesicle *parseTrisurfTag(xmlDocPtr doc, xmlNodePtr cur){
SP 87     fprintf(stderr,"Parsing trisurf tag\n");
487968 88     xmlNodePtr child;
SP 89
90 #ifdef COMPRESS
8c1bb1 91     /* base64decode */
SP 92     size_t cLen;
93     /*size_t tLen;
94     const unsigned char test[]="Test";
95     char *cTest=base64_encode(test, 4,&tLen);
96     unsigned char *cuTest=base64_decode((char *)cTest,tLen,&tLen);
97     cuTest[tLen]=0;
98     fprintf(stderr,"%s\n",cuTest);
99     */
100     xmlChar *b64=xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
101     unsigned char *compressed=base64_decode((char *)b64,strlen((char *)b64)-1,&cLen);
102     /* uncompress */
103     unsigned char *subtree=(unsigned char *)malloc(512000*sizeof(unsigned char)); /* TODO: again, the uncompressed string must not exceed this */
104     z_stream infstream;
105     infstream.zalloc = Z_NULL;
106     infstream.zfree = Z_NULL;
107     infstream.opaque = Z_NULL;
108     infstream.avail_in = (ts_uint)cLen; // size of input
109         infstream.next_in = compressed; // input char array
110         infstream.avail_out = (ts_uint)512000; // size of output
111         infstream.next_out = subtree; // output char array
112      
113         // the actual DE-compression work.
114         inflateInit(&infstream);
115         inflate(&infstream, Z_NO_FLUSH);
116         inflateEnd(&infstream);    
117     fprintf(stderr,"%lu\n",cLen);
118     subtree[infstream.total_out]='\0'; //zero terminate string    
119     fprintf(stderr,"%s\n",subtree);
120     
121     free(subtree);
487968 122 #endif
8c1bb1 123     /*parse xml subtree */
SP 124     xmlChar *nvtx, *npoly, *nfono;
125     nvtx = xmlGetProp(cur, (xmlChar *)"nvtx");
126     npoly=xmlGetProp(cur, (xmlChar *)"npoly");
127     nfono=xmlGetProp(cur, (xmlChar *)"nfono");
128     fprintf(stderr,"nvtx=%u\n",atoi((char *)nvtx));
129     ts_vesicle *vesicle=init_vesicle(atoi((char *)nvtx),10,10,10,0.1);
130     //vesicle->poly_list=init_poly_list(atoi((char *)npoly),atoi((char *)nmono), vesicle->vlist, vesicle);
131     xmlFree(nvtx);
132     xmlFree(npoly);
133     xmlFree(nfono);
487968 134
SP 135     child = cur->xmlChildrenNode;
136     while (child != NULL) {
137         if ((!xmlStrcmp(child->name, (const xmlChar *)"vtxn"))){
138             parseTrisurfVtxn(vesicle->vlist, doc, child);
139         }
a011d2 140         if ((!xmlStrcmp(child->name, (const xmlChar *)"tria"))){
SP 141             parseTrisurfTria(vesicle, doc, child);
142         }
e9ac7b 143         if ((!xmlStrcmp(child->name, (const xmlChar *)"trianeigh"))){
SP 144             parseTrisurfTriaNeigh(vesicle, doc, child);
145         }
a011d2 146          if ((!xmlStrcmp(child->name, (const xmlChar *)"tristar"))){
SP 147             parseTrisurfTristar(vesicle, doc, child);
148         }
149
487968 150     child = child->next;
SP 151     }
152
153
154
8c1bb1 155     return vesicle;
SP 156 }
487968 157
a011d2 158
SP 159
160 /* Low level tags parsers */
161
162 ts_bool parseTrisurfVtxn(ts_vertex_list *vlist, xmlDocPtr doc, xmlNodePtr cur){
487968 163
SP 164     xmlChar *chari;
165     xmlChar *neighs;
166     char *n;
167     char *token;
168     ts_uint neighi;
169     ts_uint i;
170     chari = xmlGetProp(cur, (xmlChar *)"idx");
171     i=atoi((char *)chari);
172     xmlFree(chari);
173     ts_vertex *vtx=vlist->vtx[i];
174     neighs = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
175     //fprintf(stderr,"Found neigh for vtx %u that seems to have index %u with neighs=%s\n",i,vtx->idx,neighs);
176
177     n=(char *)neighs;
178     token=strtok(n," ");
179     while(token!=NULL){
180         neighi=atoi(token);
181         //fprintf(stderr,"%u", neighi);
182         vtx_add_neighbour(vtx,vlist->vtx[neighi]);
183         token=strtok(NULL," ");
184     }    
185     xmlFree(neighs);
186     return TS_SUCCESS;
187 }
188
a011d2 189 ts_bool parseTrisurfTria(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur){
SP 190     xmlChar *triangles;
191     char *tria;
192     char *vtx[3];
193     
e87e4e 194     ts_uint i,j;
a011d2 195     triangles = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
SP 196     tria=(char *)triangles;
e87e4e 197     vtx[0]=strtok(tria," ");
SP 198     for(i=1;i<3;i++)    vtx[i]=strtok(NULL," ");
199     j=0;
a011d2 200     while(vtx[2]!=NULL){
SP 201         triangle_add(vesicle->tlist, vesicle->vlist->vtx[atoi(vtx[0])],vesicle->vlist->vtx[atoi(vtx[1])],vesicle->vlist->vtx[atoi(vtx[2])]);
202         for(i=0;i<3;i++)    vtx[i]=strtok(NULL," ");
e87e4e 203         j++;
a011d2 204     }    
e87e4e 205     fprintf(stderr,"Parsing triangles %s j=%d\n",triangles,j);    
a011d2 206
SP 207     xmlFree(triangles);
208     return TS_SUCCESS;
209 }
210
211
e9ac7b 212 ts_bool parseTrisurfTriaNeigh(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur){
SP 213     xmlChar *triangles;
214     char *tria;
215     char *ntria[3];
216     ts_uint i,j;
217     triangles = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
218     tria=(char *)triangles;
e87e4e 219     ntria[0]=strtok(tria," ");
SP 220     for(i=1;i<3;i++)    ntria[i]=strtok(NULL," ");
e9ac7b 221     j=0;
SP 222     while(ntria[2]!=NULL){
223         triangle_add_neighbour(vesicle->tlist->tria[j],vesicle->tlist->tria[atoi(ntria[0])]);
224         triangle_add_neighbour(vesicle->tlist->tria[j],vesicle->tlist->tria[atoi(ntria[1])]);
225         triangle_add_neighbour(vesicle->tlist->tria[j],vesicle->tlist->tria[atoi(ntria[2])]);
226         j++;
227         for(i=0;i<3;i++)    ntria[i]=strtok(NULL," ");
228     }    
229     fprintf(stderr,"Parsing triangle neighbors j=%d\n",j);    
230
231     xmlFree(triangles);
232     return TS_SUCCESS;
233 }
234
235
a011d2 236 ts_bool parseTrisurfTristar(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur){
SP 237
238     xmlChar *chari;
239     xmlChar *tristar;
240     char *t;
241     char *token;
242     ts_uint neighi;
243     ts_uint i;
244     chari = xmlGetProp(cur, (xmlChar *)"idx");
245     i=atoi((char *)chari);
246     xmlFree(chari);
247     ts_vertex *vtx=vesicle->vlist->vtx[i];
248     tristar = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
249 //    fprintf(stderr,"Found tristar for vtx %u that seems to have index %u with tristar=%s\n",i,vtx->idx,tristar);
250
251     t=(char *)tristar;
252     token=strtok(t," ");
253     while(token!=NULL){
254         neighi=atoi(token);
255         //fprintf(stderr,"%u", neighi);
256         vertex_add_tristar(vtx,vesicle->tlist->tria[neighi]);
257         token=strtok(NULL," ");
258     }    
259     xmlFree(tristar);
260     return TS_SUCCESS;
261 }
262
8fa1c0 263
SP 264 /* this is a parser of vertex positions and bonds from main xml data */
265 ts_bool parseXMLVertexPosition(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
266     xmlNodePtr child = cur->xmlChildrenNode;
267     xmlChar *points;
268     char *pts;
269     int i, idx;
270     char *token[3];
271     while (child != NULL) {
272         if ((!xmlStrcmp(child->name, (const xmlChar *)"DataArray"))){
273             points = xmlNodeListGetString(doc, child->xmlChildrenNode, 1);
274             pts=(char *)points;
e87e4e 275             token[0]=strtok(pts," ");
ecdd71 276             token[1]=strtok(NULL," ");
SP 277             token[2]=strtok(NULL,"\n");
8fa1c0 278             idx=0;
SP 279             while(token[0]!=NULL){
280                 vesicle->vlist->vtx[idx]->x=atof(token[0]);
281                 vesicle->vlist->vtx[idx]->y=atof(token[1]);
282                 vesicle->vlist->vtx[idx]->z=atof(token[2]);
ecdd71 283                 for(i=0;i<2;i++)    token[i]=strtok(NULL," ");    
SP 284                 token[2]=strtok(NULL,"\n");
8fa1c0 285                 idx++;
SP 286             }
287             xmlFree(points);
288         }
289         child=child->next;
290     }
28efdb 291     fprintf(stderr,"Vertices position j=%d\n",idx);    
SP 292
8fa1c0 293     return TS_SUCCESS;
SP 294 }
a3dbf0 295 ts_bool parseXMLBonds(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
SP 296     xmlNodePtr child = cur->xmlChildrenNode;
297     xmlChar *bonds;
298     char *b;
ecdd71 299     int idx;
a3dbf0 300     char *token[2];
SP 301     while (child != NULL) {
28efdb 302         if ((!xmlStrcmp(child->name, (const xmlChar *)"DataArray")) && !xmlStrcmp(xmlGetProp(child, (xmlChar *)"Name"), (const xmlChar *)"connectivity") ){
a3dbf0 303             bonds = xmlNodeListGetString(doc, child->xmlChildrenNode, 1);
SP 304             b=(char *)bonds;
e87e4e 305             token[0]=strtok(b," ");
ecdd71 306             token[1]=strtok(NULL,"\n");
a3dbf0 307             idx=0;
SP 308             while(token[0]!=NULL){
309                 bond_add(vesicle->blist, vesicle->vlist->vtx[atoi(token[0])], vesicle->vlist->vtx[atoi(token[1])]);
ecdd71 310                 token[0]=strtok(NULL," ");    
SP 311                 token[1]=strtok(NULL,"\n");    
a3dbf0 312                 idx++;
SP 313             }
314             xmlFree(bonds);
315         }
316         child=child->next;
317     }
28efdb 318     fprintf(stderr,"Bond data j=%d\n",idx);    
a3dbf0 319     return TS_SUCCESS;
SP 320 }
321
8fa1c0 322
SP 323