From 36387714262ca674e86df01867f0793508ba7d4c Mon Sep 17 00:00:00 2001 From: Sylvain Thery Date: Fri, 26 Apr 2013 17:09:32 +0200 Subject: [PATCH] add import NAS volume files --- include/Algo/Import/import.h | 4 + include/Algo/Import/importMesh.hpp | 3 + include/Algo/Import/importNAS.hpp | 436 +++++++++++++++++++++++++++++ 3 files changed, 443 insertions(+) create mode 100644 include/Algo/Import/importNAS.hpp diff --git a/include/Algo/Import/import.h b/include/Algo/Import/import.h index ffdbc2e46..5ec882c8b 100644 --- a/include/Algo/Import/import.h +++ b/include/Algo/Import/import.h @@ -124,6 +124,9 @@ bool importMSH(typename PFP::MAP& the_map, const std::string& filename, std::vec template bool importVTU(typename PFP::MAP& the_map, const std::string& filename, std::vector& attrNames, float scaleFactor = 1.0f); +template +bool importNAS(typename PFP::MAP& the_map, const std::string& filename, std::vector& attrNames, float scaleFactor = 1.0f); + } // Import @@ -143,6 +146,7 @@ bool importVTU(typename PFP::MAP& the_map, const std::string& filename, std::vec #include "Algo/Import/importNodeEle.hpp" #include "Algo/Import/importMSH.hpp" #include "Algo/Import/importVTU.hpp" +#include "Algo/Import/importNAS.hpp" #include "Algo/Import/importChoupi.hpp" diff --git a/include/Algo/Import/importMesh.hpp b/include/Algo/Import/importMesh.hpp index 4e9e680a0..1fe6cd4b5 100644 --- a/include/Algo/Import/importMesh.hpp +++ b/include/Algo/Import/importMesh.hpp @@ -552,6 +552,9 @@ bool importMesh(typename PFP::MAP& map, const std::string& filename, std::vector return importVTU(map, filename, attrNames, 1.0f); break; + case NAS: + return importNAS(map, filename, attrNames, 1.0f); + break; case OFF: { diff --git a/include/Algo/Import/importNAS.hpp b/include/Algo/Import/importNAS.hpp new file mode 100644 index 000000000..b1b4d300c --- /dev/null +++ b/include/Algo/Import/importNAS.hpp @@ -0,0 +1,436 @@ +/******************************************************************************* +* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps * +* version 0.1 * +* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg * +* * +* This library is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by the * +* Free Software Foundation; either version 2.1 of the License, or (at your * +* option) any later version. * +* * +* This library is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this library; if not, write to the Free Software Foundation, * +* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. * +* * +* Web site: http://cgogn.unistra.fr/ * +* Contact information: cgogn@unistra.fr * +* * +*******************************************************************************/ + +#include "Algo/Modelisation/polyhedron.h" +#include "Geometry/orientation.h" +#include + +namespace CGoGN +{ + +namespace Algo +{ + +namespace Volume +{ + +namespace Import +{ + +float floatFromNas(std::string& s_v) +{ + float x = 0.0f; + + std::size_t pos1 = s_v.find_last_of('-'); + if ((pos1!=std::string::npos) && (pos1!=0)) + { + std::string res = s_v.substr(0,pos1) + "e" + s_v.substr(pos1,8-pos1); + x = atof(res.c_str()); + } + else + { + std::size_t pos2 = s_v.find_last_of('+'); + if ((pos2!=std::string::npos) && (pos2!=0)) + { + std::string res = s_v.substr(0,pos2) + "e" + s_v.substr(pos2,8-pos2); + x = atof(res.c_str()); + } + else + { + x = atof(s_v.c_str()); + } + } + return x; +} + + +template +bool importNAS(typename PFP::MAP& map, const std::string& filename, std::vector& attrNames, float scaleFactor) +{ + typedef typename PFP::VEC3 VEC3; + + VertexAttribute position = map.template addAttribute("position") ; + attrNames.push_back(position.name()) ; + + AttributeContainer& container = map.template getAttributeContainer() ; + + + VertexAutoAttribute< NoMathIONameAttribute< std::vector > > vecDartsPerVertex(map, "incidents"); + + //open file + std::ifstream fp(filename.c_str(), std::ios::in); + if (!fp.good()) + { + CGoGNerr << "Unable to open file " << filename << CGoGNendl; + return false; + } + + std::string ligne; + std::string tag; + + std::getline (fp, ligne); + do + { + std::getline (fp, ligne); + tag = ligne.substr(0,4); + } while (tag !="GRID"); + + unsigned int m_nbVertices = 0; + //reading vertices + std::map verticesMapID; + do + { + std::string s_v = ligne.substr(8,8); + unsigned int ind = atoi(s_v.c_str()); + + s_v = ligne.substr(24,8); + float x = floatFromNas(s_v); + s_v = ligne.substr(32,8); + float y = floatFromNas(s_v); + s_v = ligne.substr(40,8); + float z = floatFromNas(s_v); + + VEC3 pos(x*scaleFactor,y*scaleFactor,z*scaleFactor); + unsigned int id = container.insertLine(); + position[id] = pos; + verticesMapID.insert(std::pair(ind,id)); +// std::cout << "P: "<< ind << " "< tet; + tet.reserve(1000); + std::vector hexa; + tet.reserve(1000); + + do + { + std::string s_v = ligne.substr(0,6); + + if (s_v == "CHEXA ") + { + s_v = ligne.substr(24,8); + unsigned int ind1 = atoi(s_v.c_str()); + s_v = ligne.substr(32,8); + unsigned int ind2 = atoi(s_v.c_str()); + s_v = ligne.substr(40,8); + unsigned int ind3 = atoi(s_v.c_str()); + s_v = ligne.substr(48,8); + unsigned int ind4 = atoi(s_v.c_str()); + s_v = ligne.substr(56,8); + unsigned int ind5 = atoi(s_v.c_str()); + s_v = ligne.substr(64,8); + unsigned int ind6 = atoi(s_v.c_str()); + + std::getline (fp, ligne); + s_v = ligne.substr(8,8); + unsigned int ind7 = atoi(s_v.c_str()); + s_v = ligne.substr(16,8); + unsigned int ind8 = atoi(s_v.c_str()); + + typename PFP::VEC3 P = position[verticesMapID[ind5]]; + typename PFP::VEC3 A = position[verticesMapID[ind1]]; + typename PFP::VEC3 B = position[verticesMapID[ind2]]; + typename PFP::VEC3 C = position[verticesMapID[ind3]]; + + Geom::Vec4ui v; + + if (Geom::testOrientation3D(P,A,B,C) == Geom::OVER) + { + v[0] = ind4; + v[1] = ind3; + v[2] = ind2; + v[3] = ind1; + hexa.push_back(v); + v[0] = ind8; + v[1] = ind7; + v[2] = ind6; + v[3] = ind5; + hexa.push_back(v); + } + else + { + v[0] = ind1; + v[1] = ind2; + v[2] = ind3; + v[3] = ind4; + hexa.push_back(v); + v[0] = ind5; + v[1] = ind6; + v[2] = ind7; + v[3] = ind8; + hexa.push_back(v); + } + } + if (s_v == "CTETRA") + { + s_v = ligne.substr(24,8); + unsigned int ind1 = atoi(s_v.c_str()); + s_v = ligne.substr(32,8); + unsigned int ind2 = atoi(s_v.c_str()); + s_v = ligne.substr(40,8); + unsigned int ind3 = atoi(s_v.c_str()); + s_v = ligne.substr(48,8); + unsigned int ind4 = atoi(s_v.c_str()); + + typename PFP::VEC3 P = position[verticesMapID[ind1]]; + typename PFP::VEC3 A = position[verticesMapID[ind2]]; + typename PFP::VEC3 B = position[verticesMapID[ind3]]; + typename PFP::VEC3 C = position[verticesMapID[ind4]]; + + Geom::Vec4ui v; + + if (Geom::testOrientation3D(P,A,B,C) == Geom::OVER) + { + v[0] = ind4; + v[1] = ind3; + v[2] = ind2; + v[3] = ind1; + } + else + { + v[0] = ind1; + v[1] = ind2; + v[2] = ind3; + v[3] = ind4; + } + tet.push_back(v); + } + std::getline (fp, ligne); + tag = ligne.substr(0,4); + } while (!fp.eof()); + + + + CGoGNout << "nb points = " << m_nbVertices ; + + + unsigned int m_nbVolumes = 0; + + DartMarkerNoUnmark m(map) ; + + if (tet.size() > 0) + { + m_nbVolumes += tet.size(); + + //Read and embed all tetrahedrons + for(unsigned int i = 0; i < tet.size() ; ++i) + { + //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(map,false); + + + // Embed three "base" vertices + for(unsigned int j = 0 ; j < 3 ; ++j) + { + FunctorSetEmb fsetemb(map, verticesMapID[pt[2-j]]); + map.template foreach_dart_of_orbit(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); + } + + //Embed the last "top" vertex + d = map.phi_1(map.phi2(d)); + + FunctorSetEmb fsetemb(map, verticesMapID[pt[3]]); + map.template foreach_dart_of_orbit(d, fsetemb); + + //store darts per vertices to optimize reconstruction + Dart dd = d; + do + { + 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; + + //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(map,false); + + FunctorSetEmb fsetemb(map, verticesMapID[pt[0]]); + + map.template foreach_dart_of_orbit(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); + + d = map.phi1(d); + fsetemb.changeEmb(verticesMapID[pt[1]]); + map.template foreach_dart_of_orbit(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); + + + d = map.phi1(d); + fsetemb.changeEmb(verticesMapID[pt[2]]); + map.template foreach_dart_of_orbit(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); + + + d = map.phi1(d); + fsetemb.changeEmb(verticesMapID[pt[3]]); + map.template foreach_dart_of_orbit(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); + + d = map.template phi<2112>(d); + fsetemb.changeEmb(verticesMapID[ppt[0]]); + + map.template foreach_dart_of_orbit(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); + + d = map.phi_1(d); + fsetemb.changeEmb(verticesMapID[ppt[1]]); + map.template foreach_dart_of_orbit(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); + + d = map.phi_1(d); + fsetemb.changeEmb(verticesMapID[ppt[2]]); + map.template foreach_dart_of_orbit(d, fsetemb); + dd = d; + vecDartsPerVertex[verticesMapID[ppt[2]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd)); + vecDartsPerVertex[verticesMapID[ppt[2]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd)); + vecDartsPerVertex[verticesMapID[ppt[2]]].push_back(dd); m.mark(dd); + + d = map.phi_1(d); + fsetemb.changeEmb(verticesMapID[ppt[3]]); + map.template foreach_dart_of_orbit(d, fsetemb); + dd = d; + vecDartsPerVertex[verticesMapID[ppt[3]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd)); + vecDartsPerVertex[verticesMapID[ppt[3]]].push_back(dd); m.mark(dd); dd = map.phi1(map.phi2(dd)); + vecDartsPerVertex[verticesMapID[ppt[3]]].push_back(dd); m.mark(dd); + + //end of hexa + } + CGoGNout << " / nb hexa = " << hexa.size()/2 << CGoGNendl; + + } + + //Association des phi3 + unsigned int nbBoundaryFaces = 0 ; + for (Dart d = map.begin(); d != map.end(); map.next(d)) + { + if (m.isMarked(d)) + { + std::vector& vec = vecDartsPerVertex[map.phi1(d)]; + + Dart good_dart = NIL; + for(typename std::vector::iterator it = vec.begin(); it != vec.end() && good_dart == NIL; ++it) + { + if(map.template getEmbedding(map.phi1(*it)) == map.template getEmbedding(d) && + map.template getEmbedding(map.phi_1(*it)) == map.template getEmbedding(map.template phi<11>(d))) + { + good_dart = *it ; + } + } + + if (good_dart != NIL) + { + map.sewVolumes(d, good_dart, false); + m.template unmarkOrbit(d); + } + else + { + m.unmarkOrbit(d); + ++nbBoundaryFaces; + } + } + } + + if (nbBoundaryFaces > 0) + { + std::cout << "closing" << std::endl ; + map.closeMap(); + CGoGNout << "Map closed (" << nbBoundaryFaces << " boundary faces)" << CGoGNendl; + } + + fp.close(); + return true; +} + +} // namespace Import + +} + +} // namespace Algo + +} // namespace CGoGN -- GitLab