Skip to content

GitLab

C
CGoGN
Commit a729c433 authored by Sylvain Thery's avatar Sylvain Thery
Browse files
Options

import VTU ok

parent 45f74840
Hide whitespace changes
Inline Side-by-side
Showing with 216 additions and 345 deletions
include/Algo/Import/importVTU.hpp
View file @ a729c433
......@@ -53,455 +53,327 @@ bool importVTU(typename PFP::MAP& map, const std::string& filename, std::vector<
AttributeContainer& container = map.template getAttributeContainer<VERTEX>() ;
unsigned int m_nbVertices = 0, m_nbVolumes = 0;
VertexAutoAttribute< NoMathIONameAttribute< std::vector<Dart> > > vecDartsPerVertex(map, "incidents");
xmlDocPtr doc = xmlReadFile(filename.c_str(), NULL, 0);
xmlNodePtr vtu_node = xmlDocGetRootElement(doc);
std::cout << " NAME "<<vtu_node->name << std::endl;
// std::cout << " NAME "<<vtu_node->name << std::endl;
xmlChar *prop = xmlGetProp(vtu_node, BAD_CAST "type");
std::cout << "type = "<< prop << std::endl;
// std::cout << "type = "<< prop << std::endl;
xmlNode* cur_node = vtu_node->children;
while (strcmp((char*)(cur_node->name),(char*)"UnstructuredGrid")!=0)
cur_node = cur_node->next;
xmlNode* grid_node = vtu_node->children;
while (strcmp((char*)(grid_node->name),(char*)"UnstructuredGrid")!=0)
grid_node = grid_node->next;
cur_node = cur_node->children;
while (strcmp((char*)(cur_node->name),(char*)"Piece")!=0)
cur_node = cur_node->next;
xmlNode* piece_node = grid_node->children;
while (strcmp((char*)(piece_node->name),(char*)"Piece")!=0)
piece_node = piece_node->next;
prop = xmlGetProp(cur_node, BAD_CAST "NumberOfPoints");
unsigned int nbp = atoi((char*)(prop));
prop = xmlGetProp(piece_node, BAD_CAST "NumberOfPoints");
unsigned int nbVertices = atoi((char*)(prop));
prop = xmlGetProp(cur_node, BAD_CAST "NumberOfCells");
unsigned int nbc = atoi((char*)(prop));
prop = xmlGetProp(piece_node, BAD_CAST "NumberOfCells");
unsigned int nbVolumes = atoi((char*)(prop));
std::cout << "Number of points = "<< nbp<< std::endl;
std::cout << "Number of cells = "<< nbc<< std::endl;
std::cout << "Number of points = "<< nbVertices<< std::endl;
std::cout << "Number of cells = "<< nbVolumes << std::endl;
cur_node = cur_node->children;
while (strcmp((char*)(cur_node->name),(char*)"Points")!=0)
cur_node = cur_node->next;
xmlNode* points_node = piece_node->children;
while (strcmp((char*)(points_node->name),(char*)"Points")!=0)
points_node = points_node->next;
cur_node = cur_node->children;
while (strcmp((char*)(cur_node->name),(char*)"DataArray")!=0)
cur_node = cur_node->next;
points_node = points_node->children;
while (strcmp((char*)(points_node->name),(char*)"DataArray")!=0)
points_node = points_node->next;
std::stringstream ss((char*)(xmlNodeGetContent(cur_node->children)));
std::vector<unsigned int> verticesID;
verticesID.reserve(nbVertices);
for (unsigned int i=0; i< nbp; ++i)
std::stringstream ss((char*)(xmlNodeGetContent(points_node->children)));
for (unsigned int i=0; i< nbVertices; ++i)
{
typename PFP::VEC3 P;
ss >> P[0]; ss >> P[1]; ss >> P[2];
std::cout << P << std::endl;
// std::cout << P << std::endl;
unsigned int id = container.insertLine();
position[id] = P;
verticesID.push_back(id);
}
// for (xmlNode* x_node = cur_node->children; x_node!=NULL; x_node = x_node->next)
// {
// // float f = atof((char*)(xmlNodeGetContent(x_node)));
// std::cout << xmlNodeGetContent(x_node) << std::endl;
// }
// for (xmlNode* cur_node = vtu_node->children; cur_node; cur_node = cur_node->next)
// {
// std::cout << " NAME "<<cur_node->name << std::endl;
// }
// cur_node = cur_node->children;
// std::cout << " NAME "<<cur_node->name << std::endl;
xmlNode* cell_node = piece_node->children;
while (strcmp((char*)(cell_node->name),(char*)"Cells")!=0)
cell_node = cell_node->next;
xmlFreeDoc(doc);
return true;
std::cout <<"CELL NODE = "<< cell_node->name << std::endl;
/*
if (strcmp((char*)(map_node->name),(char*)"CGoGN_Map")!=0)
{
std::cerr << "Wrong xml format: Root node != CGoGN_Map"<< std::endl;
return false;
}
*/
// cell_node = cell_node->children;
// while (strcmp((char*)(cell_node->name),(char*)"DataArray")!=0)
// cell_node = cell_node->next;
//open file
std::ifstream fp(filename.c_str(), std::ios::in);
if (!fp.good())
{
CGoGNerr << "Unable to open file " << filename << CGoGNendl;
return false;
}
std::vector<unsigned char> typeVols;
typeVols.reserve(nbVolumes);
std::vector<unsigned int> offsets;
offsets.reserve(nbVolumes);
std::vector<unsigned int> indices;
indices.reserve(nbVolumes*4);
std::string ligne;
unsigned int nbv=0;
//read $NODE
std::getline (fp, ligne);
// reading number of vertices
std::getline (fp, ligne);
std::stringstream oss(ligne);
oss >> nbv;
// cell_node = cell_node->children;
// while (strcmp((char*)(cell_node->name),(char*)"DataArray")!=0)
// cell_node = cell_node->next;
// std::cout <<"1NODE = "<< (char*)(cell_node->name) << std::endl;
// xmlChar* type = xmlGetProp(cell_node, BAD_CAST "Name");
// std::cout <<"1NAME = "<< (char*)(type) << std::endl;
// while (strcmp((char*)(cell_node->name),(char*)"DataArray")!=0)
// cell_node = cell_node->next;
// std::cout <<"2NODE = "<< (char*)(cell_node->name) << std::endl;
// type = xmlGetProp(cell_node, BAD_CAST "Name");
// std::cout <<"2NAME = "<< (char*)(type) << std::endl;
//reading vertices
// std::vector<unsigned int> verticesID;
std::map<unsigned int, unsigned int> verticesMapID;
// while (strcmp((char*)(cell_node->name),(char*)"DataArray")!=0)
// cell_node = cell_node->next;
// std::cout <<"3NODE = "<< (char*)(cell_node->name) << std::endl;
// type = xmlGetProp(cell_node, BAD_CAST "Name");
// std::cout <<"3NAME = "<< (char*)(type) << std::endl;
// verticesID.reserve(nbv);
for(unsigned int i = 0; i < nbv;++i)
for (xmlNode* x_node = cell_node->children; x_node!=NULL; x_node = x_node->next)
{
do
{
std::getline (fp, ligne);
} while (ligne.size() == 0);
std::stringstream oss(ligne);
unsigned int pipo;
float x,y,z;
oss >> pipo;
oss >> x;
oss >> y;
oss >> z;
// TODO : if required read other vertices attributes here
VEC3 pos(x*scaleFactor,y*scaleFactor,z*scaleFactor);
unsigned int id = container.insertLine();
position[id] = pos;
verticesMapID.insert(std::pair<unsigned int, unsigned int>(pipo,id));
// verticesID.push_back(id);
}
// ENNODE
std::getline (fp, ligne);
m_nbVertices = nbv;
while ((x_node!=NULL) && (strcmp((char*)(x_node->name),(char*)"DataArray")!=0))
x_node = x_node->next;
// ELM
std::getline (fp, ligne);
// reading number of elements
std::getline (fp, ligne);
unsigned int nbe=0;
std::stringstream oss2(ligne);
oss2 >> nbe;
std::vector<Geom::Vec4ui> tet;
tet.reserve(1000);
std::vector<Geom::Vec4ui> hexa;
tet.reserve(1000);
bool invertVol = false;
for(unsigned int i=0; i<nbe; ++i)
{
unsigned int pipo,type_elm,nb;
fp >> pipo;
fp >> type_elm;
fp >> pipo;
fp >> pipo;
fp >> nb;
if ((type_elm==4) && (nb==4))
if (x_node == NULL)
break;
else
{
Geom::Vec4ui v;
xmlChar* type = xmlGetProp(x_node, BAD_CAST "Name");
// std::cout <<"NODE = "<< (char*)(x_node->name) << std::endl;
// std::cout <<"NAME = "<< (char*)(type) << std::endl;
// test orientation of first tetra
if (i==0)
if (strcmp((char*)(type),(char*)"connectivity")==0)
{
fp >> v[0];
fp >> v[1];
fp >> v[2];
fp >> v[3];
typename PFP::VEC3 P = position[verticesMapID[v[0]]];
typename PFP::VEC3 A = position[verticesMapID[v[1]]];
typename PFP::VEC3 B = position[verticesMapID[v[2]]];
typename PFP::VEC3 C = position[verticesMapID[v[3]]];
if (Geom::testOrientation3D<typename PFP::VEC3>(P,A,B,C) == Geom::OVER)
// c_indices =(char*)(xmlNodeGetContent(x_node->children));
std::stringstream ss((char*)(xmlNodeGetContent(x_node->children)));
while (!ss.eof())
{
invertVol=true;
unsigned int ui=v[0];
v[0] = v[3];
v[3] = v[2];
v[2] = v[1];
v[1] = ui;
unsigned int ind;
ss >> ind;
indices.push_back(ind);
}
}
else
if (strcmp((char*)(type),(char*)"offsets")==0)
{
if (invertVol)
{
fp >> v[1];
fp >> v[2];
fp >> v[3];
fp >> v[0];
}
else
std::stringstream ss((char*)(xmlNodeGetContent(x_node->children)));
for (unsigned int i=0; i< nbVolumes; ++i)
{
fp >> v[0];
fp >> v[1];
fp >> v[2];
fp >> v[3];
unsigned int o;
ss >> o;
offsets.push_back(o);
}
}
tet.push_back(v);
}
else
{
if ((type_elm==5) && (nb==8))
if (strcmp((char*)(type),(char*)"types")==0)
{
std::cout << "HEXA: "<< i << std::endl;
Geom::Vec4ui v;
if (i==0)
bool unsupported = false;
std::stringstream ss((char*)(xmlNodeGetContent(x_node->children)));
for (unsigned int i=0; i< nbVolumes; ++i)
{
unsigned int last;
fp >> v[0];
fp >> v[1];
fp >> v[2];
fp >> v[3];
fp >> last;
typename PFP::VEC3 P = position[verticesMapID[last]];
typename PFP::VEC3 A = position[verticesMapID[v[0]]];
typename PFP::VEC3 B = position[verticesMapID[v[1]]];
typename PFP::VEC3 C = position[verticesMapID[v[2]]];
if (Geom::testOrientation3D<typename PFP::VEC3>(P,A,B,C) == Geom::OVER)
unsigned int t;
ss >> t;
if ((t != 12) && (t!= 10))
{
invertVol=true;
unsigned int ui = v[3];
v[3] = v[0];
v[0] = ui;
ui = v[2];
v[2] = v[1];
v[1] = ui;
hexa.push_back(v);
v[3] = last;
fp >> v[2];
fp >> v[1];
fp >> v[0];
hexa.push_back(v);
unsupported = true;
typeVols.push_back(0);
}
else
{
hexa.push_back(v);
v[0] = last;
fp >> v[1];
fp >> v[2];
fp >> v[3];
hexa.push_back(v);
}
}
else
{
if (invertVol)
{
fp >> v[3];
fp >> v[2];
fp >> v[1];
fp >> v[0];
hexa.push_back(v);
fp >> v[3];
fp >> v[2];
fp >> v[1];
fp >> v[0];
hexa.push_back(v);
}
else
{
fp >> v[0];
fp >> v[1];
fp >> v[2];
fp >> v[3];
hexa.push_back(v);
fp >> v[0];
fp >> v[1];
fp >> v[2];
fp >> v[3];
hexa.push_back(v);
typeVols.push_back((unsigned char)t);
}
}
}
else
{
for (unsigned int j=0; j<nb; ++j)
{
unsigned int v;
fp >> v;
}
if (unsupported)
CGoGNerr << "warning, some unsupported volume cell types"<< CGoGNendl;
}
}
}
CGoGNout << "nb points = " << m_nbVertices ;
m_nbVolumes = 0;
xmlFreeDoc(doc);
DartMarkerNoUnmark m(map) ;
if (tet.size() > 0)
unsigned int currentOffset = 0;
for (unsigned int i=0; i< nbVolumes; ++i)
{
m_nbVolumes += tet.size();
//Read and embed all tetrahedrons
for(unsigned int i = 0; i < tet.size() ; ++i)
if (typeVols[i]==12)
{
//start one tetra
// Geom::Vec4ui& pt = tet[i];
Geom::Vec4ui pt;
pt[0] = tet[i][0];
pt[1] = tet[i][1];
pt[2] = tet[i][2];
pt[3] = tet[i][3];
Dart d = Surface::Modelisation::createTetrahedron<PFP>(map,false);
Dart d = Surface::Modelisation::createHexahedron<PFP>(map,false);
// Embed three "base" vertices
for(unsigned int j = 0 ; j < 3 ; ++j)
unsigned int pt[8];
pt[0] = indices[currentOffset];
pt[1] = indices[currentOffset+1];
pt[2] = indices[currentOffset+2];
pt[3] = indices[currentOffset+3];
pt[4] = indices[currentOffset+4];
typename PFP::VEC3 P = position[verticesID[indices[currentOffset+4]]];
typename PFP::VEC3 A = position[verticesID[indices[currentOffset ]]];
typename PFP::VEC3 B = position[verticesID[indices[currentOffset+1]]];
typename PFP::VEC3 C = position[verticesID[indices[currentOffset+2]]];
if (Geom::testOrientation3D<typename PFP::VEC3>(P,A,B,C) == Geom::OVER)
{
FunctorSetEmb<typename PFP::MAP, VERTEX> fsetemb(map, verticesMapID[pt[2-j]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
//store darts per vertices to optimize reconstruction
Dart dd = d;
do
{
m.mark(dd) ;
vecDartsPerVertex[verticesMapID[pt[2-j]]].push_back(dd);
dd = map.phi1(map.phi2(dd));
} while(dd != d);
d = map.phi1(d);
pt[0] = indices[currentOffset+3];
pt[1] = indices[currentOffset+2];
pt[2] = indices[currentOffset+1];
pt[3] = indices[currentOffset+0];
pt[4] = indices[currentOffset+7];
pt[5] = indices[currentOffset+6];
pt[6] = indices[currentOffset+5];
pt[7] = indices[currentOffset+4];
}
//Embed the last "top" vertex
d = map.phi_1(map.phi2(d));
FunctorSetEmb<typename PFP::MAP, VERTEX> fsetemb(map, verticesMapID[pt[3]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
//store darts per vertices to optimize reconstruction
Dart dd = d;
do
else
{
m.mark(dd) ;
vecDartsPerVertex[verticesMapID[pt[3]]].push_back(dd);
dd = map.phi1(map.phi2(dd));
} while(dd != d);
//end of tetra
}
CGoGNout << " / nb tetra = " << tet.size() << CGoGNendl;
}
if (hexa.size() > 0)
{
m_nbVolumes += hexa.size()/2;
pt[0] = indices[currentOffset+0];
pt[1] = indices[currentOffset+1];
pt[2] = indices[currentOffset+2];
pt[3] = indices[currentOffset+3];
pt[4] = indices[currentOffset+4];
pt[5] = indices[currentOffset+5];
pt[6] = indices[currentOffset+6];
pt[7] = indices[currentOffset+7];
}
//Read and embed all tetrahedrons
for(unsigned int i = 0; i < hexa.size()/2 ; ++i)
{
// one hexa
Geom::Vec4ui pt;
pt[0] = hexa[2*i][0];
pt[1] = hexa[2*i][1];
pt[2] = hexa[2*i][2];
pt[3] = hexa[2*i][3];
Geom::Vec4ui ppt;
ppt[0] = hexa[2*i+1][0];
ppt[1] = hexa[2*i+1][1];
ppt[2] = hexa[2*i+1][2];
ppt[3] = hexa[2*i+1][3];
Dart d = Surface::Modelisation::createHexahedron<PFP>(map,false);
FunctorSetEmb<typename PFP::MAP, VERTEX> fsetemb(map, verticesMapID[pt[0]]);
FunctorSetEmb<typename PFP::MAP, VERTEX> fsetemb(map, verticesID[pt[0]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
Dart dd = d;
vecDartsPerVertex[verticesMapID[pt[0]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[0]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[0]]].push_back(dd); m.mark(dd);
vecDartsPerVertex[verticesID[pt[0]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[0]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[0]]].push_back(dd); m.mark(dd);
d = map.phi1(d);
fsetemb.changeEmb(verticesMapID[pt[1]]);
fsetemb.changeEmb(verticesID[pt[1]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
dd = d;
vecDartsPerVertex[verticesMapID[pt[1]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[1]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[1]]].push_back(dd); m.mark(dd);
vecDartsPerVertex[verticesID[pt[1]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[1]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[1]]].push_back(dd); m.mark(dd);
d = map.phi1(d);
fsetemb.changeEmb(verticesMapID[pt[2]]);
fsetemb.changeEmb(verticesID[pt[2]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
dd = d;
vecDartsPerVertex[verticesMapID[pt[2]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[2]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[2]]].push_back(dd); m.mark(dd);
vecDartsPerVertex[verticesID[pt[2]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[2]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[2]]].push_back(dd); m.mark(dd);
d = map.phi1(d);
fsetemb.changeEmb(verticesMapID[pt[3]]);
fsetemb.changeEmb(verticesID[pt[3]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
dd = d;
vecDartsPerVertex[verticesMapID[pt[3]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[3]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[pt[3]]].push_back(dd); m.mark(dd);
vecDartsPerVertex[verticesID[pt[3]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[3]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[3]]].push_back(dd); m.mark(dd);
d = map.template phi<2112>(d);
fsetemb.changeEmb(verticesMapID[ppt[0]]);
fsetemb.changeEmb(verticesID[pt[4]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
dd = d;
vecDartsPerVertex[verticesMapID[ppt[0]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[ppt[0]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[ppt[0]]].push_back(dd); m.mark(dd);
vecDartsPerVertex[verticesID[pt[4]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[4]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[4]]].push_back(dd); m.mark(dd);
d = map.phi_1(d);
fsetemb.changeEmb(verticesMapID[ppt[1]]);
fsetemb.changeEmb(verticesID[pt[5]]);
map.template foreach_dart_of_orbit<PFP::MAP::VERTEX_OF_PARENT>(d, fsetemb);
dd = d;
vecDartsPerVertex[verticesMapID[ppt[1]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[ppt[1]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesMapID[ppt[1]]].push_back(dd); m.mark(dd);
vecDartsPerVertex[verticesID[pt[5]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[5]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd));
vecDartsPerVertex[verticesID[pt[5]]].push_back(dd); m.mark(dd);
d = map.phi_1(d);
fsetemb.changeEmb(verticesMapID[ppt[2]]);
fsetemb.changeEmb(verticesID[pt[6]]);
ma