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