diff --git a/FileIO/ImportFileTypes.h b/FileIO/ImportFileTypes.h index 6b7eb2ee6eba1f6af6c4915cb88831749fd1c4cc..8f51d734d8de2f8aa75cbf2c710783cc13aba8db 100644 --- a/FileIO/ImportFileTypes.h +++ b/FileIO/ImportFileTypes.h @@ -84,7 +84,7 @@ public: else if (t==ImportFileType::SHAPE) return "ESRI Shape files (*.shp)"; else if (t==ImportFileType::TETGEN) - return "TetGen node files (*.node *.poly)"; + return "TetGen node files (*.node *.poly *.smesh)"; else if (t==ImportFileType::VTK) return "VTK files (*.vtk *.vti *.vtr *.vts *.vtp *.vtu)"; else return "All files (*.*)"; diff --git a/FileIO/TetGenInterface.cpp b/FileIO/TetGenInterface.cpp index ad009c5e93b724430ddd8995a4ba680cfb4eb6e2..ff86aedd9af9c48dc4bebd10e24e75cdd768e388 100644 --- a/FileIO/TetGenInterface.cpp +++ b/FileIO/TetGenInterface.cpp @@ -26,6 +26,9 @@ // FileIO #include "TetGenInterface.h" +// MathLib +#include "TemplatePoint.h" + // MeshLib #include "Mesh.h" #include "Node.h" @@ -35,7 +38,7 @@ namespace FileIO { TetGenInterface::TetGenInterface() : - _zero_based_idx (false) + _zero_based_idx (false), _boundary_markers (false) { } @@ -43,14 +46,20 @@ TetGenInterface::~TetGenInterface() { } -bool TetGenInterface::readTetGenPoly (std::string const& poly_fname, - GeoLib::GEOObjects &geo_objects) +bool TetGenInterface::readTetGenGeometry (std::string const& geo_fname, + GeoLib::GEOObjects &geo_objects) { - std::ifstream poly_stream (poly_fname.c_str()); + std::ifstream poly_stream (geo_fname.c_str()); if (!poly_stream) { - ERR ("TetGenInterface::readTetGenPoly failed to open %s", poly_fname.c_str()); + ERR ("TetGenInterface::readTetGenPoly() failed to open %s", geo_fname.c_str()); + return false; + } + std::string ext (BaseLib::getFileExtension(geo_fname)); + if (ext.compare("poly") && ext.compare("smesh")) + { + ERR ("TetGenInterface::readTetGenPoly() - unknown file type (only *.poly or *.smesh supported)."); return false; } @@ -71,7 +80,12 @@ bool TetGenInterface::readTetGenPoly (std::string const& poly_fname, delete nodes[k]; } std::vector<GeoLib::Surface*> *surfaces = new std::vector<GeoLib::Surface*>; - if (!parseFacets(poly_stream, *surfaces, *points)) + bool failed (true); + if (ext.compare("poly") == 0) + failed = !parsePolyFacets(poly_stream, *surfaces, *points); + else + failed = !parseSmeshFacets(poly_stream, *surfaces, *points); + if (failed) { // remove surfaces read until now but keep the points for (std::size_t k=0; k<surfaces->size(); k++) @@ -80,14 +94,14 @@ bool TetGenInterface::readTetGenPoly (std::string const& poly_fname, surfaces = nullptr; } - std::string geo_name (BaseLib::extractBaseNameWithoutExtension(poly_fname)); + std::string geo_name (BaseLib::extractBaseNameWithoutExtension(geo_fname)); geo_objects.addPointVec(points, geo_name); if (surfaces) geo_objects.addSurfaceVec(surfaces, geo_name); return true; } -std::size_t TetGenInterface::getNFacets(std::ifstream &input) const +std::size_t TetGenInterface::getNFacets(std::ifstream &input) { std::string line; while (!input.fail()) @@ -104,15 +118,18 @@ std::size_t TetGenInterface::getNFacets(std::ifstream &input) const continue; const std::list<std::string> fields = BaseLib::splitString(line, ' '); - return BaseLib::str2number<size_t> (*fields.begin()); - // here this line also includes a flag for boundary markers which we ignore for now + std::list<std::string>::const_iterator it = fields.begin(); + const std::size_t nFacets (BaseLib::str2number<size_t> (*it)); + if (fields.size() > 1) + _boundary_markers = (BaseLib::str2number<size_t> (*(++it)) == 0) ? false : true; + return nFacets; } return false; } -bool TetGenInterface::parseFacets(std::ifstream &input, - std::vector<GeoLib::Surface*> &surfaces, - std::vector<GeoLib::Point*> &points) const +bool TetGenInterface::parsePolyFacets(std::ifstream &input, + std::vector<GeoLib::Surface*> &surfaces, + std::vector<GeoLib::Point*> &points) { const std::size_t nFacets (this->getNFacets(input)); std::size_t nMultPolys (0); @@ -188,6 +205,79 @@ bool TetGenInterface::parseFacets(std::ifstream &input, return false; } +bool TetGenInterface::parseSmeshFacets(std::ifstream &input, + std::vector<GeoLib::Surface*> &surfaces, + std::vector<GeoLib::Point*> &points) +{ + const std::size_t nFacets (this->getNFacets(input)); + std::size_t nMultPolys (0); + std::string line; + surfaces.reserve(nFacets); + std::list<std::string>::const_iterator it; + + const unsigned offset = (_zero_based_idx) ? 0 : 1; + MathLib::TemplatePoint<std::size_t, 3> tri; + std::vector<std::size_t> idx_map; + + for (std::size_t k(0); k<nFacets && !input.fail(); k++) + { + getline (input, line); + if (input.fail()) + { + ERR("TetGenInterface::parseFacets(): Error reading facet %d.", k); + return false; + } + + BaseLib::simplify(line); + if (line.empty() || line.compare(0,1,"#") == 0) + { + k--; + continue; + } + + // read facets + const std::list<std::string> point_fields = BaseLib::splitString(line, ' '); + it = point_fields.begin(); + const std::size_t nPoints = BaseLib::str2number<std::size_t>(*it); + if (nPoints != 3) + { + ERR ("Smesh-files are currently only supported for triangle meshes."); + return false; + } + const std::size_t point_field_size = (_boundary_markers) ? nPoints+1 : nPoints; + if (point_fields.size() > point_field_size) + { + for (std::size_t j(0); j<nPoints; ++j) + tri[j] = (BaseLib::str2number<std::size_t>(*(++it))-offset); + + const std::size_t sfc_marker = (_boundary_markers) ? BaseLib::str2number<std::size_t>(*(++it)) : 0; + const std::size_t idx = std::find(idx_map.begin(), idx_map.end(), sfc_marker) - idx_map.begin(); + if (idx >= surfaces.size()) + { + idx_map.push_back(sfc_marker); + surfaces.push_back(new GeoLib::Surface(points)); + } + surfaces[idx]->addTriangle(tri[0], tri[1], tri[2]); + } + else + { + ERR("TetGenInterface::parseFacets(): Error reading points for facet %d.", k); + return false; + } + } + // here the poly-file potentially defines a number of points to mark holes within the volumes defined by the facets, these are ignored for now + // here the poly-file potentially defines a number of region attributes, these are ignored for now + + std::size_t nTotalTriangles (0); + for (std::size_t i=0; i<surfaces.size(); ++i) + nTotalTriangles += surfaces[i]->getNTriangles(); + if (nTotalTriangles == nFacets) + return true; + + ERR ("TetGenInterface::parseFacets(): Number of expected total triangles (%d) does not match number of found triangles (%d).", surfaces.size(), nTotalTriangles); + return false; +} + MeshLib::Mesh* TetGenInterface::readTetGenMesh (std::string const& nodes_fname, std::string const& ele_fname) { @@ -234,18 +324,16 @@ bool TetGenInterface::readNodesFromStream (std::ifstream &ins, { std::string line; getline (ins, line); - size_t pos_beg (line.find_first_not_of(" ")); size_t n_nodes, dim, n_attributes; bool boundary_markers; while (!ins.fail()) { - line = line.substr(pos_beg); - if (line.compare(0,1,"#") == 0) + BaseLib::simplify(line); + if (line.compare(0,1,"#") == 0 || line.empty()) { // this line is a comment - skip getline (ins, line); - pos_beg = line.find_first_not_of(" "); continue; } // read header line @@ -265,35 +353,25 @@ bool TetGenInterface::parseNodesFileHeader(std::string &line, std::size_t &n_attributes, bool &boundary_markers) const { - std::size_t pos_beg, pos_end; - - // number of nodes - pos_beg = line.find_first_not_of (" "); - pos_end = line.find_first_of(" ", pos_beg); - if (pos_beg != std::string::npos && pos_end != std::string::npos) - n_nodes = BaseLib::str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg)); - else + std::list<std::string> pnt_header = BaseLib::splitString(line, ' '); + if (pnt_header.empty()) { ERR("TetGenInterface::parseNodesFileHeader(): could not read number of nodes specified in header."); return false; } - // dimension - pos_beg = line.find_first_not_of (" ", pos_end); - pos_end = line.find_first_of(" ", pos_beg); - dim = BaseLib::str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg)); - // number of attributes - pos_beg = line.find_first_not_of (" ", pos_end); - pos_end = line.find_first_of(" ", pos_beg); - n_attributes = BaseLib::str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg)); - // boundary marker at nodes? - pos_beg = line.find_first_not_of (" ", pos_end); - pos_end = line.find_first_of(" ", pos_beg); - if (pos_end == std::string::npos) - pos_end = line.size(); - if ((line.substr(pos_beg, pos_end - pos_beg)).compare("1") == 0) - boundary_markers = true; - else + auto it = pnt_header.begin(); + n_nodes = BaseLib::str2number<size_t> (*it); + dim = (pnt_header.size()==1) ? 3 : BaseLib::str2number<size_t> (*(++it)); + + if (pnt_header.size()<4) + { + n_attributes = 0; boundary_markers = false; + return true; + } + + n_attributes = BaseLib::str2number<size_t> (*(++it)); + boundary_markers = ((++it)->compare("1") == 0) ? true : false; return true; } @@ -365,18 +443,16 @@ bool TetGenInterface::readElementsFromStream(std::ifstream &ins, { std::string line; getline (ins, line); - std::size_t pos_beg (line.find_first_not_of(" ")); std::size_t n_tets, n_nodes_per_tet; bool region_attributes; while (!ins.fail()) { - line = line.substr(pos_beg); - if (line.compare(0,1,"#") == 0) + BaseLib::simplify(line); + if (line.compare(0,1,"#") == 0 || line.empty()) { // this line is a comment - skip getline (ins, line); - pos_beg = line.find_first_not_of(" "); continue; } @@ -500,9 +576,9 @@ bool TetGenInterface::parseElements(std::ifstream& ins, return true; } -bool TetGenInterface::writeTetGenPoly(const std::string &file_name, - const GeoLib::GEOObjects &geo_objects, - const std::string &geo_name) const +bool TetGenInterface::writeTetGenSmesh(const std::string &file_name, + const GeoLib::GEOObjects &geo_objects, + const std::string &geo_name) const { std::vector<GeoLib::Point*> const*const points = geo_objects.getPointVec(geo_name); std::vector<GeoLib::Surface*> const*const surfaces = geo_objects.getSurfaceVec(geo_name); @@ -518,19 +594,20 @@ bool TetGenInterface::writeTetGenPoly(const std::string &file_name, std::ofstream out( file_name.c_str(), std::ios::out ); // the points header const std::size_t nPoints (points->size()); - out << nPoints << " 3 0 0\n"; + out << nPoints << " 3\n"; // the point list for (std::size_t i=0; i<nPoints; ++i) out << i << " " << (*(*points)[i])[0] << " " << (*(*points)[i])[1] << " " << (*(*points)[i])[2] << "\n"; // the surfaces header const std::size_t nSurfaces = (surfaces) ? surfaces->size() : 0; - out << nSurfaces << " 0\n"; - // the facets list + std::size_t nTotalTriangles (0); + for (std::size_t i=0; i<nSurfaces; ++i) + nTotalTriangles += (*surfaces)[i]->getNTriangles(); + out << nSurfaces << " 1\n"; + for (std::size_t i=0; i<nSurfaces; ++i) { - // the number of polys per facet const std::size_t nTriangles ((*surfaces)[i]->getNTriangles()); - out << nTriangles << "\n"; // the poly list for (std::size_t j=0; j<nTriangles; ++j) { diff --git a/FileIO/TetGenInterface.h b/FileIO/TetGenInterface.h index 72ca2c04618a2dabfe7e2134e9633a2879395eea..38b677277e7fa7f66c58ad8b3bfa2860df270889 100644 --- a/FileIO/TetGenInterface.h +++ b/FileIO/TetGenInterface.h @@ -43,13 +43,13 @@ public: virtual ~TetGenInterface(); /** - * Method reads the TetGen mesh from node file and element file. - * @param poly_fname file name of the poly file + * Method reads geometry from a TetGen poly or smesh file. + * @param geo_fname file name of the poly file * @param geo_objects where the geometry is written to * @return on success the method returns true, otherwise it returns false */ - bool readTetGenPoly (std::string const& poly_fname, - GeoLib::GEOObjects &geo_objects); + bool readTetGenGeometry (std::string const& geo_fname, + GeoLib::GEOObjects &geo_objects); /** * Method reads the TetGen mesh from node file and element file. @@ -61,19 +61,19 @@ public: std::string const& ele_fname); /** - * Writes the geometry of a given name to TetGen poly-file. - * @param file_name file name of the new poly file + * Writes the geometry of a given name to TetGen smesh-file. + * @param file_name file name of the new smesh file * @param geo_objects the container for the geometry. * @param geo_name the name for the geometry. * @return returns true on success and false otherwise. */ - bool writeTetGenPoly(const std::string &file_name, - const GeoLib::GEOObjects &geo_objects, - const std::string &geo_name) const; + bool writeTetGenSmesh(const std::string &file_name, + const GeoLib::GEOObjects &geo_objects, + const std::string &geo_name) const; private: /// Returns the declared number of facets in the poly file. - std::size_t getNFacets(std::ifstream &input) const; + std::size_t getNFacets(std::ifstream &input); /** * Method parses the lines reading the facets from TetGen poly file @@ -82,9 +82,20 @@ private: * @param points the point vector needed for creating surfaces (input) * @return true, if the facets have been read correctly, false if the method detects an error */ - bool parseFacets(std::ifstream &input, - std::vector<GeoLib::Surface*> &surfaces, - std::vector<GeoLib::Point*> &points) const; + bool parsePolyFacets(std::ifstream &input, + std::vector<GeoLib::Surface*> &surfaces, + std::vector<GeoLib::Point*> &points); + + /** + * Method parses the lines reading the facets from TetGen smesh file + * @param input the input stream (input) + * @param surfaces the vector of surfaces to be filled (output) + * @param points the point vector needed for creating surfaces (input) + * @return true, if the facets have been read correctly, false if the method detects an error + */ + bool parseSmeshFacets(std::ifstream &input, + std::vector<GeoLib::Surface*> &surfaces, + std::vector<GeoLib::Point*> &points); /** * Method reads the nodes from stream and stores them in a node vector. @@ -162,10 +173,12 @@ private: std::size_t n_nodes_per_tet, bool region_attribute) const; - /** - * the value is true if the indexing is zero based, else false - */ + + /// the value is true if the indexing is zero based, else false bool _zero_based_idx; + + /// true if boundary markers are set, false otherwise + bool _boundary_markers; }; }