Trisurf Monte Carlo simulator
Samo Penic
2016-06-21 032273a4cf21e1a43038b60b8744021930c97a5b
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
ee84bd 2 #include <string.h>
SP 3 #include "general.h"
4 #include <stdio.h>
5 #include "initial_distribution.h"
6 #include "vesicle.h"
7 #include "dumpstate.h"
8 #include <libxml/parser.h>
9 #include <libxml/tree.h>
43e534 10 #include "vertex.h"
SP 11 #include "energy.h"
ee84bd 12
SP 13 ts_vesicle *vtk2vesicle(char *filename, ts_tape *tape){
14
15     ts_uint nshell=tape->nshell;
16     ts_uint ncmax1=tape->ncxmax;
17     ts_uint ncmax2=tape->ncymax;
18     ts_uint ncmax3=tape->nczmax;
19     ts_double stepsize=tape->stepsize;
20
21     ts_uint no_vertices=5*nshell*nshell+2;
22     ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize);
23     vesicle->nshell=nshell;
24     parse_vtk(filename, vesicle);
25     exit(1);
26     return vesicle;
27 }
28
29
30 ts_bool parse_vtk(char *filename, ts_vesicle *vesicle){
31     xmlDoc *doc;
32     xmlNode *root_element=NULL;
33     xmlNode *cur_node = NULL;
34     doc = xmlReadFile(filename, NULL, 0);
35     root_element=xmlDocGetRootElement(doc);
36     cur_node=root_element->children;
37     while(cur_node!=NULL){
38 //        fprintf(stderr,"Node name is: %s\n",cur_node->name);
39         if(strcmp((char *)cur_node->name,"UnstructuredGrid")==0) break;
40         cur_node=cur_node->next;
41     }
42
43     cur_node=cur_node->children;
44     while(cur_node!=NULL){
45 //        fprintf(stderr,"Node name is: %s\n",cur_node->name);
46         cur_node=cur_node->next;
47         if(strcmp((char *)cur_node->name,"Piece")==0) break;
48     }
49
50     cur_node=cur_node->children;
51     while(cur_node!=NULL){
43e534 52         //fprintf(stderr,"Node name is: %s\n",cur_node->name);
ee84bd 53         cur_node=cur_node->next;
SP 54         if(strcmp((char *)cur_node->name,"PointData")==0) vtk_index2vesicle(cur_node->children->next->children, vesicle);
43e534 55         if(strcmp((char *)cur_node->name,"Points")==0) vtk_coordinates(cur_node->children->next->children,vesicle);
SP 56         if(strcmp((char *)cur_node->name,"Cells")==0) {
57         vtk_neighbours(cur_node->children->next->children,vesicle);
58             break; //segfaults, because it finds another cells. Why?
59     }    
ee84bd 60     }
43e534 61     //we have vertices and neighbour relations all set, but unordered. Let's do bonds, ordering, triangles...
SP 62     vesicle->vlist = vtk_sort_neighbours(vesicle->blist,vesicle->vlist);
63     init_triangles(vesicle);
64     init_triangle_neighbours(vesicle);
65     init_common_vertex_triangle_neighbours(vesicle);
66     init_normal_vectors(vesicle->tlist);
67     mean_curvature_and_energy(vesicle);
68     ts_fprintf(stdout,"restoring from vtk dump finished!\n");
ee84bd 69     return TS_SUCCESS;
SP 70 }
71
72
73 ts_bool vtk_index2vesicle(xmlNode *node, ts_vesicle *vesicle){
74     //fprintf(stderr, "vsebina: %s\n",node->content);
75     ts_uint i;
76     char *token;
77     token = strtok((char *)node->content, " ");
78     for(i=0;i<vesicle->vlist->n;i++){
79         vesicle->vlist->vtx[i]->idx=atoi(token);
80         token=strtok(NULL," ");
81     }
82     //fprintf(stderr,"idx[11]=%d\n",vesicle->vlist->vtx[11]->idx);
83     return TS_SUCCESS;
84 }
43e534 85
SP 86 ts_bool vtk_coordinates(xmlNode *node, ts_vesicle *vesicle){
87     //fprintf(stderr, "vsebina: %s\n",node->content);
88     ts_uint i;
89     char *token;
90     token = strtok((char *)node->content, " ");
91     for(i=0;i<vesicle->vlist->n;i++){
92         vesicle->vlist->vtx[i]->x=atof(token);
93 //    fprintf(stderr,"%e, ",atof(token));
94         token=strtok(NULL," ");
95         vesicle->vlist->vtx[i]->y=atof(token);
96         token=strtok(NULL,"\n");
97         vesicle->vlist->vtx[i]->z=atof(token);
98         if(i>0) token=strtok(NULL," ");
99     }
100     return TS_SUCCESS;
101
102 }
103
104 ts_bool vtk_neighbours(xmlNode *node, ts_vesicle *vesicle){
105  //    fprintf(stderr, "vsebina: %s\n",node->content);
106     ts_uint i;
107     ts_uint idx1, idx2;
108     char *token;
109     token = strtok((char *)node->content, " ");
110     for(i=0;i<3*(vesicle->vlist->n-2);i++){
111     idx1=atoi(token);
112         token=strtok(NULL,"\n");
113     idx2=atoi(token);
114         vtx_add_neighbour(vesicle->vlist->vtx[idx1], vesicle->vlist->vtx[idx2]);
115         vtx_add_neighbour(vesicle->vlist->vtx[idx2], vesicle->vlist->vtx[idx1]);
116     fprintf(stderr, "%d. povezujem %d in %d\n",i,idx1,idx2);
117         if(i>0) token=strtok(NULL," ");
118     }
119     return TS_SUCCESS;
120 }
121
122 /* this is almost exact copy of init_sort_neighbours */
123 ts_vertex_list *vtk_sort_neighbours(ts_bond_list *blist,ts_vertex_list *vlist){
124     ts_vertex **vtx=vlist->vtx -1; // take a look at dipyramid function for comment.
125     ts_uint i,l,j,jj,jjj,k=0;   
126     ts_double eps=0.001; // Take a look if EPS from math.h can be used
127
128 /*lets initialize memory for temporary vertex_list. Should we write a function instead */
129     ts_vertex_list *tvlist=vertex_list_copy(vlist);
130     ts_vertex **tvtx=tvlist->vtx -1;  /* again to compensate for 0-indexing */
131
132     ts_double dist2; // Square of distance of neighbours
133     ts_double direct; // Something, dont know what, but could be normal of some kind
134     for(i=1;i<=vlist->n;i++){
135         k++; // WHY i IS NOT GOOD??
136            vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->neigh[0]->idx+1]); //always add 1st
137            jjj=1;
138            jj=1;
139            for(l=2;l<=vtx[i]->neigh_no;l++){
140                for(j=2;j<=vtx[i]->neigh_no;j++){
141                    dist2=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
142                    direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
143 // TODO: check if fabs can be used with all floating point types!!
144                    if( (direct>0.0) && (j!=jjj) ){
145                        vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->neigh[j-1]->idx+1]);
146                        jjj=jj;
147                        jj=j;
148                    //    break;
149                    }
150                }
151            }
152     if(vtx[i]->neigh_no!=tvtx[i]->neigh_no){
153         fprintf(stderr,"Hej! pri vtx=%d se razikuje stevilo sosedov (prej=%d, potem=%d)\n",i, vtx[i]->neigh_no, tvtx[i]->neigh_no);
154         }
155     }
156     if(eps>dist2);
157 /* We use the temporary vertex for our main vertices and we abandon main
158  * vertices, because their neighbours are not correctly ordered */
159    // tvtx=vlist->vtx;
160    // vlist->vtx=tvtx;
161    // tvlist->vtx=vtx;
162     vtx_list_free(vlist);
163 /* Let's make a check if the number of bonds is correct */
164     if((blist->n)!=3*(tvlist->n-2)){
165         ts_fprintf(stderr,"Number of bonds is %u should be %u!\n", blist->n, 3*(tvlist->n-2));
166         fatal("Number of bonds is not 3*(no_vertex-2).",4);
167     }
168
169     return tvlist;
170 }
171
172