#include <string.h>
|
#include "general.h"
|
#include <stdio.h>
|
#include "initial_distribution.h"
|
#include "vesicle.h"
|
#include "dumpstate.h"
|
#include <libxml/parser.h>
|
#include <libxml/tree.h>
|
#include "vertex.h"
|
#include "energy.h"
|
|
ts_vesicle *vtk2vesicle(char *filename, ts_tape *tape){
|
|
ts_uint nshell=tape->nshell;
|
ts_uint ncmax1=tape->ncxmax;
|
ts_uint ncmax2=tape->ncymax;
|
ts_uint ncmax3=tape->nczmax;
|
ts_double stepsize=tape->stepsize;
|
|
ts_uint no_vertices=5*nshell*nshell+2;
|
ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize);
|
vesicle->nshell=nshell;
|
parse_vtk(filename, vesicle);
|
exit(1);
|
return vesicle;
|
}
|
|
|
ts_bool parse_vtk(char *filename, ts_vesicle *vesicle){
|
xmlDoc *doc;
|
xmlNode *root_element=NULL;
|
xmlNode *cur_node = NULL;
|
doc = xmlReadFile(filename, NULL, 0);
|
root_element=xmlDocGetRootElement(doc);
|
cur_node=root_element->children;
|
while(cur_node!=NULL){
|
// fprintf(stderr,"Node name is: %s\n",cur_node->name);
|
if(strcmp((char *)cur_node->name,"UnstructuredGrid")==0) break;
|
cur_node=cur_node->next;
|
}
|
|
cur_node=cur_node->children;
|
while(cur_node!=NULL){
|
// fprintf(stderr,"Node name is: %s\n",cur_node->name);
|
cur_node=cur_node->next;
|
if(strcmp((char *)cur_node->name,"Piece")==0) break;
|
}
|
|
cur_node=cur_node->children;
|
while(cur_node!=NULL){
|
//fprintf(stderr,"Node name is: %s\n",cur_node->name);
|
cur_node=cur_node->next;
|
if(strcmp((char *)cur_node->name,"PointData")==0) vtk_index2vesicle(cur_node->children->next->children, vesicle);
|
if(strcmp((char *)cur_node->name,"Points")==0) vtk_coordinates(cur_node->children->next->children,vesicle);
|
if(strcmp((char *)cur_node->name,"Cells")==0) {
|
vtk_neighbours(cur_node->children->next->children,vesicle);
|
break; //segfaults, because it finds another cells. Why?
|
}
|
}
|
//we have vertices and neighbour relations all set, but unordered. Let's do bonds, ordering, triangles...
|
vesicle->vlist = vtk_sort_neighbours(vesicle->blist,vesicle->vlist);
|
init_triangles(vesicle);
|
init_triangle_neighbours(vesicle);
|
init_common_vertex_triangle_neighbours(vesicle);
|
init_normal_vectors(vesicle->tlist);
|
mean_curvature_and_energy(vesicle);
|
ts_fprintf(stdout,"restoring from vtk dump finished!\n");
|
return TS_SUCCESS;
|
}
|
|
|
ts_bool vtk_index2vesicle(xmlNode *node, ts_vesicle *vesicle){
|
//fprintf(stderr, "vsebina: %s\n",node->content);
|
ts_uint i;
|
char *token;
|
token = strtok((char *)node->content, " ");
|
for(i=0;i<vesicle->vlist->n;i++){
|
vesicle->vlist->vtx[i]->idx=atoi(token);
|
token=strtok(NULL," ");
|
}
|
//fprintf(stderr,"idx[11]=%d\n",vesicle->vlist->vtx[11]->idx);
|
return TS_SUCCESS;
|
}
|
|
ts_bool vtk_coordinates(xmlNode *node, ts_vesicle *vesicle){
|
//fprintf(stderr, "vsebina: %s\n",node->content);
|
ts_uint i;
|
char *token;
|
token = strtok((char *)node->content, " ");
|
for(i=0;i<vesicle->vlist->n;i++){
|
vesicle->vlist->vtx[i]->x=atof(token);
|
// fprintf(stderr,"%e, ",atof(token));
|
token=strtok(NULL," ");
|
vesicle->vlist->vtx[i]->y=atof(token);
|
token=strtok(NULL,"\n");
|
vesicle->vlist->vtx[i]->z=atof(token);
|
if(i>0) token=strtok(NULL," ");
|
}
|
return TS_SUCCESS;
|
|
}
|
|
ts_bool vtk_neighbours(xmlNode *node, ts_vesicle *vesicle){
|
// fprintf(stderr, "vsebina: %s\n",node->content);
|
ts_uint i;
|
ts_uint idx1, idx2;
|
char *token;
|
token = strtok((char *)node->content, " ");
|
for(i=0;i<3*(vesicle->vlist->n-2);i++){
|
idx1=atoi(token);
|
token=strtok(NULL,"\n");
|
idx2=atoi(token);
|
vtx_add_neighbour(vesicle->vlist->vtx[idx1], vesicle->vlist->vtx[idx2]);
|
vtx_add_neighbour(vesicle->vlist->vtx[idx2], vesicle->vlist->vtx[idx1]);
|
fprintf(stderr, "%d. povezujem %d in %d\n",i,idx1,idx2);
|
if(i>0) token=strtok(NULL," ");
|
}
|
return TS_SUCCESS;
|
}
|
|
/* this is almost exact copy of init_sort_neighbours */
|
ts_vertex_list *vtk_sort_neighbours(ts_bond_list *blist,ts_vertex_list *vlist){
|
ts_vertex **vtx=vlist->vtx -1; // take a look at dipyramid function for comment.
|
ts_uint i,l,j,jj,jjj,k=0;
|
ts_double eps=0.001; // Take a look if EPS from math.h can be used
|
|
/*lets initialize memory for temporary vertex_list. Should we write a function instead */
|
ts_vertex_list *tvlist=vertex_list_copy(vlist);
|
ts_vertex **tvtx=tvlist->vtx -1; /* again to compensate for 0-indexing */
|
|
ts_double dist2; // Square of distance of neighbours
|
ts_double direct; // Something, dont know what, but could be normal of some kind
|
for(i=1;i<=vlist->n;i++){
|
k++; // WHY i IS NOT GOOD??
|
vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->neigh[0]->idx+1]); //always add 1st
|
jjj=1;
|
jj=1;
|
for(l=2;l<=vtx[i]->neigh_no;l++){
|
for(j=2;j<=vtx[i]->neigh_no;j++){
|
dist2=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
|
direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]);
|
// TODO: check if fabs can be used with all floating point types!!
|
if( (direct>0.0) && (j!=jjj) ){
|
vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->neigh[j-1]->idx+1]);
|
jjj=jj;
|
jj=j;
|
// break;
|
}
|
}
|
}
|
if(vtx[i]->neigh_no!=tvtx[i]->neigh_no){
|
fprintf(stderr,"Hej! pri vtx=%d se razikuje stevilo sosedov (prej=%d, potem=%d)\n",i, vtx[i]->neigh_no, tvtx[i]->neigh_no);
|
}
|
}
|
if(eps>dist2);
|
/* We use the temporary vertex for our main vertices and we abandon main
|
* vertices, because their neighbours are not correctly ordered */
|
// tvtx=vlist->vtx;
|
// vlist->vtx=tvtx;
|
// tvlist->vtx=vtx;
|
vtx_list_free(vlist);
|
/* Let's make a check if the number of bonds is correct */
|
if((blist->n)!=3*(tvlist->n-2)){
|
ts_fprintf(stderr,"Number of bonds is %u should be %u!\n", blist->n, 3*(tvlist->n-2));
|
fatal("Number of bonds is not 3*(no_vertex-2).",4);
|
}
|
|
return tvlist;
|
}
|
