38 std::ifstream poly_stream(geo_fname.c_str());
42 ERR(
"TetGenInterface::readTetGenGeometry() failed to open {:s}",
49 ERR(
"TetGenInterface::readTetGenGeometry() - unknown file type (only "
50 "*.smesh is supported).");
54 std::vector<MeshLib::Node*> nodes;
58 for (
auto& node : nodes)
64 const std::size_t nNodes(nodes.size());
65 auto points = std::make_unique<std::vector<GeoLib::Point*>>();
66 points->reserve(nNodes);
67 for (std::size_t k(0); k < nNodes; ++k)
69 points->push_back(
new GeoLib::Point(*(nodes[k]), nodes[k]->getID()));
73 geo_objects.
addPointVec(std::move(points), geo_name);
74 const std::vector<std::size_t>& id_map(
77 auto surfaces = std::make_unique<std::vector<GeoLib::Surface*>>();
79 poly_stream, *surfaces, *geo_objects.
getPointVec(geo_name), id_map))
82 for (std::size_t k = 0; k < surfaces->size(); k++)
84 delete (*surfaces)[k];
100 ERR(
"TetGenInterface::getNFacets(): Error reading number of "
106 if (line.empty() || line.compare(0, 1,
"#") == 0)
112 auto it = fields.begin();
113 const auto nFacets(BaseLib::str2number<std::size_t>(*it));
114 if (fields.size() > 1)
124 std::ifstream& input,
125 std::vector<GeoLib::Surface*>& surfaces,
126 const std::vector<GeoLib::Point*>& points,
127 const std::vector<std::size_t>& pnt_id_map)
129 const std::size_t nFacets(this->
getNFacets(input));
131 surfaces.reserve(nFacets);
132 std::list<std::string>::const_iterator it;
135 std::vector<std::size_t> idx_map;
138 while (k < nFacets && !input.fail())
140 getline(input, line);
143 ERR(
"TetGenInterface::parseFacets(): Error reading facet {:d}.", k);
148 if (line.empty() || line.compare(0, 1,
"#") == 0)
154 const std::list<std::string> point_fields =
156 it = point_fields.begin();
157 const auto nPoints = BaseLib::str2number<std::size_t>(*it);
160 ERR(
"Smesh-files are currently only supported for triangle "
164 std::vector<std::size_t> point_ids;
165 const std::size_t point_field_size =
167 if (point_fields.size() > point_field_size)
169 for (std::size_t j(0); j < nPoints; ++j)
172 pnt_id_map[BaseLib::str2number<std::size_t>(*(++it)) -
176 const std::size_t sfc_marker =
179 const std::size_t idx =
180 std::find(idx_map.begin(), idx_map.end(), sfc_marker) -
182 if (idx >= surfaces.size())
184 idx_map.push_back(sfc_marker);
187 surfaces[idx]->addTriangle(
188 point_ids[0], point_ids[1], point_ids[2]);
192 ERR(
"TetGenInterface::parseFacets(): Error reading points for "
204 std::size_t nTotalTriangles(0);
205 for (
auto& surface : surfaces)
207 nTotalTriangles += surface->getNumberOfTriangles();
209 if (nTotalTriangles == nFacets)
214 ERR(
"TetGenInterface::parseFacets(): Number of expected total triangles "
215 "({:d}) does not match number of found triangles ({:d}).",
222 std::string
const& ele_fname)
224 std::ifstream ins_nodes(nodes_fname.c_str());
225 std::ifstream ins_ele(ele_fname.c_str());
227 if (!ins_nodes || !ins_ele)
231 ERR(
"TetGenInterface::readTetGenMesh failed to open {:s}",
236 ERR(
"TetGenInterface::readTetGenMesh failed to open {:s}",
242 std::vector<MeshLib::Node*> nodes;
246 for (
auto& node : nodes)
253 std::vector<MeshLib::Element*> elements;
254 std::vector<int> materials;
264 materials.cbegin(), materials.cend(), [](
int m) { return m != 0; }))
268 mat_props->reserve(elements.size());
271 std::back_inserter(*mat_props));
274 const std::string mesh_name(
276 return new MeshLib::Mesh(mesh_name, nodes, elements, properties);
280 std::vector<MeshLib::Node*>& nodes)
286 std::size_t n_attributes;
287 bool boundary_markers;
292 if (line.empty() || line.compare(0, 1,
"#") == 0)
300 line, n_nodes, dim, n_attributes, boundary_markers);
315 std::size_t& n_nodes,
317 std::size_t& n_attributes,
318 bool& boundary_markers)
321 if (pnt_header.empty())
323 ERR(
"TetGenInterface::parseNodesFileHeader(): could not read number of "
324 "nodes specified in header.");
327 auto it = pnt_header.begin();
328 n_nodes = BaseLib::str2number<std::size_t>(*it);
329 dim = (pnt_header.size() == 1) ? 3
330 : BaseLib::str2number<std::size_t>(*(++it));
332 if (pnt_header.size() < 4)
335 boundary_markers =
false;
339 n_attributes = BaseLib::str2number<std::size_t>(*(++it));
340 boundary_markers = *(++it) ==
"1";
346 std::vector<MeshLib::Node*>& nodes,
351 nodes.reserve(n_nodes);
354 while (k < n_nodes && !ins.fail())
356 std::vector<double> coordinates(dim);
357 std::getline(ins, line);
360 ERR(
"TetGenInterface::parseNodes(): Error reading node {:d}.", k);
365 std::size_t pos_end = 0;
366 std::size_t pos_beg = line.find_first_not_of(
' ', pos_end);
367 pos_end = line.find_first_of(
" \n", pos_beg);
369 if (line.empty() || pos_beg == pos_end ||
370 line.compare(pos_beg, 1,
"#") == 0)
375 if (pos_beg != std::string::npos && pos_end != std::string::npos)
377 id = BaseLib::str2number<std::size_t>(
378 line.substr(pos_beg, pos_end - pos_beg));
379 if (k == 0 &&
id == 0)
386 ERR(
"TetGenInterface::parseNodes(): Error reading ID of node {:d}.",
392 for (std::size_t i(0); i < dim; i++)
394 pos_beg = line.find_first_not_of(
' ', pos_end);
395 pos_end = line.find_first_of(
" \n", pos_beg);
396 if (pos_end == std::string::npos)
398 pos_end = line.size();
400 if (pos_beg != std::string::npos)
402 coordinates[i] = BaseLib::str2number<double>(
403 line.substr(pos_beg, pos_end - pos_beg));
407 ERR(
"TetGenInterface::parseNodes(): error reading coordinate "
408 "{:d} of node {:d}.",
415 nodes.push_back(
new MeshLib::Node(coordinates.data(),
id - offset));
426 std::vector<MeshLib::Element*>& elements,
427 std::vector<int>& materials,
428 const std::vector<MeshLib::Node*>& nodes)
const
431 std::getline(ins, line);
433 std::size_t n_nodes_per_tet;
434 bool region_attributes;
439 if (line.empty() || line.compare(0, 1,
"#") == 0)
442 std::getline(ins, line);
448 line, n_tets, n_nodes_per_tet, region_attributes);
470 std::size_t& n_nodes_per_tet,
471 bool& region_attribute)
477 pos_beg = line.find_first_not_of(
' ');
478 pos_end = line.find_first_of(
' ', pos_beg);
479 if (pos_beg != std::string::npos && pos_end != std::string::npos)
481 n_tets = BaseLib::str2number<std::size_t>(
482 line.substr(pos_beg, pos_end - pos_beg));
486 ERR(
"TetGenInterface::parseElementsFileHeader(): Could not read number "
487 "of tetrahedra specified in header.");
491 pos_beg = line.find_first_not_of(
" \t", pos_end);
492 pos_end = line.find_first_of(
" \t", pos_beg);
493 n_nodes_per_tet = BaseLib::str2number<std::size_t>(
494 line.substr(pos_beg, pos_end - pos_beg));
496 pos_beg = line.find_first_not_of(
" \t", pos_end);
497 pos_end = line.find_first_of(
" \t\n", pos_beg);
498 if (pos_end == std::string::npos)
500 pos_end = line.size();
502 region_attribute = line.substr(pos_beg, pos_end - pos_beg) ==
"1";
508 std::vector<MeshLib::Element*>& elements,
509 std::vector<int>& materials,
510 const std::vector<MeshLib::Node*>& nodes,
512 std::size_t n_nodes_per_tet,
513 bool region_attribute)
const
516 std::vector<std::size_t> ids(n_nodes_per_tet);
518 elements.reserve(n_tets);
519 materials.reserve(n_tets);
522 for (std::size_t k(0); k < n_tets && !ins.fail(); k++)
527 ERR(
"TetGenInterface::parseElements(): Error reading tetrahedron "
533 std::size_t pos_end = 0;
534 std::size_t pos_beg = line.find_first_not_of(
' ', pos_end);
535 pos_end = line.find_first_of(
" \n", pos_beg);
537 if (line.empty() || pos_beg == pos_end ||
538 line.compare(pos_beg, 1,
"#") == 0)
544 if (pos_beg == std::string::npos || pos_end == std::string::npos)
546 ERR(
"TetGenInterface::parseElements(): Error reading id of "
553 for (std::size_t i(0); i < n_nodes_per_tet; i++)
555 pos_beg = line.find_first_not_of(
' ', pos_end);
556 pos_end = line.find_first_of(
' ', pos_beg);
557 if (pos_end == std::string::npos)
559 pos_end = line.size();
561 if (pos_beg != std::string::npos && pos_end != std::string::npos)
563 ids[i] = BaseLib::str2number<std::size_t>(
564 line.substr(pos_beg, pos_end - pos_beg)) -
569 ERR(
"TetGenInterface::parseElements(): Error reading node {:d} "
570 "of tetrahedron {:d}.",
579 if (region_attribute)
581 pos_beg = line.find_first_not_of(
' ', pos_end);
582 pos_end = line.find_first_of(
' ', pos_beg);
583 if (pos_end == std::string::npos)
585 pos_end = line.size();
587 if (pos_beg != std::string::npos && pos_end != std::string::npos)
589 region = BaseLib::str2number<int>(
590 line.substr(pos_beg, pos_end - pos_beg));
594 ERR(
"TetGenInterface::parseElements(): Error reading region "
595 "attribute of tetrahedron {:d}.",
602 for (
int n = 0; n < 4; n++)
604 tet_nodes[n] = nodes[ids[n]];
607 materials.push_back(region);
614 const std::string& file_name,
616 const std::string& geo_name,
617 const std::vector<GeoLib::Point>& attribute_points)
619 std::vector<GeoLib::Point*>
const*
const points =
621 std::vector<GeoLib::Surface*>
const*
const surfaces =
624 if (points ==
nullptr)
626 ERR(
"Geometry {:s} not found.", geo_name);
629 if (surfaces ==
nullptr)
631 WARN(
"No surfaces found for geometry {:s}. Writing points only.",
635 std::ofstream out(file_name.c_str(), std::ios::out);
636 out.precision(std::numeric_limits<double>::digits10);
638 const std::size_t nPoints(points->size());
639 out << nPoints <<
" 3\n";
641 for (std::size_t i = 0; i < nPoints; ++i)
643 out << i <<
" " << (*(*points)[i])[0] <<
" " << (*(*points)[i])[1]
644 <<
" " << (*(*points)[i])[2] <<
"\n";
647 const std::size_t nSurfaces = (surfaces) ? surfaces->size() : 0;
648 std::size_t nTotalTriangles(0);
649 for (std::size_t i = 0; i < nSurfaces; ++i)
651 nTotalTriangles += (*surfaces)[i]->getNumberOfTriangles();
653 out << nTotalTriangles <<
" 1\n";
655 for (std::size_t i = 0; i < nSurfaces; ++i)
657 const std::size_t nTriangles((*surfaces)[i]->getNumberOfTriangles());
658 const std::size_t marker(i + 1);
660 for (std::size_t j = 0; j < nTriangles; ++j)
663 out <<
"3 " << tri[0] <<
" " << tri[1] <<
" " << tri[2] <<
" "
669 if (attribute_points.empty())
675 const std::size_t nAttributePoints(attribute_points.size());
676 out << nAttributePoints <<
"\n";
677 for (std::size_t i = 0; i < nAttributePoints; ++i)
679 out << i + 1 <<
" " << attribute_points[i][0] <<
" "
680 << attribute_points[i][1] <<
" " << attribute_points[i][2]
681 <<
" " << 10 * attribute_points[i].getID() <<
"\n";
685 "TetGenInterface::writeTetGenSmesh() - {:d} points and {:d} surfaces "
686 "successfully written.",
694 const std::string& file_name,
696 std::vector<MeshLib::Node>& attribute_points)
const
703 const std::vector<MeshLib::Node*>& nodes = mesh.
getNodes();
705 std::ofstream out(file_name.c_str(), std::ios::out);
706 out.precision(std::numeric_limits<double>::digits10);
708 const std::size_t nPoints(nodes.size());
709 out << nPoints <<
" 3\n";
711 for (std::size_t i = 0; i < nPoints; ++i)
713 out << i <<
" " << (*nodes[i])[0] <<
" " << (*nodes[i])[1] <<
" "
714 << (*nodes[i])[2] <<
"\n";
729 if (attribute_points.empty())
735 const std::size_t nAttributePoints(attribute_points.size());
736 out << nAttributePoints <<
"\n";
737 for (std::size_t i = 0; i < nAttributePoints; ++i)
739 out << i + 1 <<
" " << attribute_points[i][0] <<
" "
740 << attribute_points[i][1] <<
" " << attribute_points[i][2]
741 <<
" " << 10 * attribute_points[i].getID() <<
"\n";
746 "TetGenInterface::writeTetGenPoly() - {:d} points and {:d} surfaces "
747 "successfully written.",
759 std::size_t
const n_tri =
763 std::size_t
const n_quad =
767 unsigned const nTotalTriangles = n_tri + 2 * n_quad;
768 out << nTotalTriangles <<
" 1\n";
770 const std::vector<MeshLib::Element*>& elements = mesh.
getElements();
775 const std::size_t nElements(elements.size());
776 unsigned element_count(0);
777 for (std::size_t i = 0; i < nElements; ++i)
779 std::string
const matId = mat_ids ? std::to_string((*mat_ids)[i]) :
"";
787 std::vector<MeshLib::Node>& attribute_points)
const
789 const std::vector<MeshLib::Element*>& elements = mesh.
getElements();
790 const std::size_t nElements(elements.size());
791 if (!attribute_points.empty())
793 attribute_points.clear();
798 const std::streamoff before_elems_pos(out.tellp());
799 const unsigned n_spaces(
800 static_cast<unsigned>(std::floor(log(nElements * 8))) + 1);
801 out << std::string(n_spaces,
' ') <<
" 1\n";
802 auto const*
const materialIds =
804 unsigned element_count(0);
805 for (std::size_t i = 0; i < nElements; ++i)
812 const unsigned nFaces(elements[i]->getNumberOfNeighbors());
813 std::string
const mat_id_str =
814 materialIds ? std::to_string((*materialIds)[i]) :
"";
815 for (std::size_t j = 0; j < nFaces; ++j)
819 if (neighbor && materialIds &&
820 (*materialIds)[i] <= (*materialIds)[neighbor->
getID()])
825 std::unique_ptr<MeshLib::Element const>
const face(
826 elements[i]->getFace(j));
831 attribute_points.emplace_back(
837 const std::streamoff after_elems_pos(out.tellp());
838 out.seekp(before_elems_pos);
839 out << element_count;
840 out.seekp(after_elems_pos);
845 unsigned& element_count,
846 std::string
const& matId)
853 <<
" " << matId <<
" # " << element_count <<
"\n";
859 <<
" " << matId <<
" # " << element_count <<
"\n";
863 <<
" " << matId <<
" # " << element_count <<
"\n";
Definition of the Element class.
void INFO(char const *fmt, Args const &... args)
void ERR(char const *fmt, Args const &... args)
void WARN(char const *fmt, Args const &... args)
Definition of the Mesh class.
Definition of the Node class.
Definition of the TetGenInterface class.
Definition of the Tet class.
bool readElementsFromStream(std::ifstream &ins, std::vector< MeshLib::Element * > &elements, std::vector< int > &materials, const std::vector< MeshLib::Node * > &nodes) const
bool _boundary_markers
true if boundary markers are set, false otherwise
MeshLib::Mesh * readTetGenMesh(std::string const &nodes_fname, std::string const &ele_fname)
bool parseNodes(std::ifstream &ins, std::vector< MeshLib::Node * > &nodes, std::size_t n_nodes, std::size_t dim)
static bool parseElementsFileHeader(std::string &line, std::size_t &n_tets, std::size_t &n_nodes_per_tet, bool ®ion_attribute)
void write3dElements(std::ofstream &out, const MeshLib::Mesh &mesh, std::vector< MeshLib::Node > &attribute_points) const
bool readNodesFromStream(std::ifstream &ins, std::vector< MeshLib::Node * > &nodes)
static bool parseNodesFileHeader(std::string const &line, std::size_t &n_nodes, std::size_t &dim, std::size_t &n_attributes, bool &boundary_markers)
bool _zero_based_idx
the value is true if the indexing is zero based, else false
bool readTetGenGeometry(std::string const &geo_fname, GeoLib::GEOObjects &geo_objects)
void write2dElements(std::ofstream &out, const MeshLib::Mesh &mesh) const
bool parseElements(std::ifstream &ins, std::vector< MeshLib::Element * > &elements, std::vector< int > &materials, const std::vector< MeshLib::Node * > &nodes, std::size_t n_tets, std::size_t n_nodes_per_tet, bool region_attribute) const
static bool writeTetGenSmesh(const std::string &file_name, const GeoLib::GEOObjects &geo_objects, const std::string &geo_name, const std::vector< GeoLib::Point > &attribute_points)
static void writeElementToFacets(std::ofstream &out, const MeshLib::Element &element, unsigned &element_count, std::string const &matId)
Writes facet information from a 2D element to the stream and increments the total element count accor...
bool parseSmeshFacets(std::ifstream &input, std::vector< GeoLib::Surface * > &surfaces, const std::vector< GeoLib::Point * > &points, const std::vector< std::size_t > &pnt_id_map)
std::size_t getNFacets(std::ifstream &input)
Returns the declared number of facets in the poly file.
Container class for geometric objects.
void addSurfaceVec(std::unique_ptr< std::vector< Surface * >> sfc, const std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> sfc_names=nullptr)
const std::vector< Point * > * getPointVec(const std::string &name) const
const PointVec * getPointVecObj(const std::string &name) const
void addPointVec(std::unique_ptr< std::vector< Point * >> points, std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> pnt_id_name_map=nullptr, double eps=std::sqrt(std::numeric_limits< double >::epsilon()))
const std::vector< Surface * > * getSurfaceVec(const std::string &name) const
Returns the surface vector with the given name as a const.
const std::vector< std::size_t > & getIDMap() const
A Surface is represented by Triangles. It consists of a reference to a vector of (pointers to) points...
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
virtual MeshElemType getGeomType() const =0
virtual std::size_t getID() const final
Returns the ID of the element.
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Properties & getProperties()
std::size_t getNumberOfElements() const
Get the number of elements.
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
PropertyVector< T > const * getPropertyVector(std::string const &name) const
bool existsPropertyVector(std::string const &name) const
PropertyVector< T > * createNewPropertyVector(std::string const &name, MeshItemType mesh_item_type, std::size_t n_components=1)
void simplify(std::string &str)
std::string getFileExtension(const std::string &path)
std::string extractBaseNameWithoutExtension(std::string const &pathname)
std::vector< std::string > splitString(std::string const &str)
void cleanupVectorElements(std::vector< T1 * > const &items, std::vector< T2 * > const &dependent_items)
void copy(PETScVector const &x, PETScVector &y)
std::size_t getNodeIndex(Element const &element, unsigned const idx)
MeshLib::Node getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.
unsigned getDimension(MeshLib::MeshElemType eleType)