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