Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2016-06-21 f3d3cb6eb7bc70c82aa0e57c693b259c26a7f84e
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
8c1bb1 2 #include <stdio.h>
SP 3 #include <string.h>
4 #include <stdlib.h>
5 #include <libxml/xmlmemory.h>
6 #include <libxml/parser.h>
7 #include <general.h>
8 #include <restore.h>
9 #include <snapshot.h>
10 #include <zlib.h>
11 #include "vesicle.h"
487968 12 #include "vertex.h"
a011d2 13 #include "triangle.h"
a3dbf0 14 #include "bond.h"
SP 15 #include "energy.h"
bb4033 16 #include "poly.h"
a3dbf0 17 #include "initial_distribution.h"
698ae1 18 #include "io.h"
a011d2 19
ab798b 20 ts_vesicle *parseDump(char *dumpfname) {
8c1bb1 21     xmlDocPtr doc;
8fa1c0 22     xmlNodePtr cur, cur1,cur2;
SP 23     ts_vesicle *vesicle=NULL;
8c1bb1 24     doc = xmlParseFile(dumpfname);
SP 25     
26     if (doc == NULL ) {
27         fatal("Dump file could not be found or parsed. It is correct file?",1);
28     }
29     
30     cur = xmlDocGetRootElement(doc);
31     
32     if (cur == NULL) {
33         fatal("Dump file is empty.",1);
34     }
35     
36     if (xmlStrcmp(cur->name, (const xmlChar *) "VTKFile")) {
37         fatal("document of the wrong type, root node != story",1);
38     }
39     
40     cur = cur->xmlChildrenNode;
41     while (cur != NULL) {
698ae1 42
SP 43         if ((!xmlStrcmp(cur->name, (const xmlChar *)"tape"))){
44             setGlobalTapeTXTfromTapeTag(doc, cur);
45         }
46
8c1bb1 47         if ((!xmlStrcmp(cur->name, (const xmlChar *)"trisurf"))){
3c772b 48             vesicle=parseTrisurfTag(doc, cur);
8c1bb1 49         }
8fa1c0 50         // START Point Position data &  Bonds
SP 51         if ((!xmlStrcmp(cur->name, (const xmlChar *)"UnstructuredGrid"))){
52             cur1 = cur->xmlChildrenNode;
53             while(cur1!=NULL){
54                 if ((!xmlStrcmp(cur1->name, (const xmlChar *)"Piece"))){
55                     cur2=cur1->xmlChildrenNode;
56                     while(cur2!=NULL){
57                         if ((!xmlStrcmp(cur2->name, (const xmlChar *)"Points"))){
720bd4 58                             //fprintf(stderr,"Found point data\n");
8fa1c0 59                             if(vesicle!=NULL)
f62178 60                                 //fprintf(stderr,"Fils: %ld, Nfono: %ld\n", vesicle->tape->nfil, vesicle->tape->nfono);
28efdb 61                                 parseXMLVertexPosition(vesicle, doc, cur2);
8fa1c0 62                         }
SP 63                         if ((!xmlStrcmp(cur2->name, (const xmlChar *)"Cells"))){
720bd4 64                         //fprintf(stderr,"Found cell(Bonds) data\n");
a3dbf0 65                             if(vesicle!=NULL)
28efdb 66                                 parseXMLBonds(vesicle, doc, cur2);
8fa1c0 67                         }
SP 68                         cur2=cur2->next;
69                     }    
70                 }
71                 cur1 = cur1->next;
72             }
73         }
74         // END Point Position data & Bonds
8c1bb1 75     cur = cur->next;
SP 76     }
77     
78     xmlFreeDoc(doc);
bb4033 79
86b69b 80 //    vesicle->poly_list=init_poly_list(0, 0, vesicle->vlist, vesicle);
bb4033 81
a3dbf0 82     init_normal_vectors(vesicle->tlist);
SP 83     mean_curvature_and_energy(vesicle);
84
ab798b 85 /* TODO: filaments */
a3dbf0 86
f4d6ca 87 //    ts_fprintf(stdout,"Restoration completed\n");
ab798b 88 //    write_vertex_xml_file(vesicle,999);
SP 89 //    vesicle_free(vesicle);
90 //    exit(0);
91     return vesicle;
8c1bb1 92 }
SP 93
698ae1 94 ts_bool setGlobalTapeTXTfromTapeTag(xmlDocPtr doc, xmlNodePtr cur){
SP 95     xmlChar *tape = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
96     strcpy(tapetxt,(char *)tape);
97     xmlFree(tape);
98     return TS_SUCCESS;
99 }
8fa1c0 100
SP 101
102 /* this is a parser of additional data in xml */
8c1bb1 103 ts_vesicle *parseTrisurfTag(xmlDocPtr doc, xmlNodePtr cur){
720bd4 104     //fprintf(stderr,"Parsing trisurf tag\n");
487968 105     xmlNodePtr child;
SP 106
107 #ifdef COMPRESS
8c1bb1 108     /* base64decode */
SP 109     size_t cLen;
110     /*size_t tLen;
111     const unsigned char test[]="Test";
112     char *cTest=base64_encode(test, 4,&tLen);
113     unsigned char *cuTest=base64_decode((char *)cTest,tLen,&tLen);
114     cuTest[tLen]=0;
115     fprintf(stderr,"%s\n",cuTest);
116     */
117     xmlChar *b64=xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
118     unsigned char *compressed=base64_decode((char *)b64,strlen((char *)b64)-1,&cLen);
119     /* uncompress */
120     unsigned char *subtree=(unsigned char *)malloc(512000*sizeof(unsigned char)); /* TODO: again, the uncompressed string must not exceed this */
121     z_stream infstream;
122     infstream.zalloc = Z_NULL;
123     infstream.zfree = Z_NULL;
124     infstream.opaque = Z_NULL;
125     infstream.avail_in = (ts_uint)cLen; // size of input
126         infstream.next_in = compressed; // input char array
127         infstream.avail_out = (ts_uint)512000; // size of output
128         infstream.next_out = subtree; // output char array
129      
130         // the actual DE-compression work.
131         inflateInit(&infstream);
132         inflate(&infstream, Z_NO_FLUSH);
133         inflateEnd(&infstream);    
720bd4 134     //fprintf(stderr,"%lu\n",cLen);
8c1bb1 135     subtree[infstream.total_out]='\0'; //zero terminate string    
720bd4 136     //fprintf(stderr,"%s\n",subtree);
8c1bb1 137     
SP 138     free(subtree);
487968 139 #endif
8c1bb1 140     /*parse xml subtree */
86b69b 141     xmlChar *nvtx, *npoly, *nmono;
8c1bb1 142     nvtx = xmlGetProp(cur, (xmlChar *)"nvtx");
SP 143     npoly=xmlGetProp(cur, (xmlChar *)"npoly");
86b69b 144     nmono=xmlGetProp(cur, (xmlChar *)"nmono");
698ae1 145     ts_tape *tape=parsetapebuffer(tapetxt);
SP 146     //fprintf(stderr,"nvtx=%u\n",atoi((char *)nvtx));
147     //TODO: check if nvtx is in agreement with nshell from tape
148     ts_vesicle *vesicle=init_vesicle(atoi((char *)nvtx),tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
8c1bb1 149     //vesicle->poly_list=init_poly_list(atoi((char *)npoly),atoi((char *)nmono), vesicle->vlist, vesicle);
SP 150     xmlFree(nvtx);
151     xmlFree(npoly);
86b69b 152     xmlFree(nmono);
487968 153
SP 154     child = cur->xmlChildrenNode;
155     while (child != NULL) {
156         if ((!xmlStrcmp(child->name, (const xmlChar *)"vtxn"))){
157             parseTrisurfVtxn(vesicle->vlist, doc, child);
158         }
a011d2 159         if ((!xmlStrcmp(child->name, (const xmlChar *)"tria"))){
SP 160             parseTrisurfTria(vesicle, doc, child);
161         }
e9ac7b 162         if ((!xmlStrcmp(child->name, (const xmlChar *)"trianeigh"))){
SP 163             parseTrisurfTriaNeigh(vesicle, doc, child);
164         }
a011d2 165          if ((!xmlStrcmp(child->name, (const xmlChar *)"tristar"))){
SP 166             parseTrisurfTristar(vesicle, doc, child);
167         }
4891eb 168          if ((!xmlStrcmp(child->name, (const xmlChar *)"nucleus"))){
SP 169             parseTrisurfNucleus(vesicle, doc, child);
170         }
49981c 171           if ((!xmlStrcmp(child->name, (const xmlChar *)"constant_volume"))){
SP 172             parseTrisurfConstantVolume(doc, child);
173         }
174           if ((!xmlStrcmp(child->name, (const xmlChar *)"constant_area"))){
175             parseTrisurfConstantArea(doc, child);
176         }
177
a011d2 178
487968 179     child = child->next;
SP 180     }
181
182
698ae1 183     vesicle->tape=tape;
SP 184     set_vesicle_values_from_tape(vesicle);
487968 185
8c1bb1 186     return vesicle;
SP 187 }
487968 188
a011d2 189
SP 190
191 /* Low level tags parsers */
49981c 192 ts_bool parseTrisurfConstantVolume(xmlDocPtr doc, xmlNodePtr cur){
SP 193     xmlChar *cvol = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
194     char *n=(char *)cvol;
195     V0=atof(n);
196     xmlFree(cvol);
197     return TS_SUCCESS;
198 }
199 ts_bool parseTrisurfConstantArea(xmlDocPtr doc, xmlNodePtr cur){
200     xmlChar *carea = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
201     char *n=(char *)carea;
202     A0=atof(n);
203     xmlFree(carea);
204     return TS_SUCCESS;
205 }
206
4891eb 207 ts_bool parseTrisurfNucleus(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur){
SP 208     xmlChar *coords = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
209     char *n=(char *)coords;
210     char *token=strtok(n," ");
211     ts_uint i;
212     for(i=0;i<3;i++){
213         vesicle->nucleus_center[i]=atof(token);
214         token=strtok(NULL," ");
215     }
216     xmlFree(coords);
217     return TS_SUCCESS;
218 }
a011d2 219
SP 220 ts_bool parseTrisurfVtxn(ts_vertex_list *vlist, xmlDocPtr doc, xmlNodePtr cur){
487968 221
SP 222     xmlChar *chari;
223     xmlChar *neighs;
224     char *n;
225     char *token;
226     ts_uint neighi;
227     ts_uint i;
228     chari = xmlGetProp(cur, (xmlChar *)"idx");
229     i=atoi((char *)chari);
230     xmlFree(chari);
231     ts_vertex *vtx=vlist->vtx[i];
232     neighs = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
233     //fprintf(stderr,"Found neigh for vtx %u that seems to have index %u with neighs=%s\n",i,vtx->idx,neighs);
234
235     n=(char *)neighs;
236     token=strtok(n," ");
237     while(token!=NULL){
238         neighi=atoi(token);
239         //fprintf(stderr,"%u", neighi);
240         vtx_add_neighbour(vtx,vlist->vtx[neighi]);
241         token=strtok(NULL," ");
242     }    
243     xmlFree(neighs);
244     return TS_SUCCESS;
245 }
246
a011d2 247 ts_bool parseTrisurfTria(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur){
SP 248     xmlChar *triangles;
249     char *tria;
250     char *vtx[3];
251     
e87e4e 252     ts_uint i,j;
a011d2 253     triangles = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
SP 254     tria=(char *)triangles;
e87e4e 255     vtx[0]=strtok(tria," ");
SP 256     for(i=1;i<3;i++)    vtx[i]=strtok(NULL," ");
257     j=0;
a011d2 258     while(vtx[2]!=NULL){
SP 259         triangle_add(vesicle->tlist, vesicle->vlist->vtx[atoi(vtx[0])],vesicle->vlist->vtx[atoi(vtx[1])],vesicle->vlist->vtx[atoi(vtx[2])]);
260         for(i=0;i<3;i++)    vtx[i]=strtok(NULL," ");
e87e4e 261         j++;
a011d2 262     }    
720bd4 263     //fprintf(stderr,"Parsing triangles %s j=%d\n",triangles,j);    
a011d2 264
SP 265     xmlFree(triangles);
266     return TS_SUCCESS;
267 }
268
269
e9ac7b 270 ts_bool parseTrisurfTriaNeigh(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur){
SP 271     xmlChar *triangles;
272     char *tria;
273     char *ntria[3];
274     ts_uint i,j;
275     triangles = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
276     tria=(char *)triangles;
e87e4e 277     ntria[0]=strtok(tria," ");
SP 278     for(i=1;i<3;i++)    ntria[i]=strtok(NULL," ");
e9ac7b 279     j=0;
SP 280     while(ntria[2]!=NULL){
281         triangle_add_neighbour(vesicle->tlist->tria[j],vesicle->tlist->tria[atoi(ntria[0])]);
282         triangle_add_neighbour(vesicle->tlist->tria[j],vesicle->tlist->tria[atoi(ntria[1])]);
283         triangle_add_neighbour(vesicle->tlist->tria[j],vesicle->tlist->tria[atoi(ntria[2])]);
284         j++;
285         for(i=0;i<3;i++)    ntria[i]=strtok(NULL," ");
286     }    
720bd4 287     //fprintf(stderr,"Parsing triangle neighbors j=%d\n",j);    
e9ac7b 288
SP 289     xmlFree(triangles);
290     return TS_SUCCESS;
291 }
292
293
a011d2 294 ts_bool parseTrisurfTristar(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur){
SP 295
296     xmlChar *chari;
297     xmlChar *tristar;
298     char *t;
299     char *token;
300     ts_uint neighi;
301     ts_uint i;
302     chari = xmlGetProp(cur, (xmlChar *)"idx");
303     i=atoi((char *)chari);
304     xmlFree(chari);
305     ts_vertex *vtx=vesicle->vlist->vtx[i];
306     tristar = xmlNodeListGetString(doc, cur->xmlChildrenNode, 1);
307 //    fprintf(stderr,"Found tristar for vtx %u that seems to have index %u with tristar=%s\n",i,vtx->idx,tristar);
308
309     t=(char *)tristar;
310     token=strtok(t," ");
311     while(token!=NULL){
312         neighi=atoi(token);
313         //fprintf(stderr,"%u", neighi);
314         vertex_add_tristar(vtx,vesicle->tlist->tria[neighi]);
315         token=strtok(NULL," ");
316     }    
317     xmlFree(tristar);
318     return TS_SUCCESS;
319 }
320
8fa1c0 321
SP 322 /* this is a parser of vertex positions and bonds from main xml data */
323 ts_bool parseXMLVertexPosition(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
324     xmlNodePtr child = cur->xmlChildrenNode;
325     xmlChar *points;
326     char *pts;
301c87 327     int i, idx, polyidx, monoidx, filidx, fonoidx;
8fa1c0 328     char *token[3];
SP 329     while (child != NULL) {
330         if ((!xmlStrcmp(child->name, (const xmlChar *)"DataArray"))){
331             points = xmlNodeListGetString(doc, child->xmlChildrenNode, 1);
332             pts=(char *)points;
e87e4e 333             token[0]=strtok(pts," ");
ecdd71 334             token[1]=strtok(NULL," ");
SP 335             token[2]=strtok(NULL,"\n");
8fa1c0 336             idx=0;
SP 337             while(token[0]!=NULL){
86b69b 338                 if(idx<vesicle->vlist->n){
SP 339                     vesicle->vlist->vtx[idx]->x=atof(token[0]);
340                     vesicle->vlist->vtx[idx]->y=atof(token[1]);
341                     vesicle->vlist->vtx[idx]->z=atof(token[2]);
301c87 342                 } else if(vesicle->tape->nmono && vesicle->tape->npoly && idx<vesicle->vlist->n+vesicle->tape->nmono*vesicle->tape->npoly) {
86b69b 343                     polyidx=(idx-vesicle->vlist->n)/vesicle->tape->nmono;
SP 344                     monoidx=(idx-vesicle->vlist->n)%vesicle->tape->nmono;
345                     vesicle->poly_list->poly[polyidx]->vlist->vtx[monoidx]->x=atof(token[0]);
346                     vesicle->poly_list->poly[polyidx]->vlist->vtx[monoidx]->y=atof(token[1]);
347                     vesicle->poly_list->poly[polyidx]->vlist->vtx[monoidx]->z=atof(token[2]);
301c87 348                 } else {
SP 349                     filidx=(idx-vesicle->vlist->n-vesicle->tape->nmono*vesicle->tape->npoly)/vesicle->tape->nfono;
350                     fonoidx=(idx-vesicle->vlist->n-vesicle->tape->nmono*vesicle->tape->npoly)%vesicle->tape->nfono;
f62178 351                     //fprintf(stderr,"filidx=%d, fonoidx=%d, coord=%s,%s,%s\n",filidx,fonoidx,token[0],token[1],token[2]);
301c87 352                     vesicle->filament_list->poly[filidx]->vlist->vtx[fonoidx]->x=atof(token[0]);
SP 353                     vesicle->filament_list->poly[filidx]->vlist->vtx[fonoidx]->y=atof(token[1]);
354                     vesicle->filament_list->poly[filidx]->vlist->vtx[fonoidx]->z=atof(token[2]);
86b69b 355                 }
ecdd71 356                 for(i=0;i<2;i++)    token[i]=strtok(NULL," ");    
SP 357                 token[2]=strtok(NULL,"\n");
8fa1c0 358                 idx++;
SP 359             }
360             xmlFree(points);
361         }
362         child=child->next;
363     }
f62178 364     //fprintf(stderr,"Came here\n");
720bd4 365     //fprintf(stderr,"Vertices position j=%d\n",idx);    
28efdb 366
8fa1c0 367     return TS_SUCCESS;
SP 368 }
a3dbf0 369 ts_bool parseXMLBonds(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
SP 370     xmlNodePtr child = cur->xmlChildrenNode;
bb4033 371     xmlChar *bonds, *conname;
a3dbf0 372     char *b;
86b69b 373     int idx, polyidx;
a3dbf0 374     char *token[2];
SP 375     while (child != NULL) {
bb4033 376         conname=xmlGetProp(child, (xmlChar *)"Name");
SP 377         if ((!xmlStrcmp(child->name, (const xmlChar *)"DataArray")) && !xmlStrcmp(conname, (const xmlChar *)"connectivity") ){
a3dbf0 378             bonds = xmlNodeListGetString(doc, child->xmlChildrenNode, 1);
SP 379             b=(char *)bonds;
e87e4e 380             token[0]=strtok(b," ");
ecdd71 381             token[1]=strtok(NULL,"\n");
a3dbf0 382             idx=0;
SP 383             while(token[0]!=NULL){
86b69b 384                 if(idx<3*(vesicle->vlist->n-2)){
SP 385                     bond_add(vesicle->blist, vesicle->vlist->vtx[atoi(token[0])], vesicle->vlist->vtx[atoi(token[1])]);
386                 }
387                 else {
388                 //find grafted vtx
301c87 389                     if(vesicle->tape->npoly && vesicle->tape->nmono && (vesicle->tape->nmono-1)==(idx-3*(vesicle->vlist->n-2))%(vesicle->tape->nmono)
SP 390                         && idx<(3*vesicle->vlist->n-2+vesicle->tape->nmono*vesicle->tape->npoly+vesicle->tape->npoly)){
86b69b 391                         polyidx=(idx-3*(vesicle->vlist->n-2))/(vesicle->tape->nmono);
720bd4 392                         //fprintf(stderr,"poly=%d, vertex=%d\n",polyidx,atoi(token[0]));
86b69b 393                         vesicle->poly_list->poly[polyidx]->grafted_vtx=vesicle->vlist->vtx[atoi(token[0])];
SP 394                         vesicle->vlist->vtx[atoi(token[0])]->grafted_poly=vesicle->poly_list->poly[polyidx];
395                     }
396                 }
ecdd71 397                 token[0]=strtok(NULL," ");    
SP 398                 token[1]=strtok(NULL,"\n");    
a3dbf0 399                 idx++;
SP 400             }
401             xmlFree(bonds);
402         }
bb4033 403         xmlFree(conname);
a3dbf0 404         child=child->next;
SP 405     }
720bd4 406     //fprintf(stderr,"Bond data j=%d\n",idx);    
a3dbf0 407     return TS_SUCCESS;
SP 408 }
409
8fa1c0 410
SP 411