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