19#include <range/v3/numeric/accumulate.hpp>
38 std::vector<GeoLib::Point*> points;
39 points.reserve(nodes.size());
40 std::transform(nodes.begin(), nodes.end(), std::back_inserter(points),
42 { return new GeoLib::Point(*node, node->getID()); });
49 std::ifstream poly_stream(geo_fname.c_str());
53 ERR(
"TetGenInterface::readTetGenGeometry() failed to open {:s}",
60 ERR(
"TetGenInterface::readTetGenGeometry() - unknown file type (only "
61 "*.smesh is supported).");
65 std::vector<MeshLib::Node*> nodes;
76 geo_objects.
addPointVec(std::move(points), geo_name);
77 const std::vector<std::size_t>& id_map(
80 std::vector<GeoLib::Surface*> surfaces;
98 std::getline(input, line);
101 ERR(
"TetGenInterface::getNFacets(): Error reading number of "
107 if (line.empty() || line.compare(0, 1,
"#") == 0)
113 auto it = fields.begin();
115 if (fields.size() > 1)
125 std::ifstream& input,
126 std::vector<GeoLib::Surface*>& surfaces,
127 const std::vector<GeoLib::Point*>& points,
128 const std::vector<std::size_t>& pnt_id_map)
130 const std::size_t nFacets(this->
getNFacets(input));
132 surfaces.reserve(nFacets);
133 std::list<std::string>::const_iterator it;
136 std::vector<std::size_t> idx_map;
139 while (k < nFacets && !input.fail())
141 std::getline(input, line);
144 ERR(
"TetGenInterface::parseFacets(): Error reading facet {:d}.", k);
149 if (line.empty() || line.compare(0, 1,
"#") == 0)
155 const std::list<std::string> point_fields =
157 it = point_fields.begin();
161 ERR(
"Smesh-files are currently only supported for triangle "
165 const std::size_t point_field_size =
167 if (point_fields.size() > point_field_size)
169 std::vector<std::size_t> point_ids;
170 for (std::size_t j(0); j < nPoints; ++j)
177 const std::size_t sfc_marker =
180 const std::size_t idx =
181 std::find(idx_map.begin(), idx_map.end(), sfc_marker) -
183 if (idx >= surfaces.size())
185 idx_map.push_back(sfc_marker);
188 surfaces[idx]->addTriangle(point_ids[0], point_ids[1],
193 ERR(
"TetGenInterface::parseFacets(): Error reading points for "
205 auto const nTotalTriangles = ranges::accumulate(
206 surfaces, std::size_t{0}, {},
207 [](
auto const* surface) {
return surface->getNumberOfTriangles(); });
208 if (nTotalTriangles == nFacets)
213 ERR(
"TetGenInterface::parseFacets(): Number of expected total triangles "
214 "({:d}) does not match number of found triangles ({:d}).",
221 std::string
const& ele_fname)
223 std::ifstream ins_nodes(nodes_fname.c_str());
224 std::ifstream ins_ele(ele_fname.c_str());
226 if (!ins_nodes || !ins_ele)
230 ERR(
"TetGenInterface::readTetGenMesh failed to open {:s}",
235 ERR(
"TetGenInterface::readTetGenMesh failed to open {:s}",
241 std::vector<MeshLib::Node*> nodes;
249 std::vector<MeshLib::Element*> elements;
250 std::vector<int> materials;
259 if (std::any_of(materials.cbegin(), materials.cend(),
260 [](
int m) { return m != 0; }))
264 mat_props->assign(materials);
267 const std::string mesh_name(
274 std::vector<MeshLib::Node*>& nodes)
277 std::getline(ins, line);
280 std::size_t n_attributes;
281 bool boundary_markers;
286 if (line.empty() || line.compare(0, 1,
"#") == 0)
289 std::getline(ins, line);
294 n_attributes, boundary_markers);
309 std::size_t& n_nodes,
311 std::size_t& n_attributes,
312 bool& boundary_markers)
315 if (pnt_header.empty())
317 ERR(
"TetGenInterface::parseNodesFileHeader(): could not read number of "
318 "nodes specified in header.");
321 auto it = pnt_header.begin();
323 dim = (pnt_header.size() == 1) ? 3
326 if (pnt_header.size() < 4)
329 boundary_markers =
false;
334 boundary_markers = *(++it) ==
"1";
340 std::vector<MeshLib::Node*>& nodes,
345 nodes.reserve(n_nodes);
348 while (k < n_nodes && !ins.fail())
350 std::vector<double> coordinates(dim);
351 std::getline(ins, line);
354 ERR(
"TetGenInterface::parseNodes(): Error reading node {:d}.", k);
359 std::size_t pos_end = 0;
360 std::size_t pos_beg = line.find_first_not_of(
' ', pos_end);
361 pos_end = line.find_first_of(
" \n", pos_beg);
363 if (line.empty() || pos_beg == pos_end ||
364 line.compare(pos_beg, 1,
"#") == 0)
369 if (pos_beg != std::string::npos && pos_end != std::string::npos)
372 line.substr(pos_beg, pos_end - pos_beg));
373 if (k == 0 &&
id == 0)
380 ERR(
"TetGenInterface::parseNodes(): Error reading ID of node {:d}.",
386 for (std::size_t i(0); i < dim; i++)
388 pos_beg = line.find_first_not_of(
' ', pos_end);
389 pos_end = line.find_first_of(
" \n", pos_beg);
390 if (pos_end == std::string::npos)
392 pos_end = line.size();
394 if (pos_beg != std::string::npos)
397 line.substr(pos_beg, pos_end - pos_beg));
401 ERR(
"TetGenInterface::parseNodes(): error reading coordinate "
402 "{:d} of node {:d}.",
409 nodes.push_back(
new MeshLib::Node(coordinates.data(),
id - offset));
420 std::vector<MeshLib::Element*>& elements,
421 std::vector<int>& materials,
422 const std::vector<MeshLib::Node*>& nodes)
const
425 std::getline(ins, line);
427 std::size_t n_nodes_per_tet;
428 bool region_attributes;
433 if (line.empty() || line.compare(0, 1,
"#") == 0)
436 std::getline(ins, line);
442 line, n_tets, n_nodes_per_tet, region_attributes);
464 std::size_t& n_nodes_per_tet,
465 bool& region_attribute)
471 pos_beg = line.find_first_not_of(
' ');
472 pos_end = line.find_first_of(
' ', pos_beg);
473 if (pos_beg != std::string::npos && pos_end != std::string::npos)
476 line.substr(pos_beg, pos_end - pos_beg));
480 ERR(
"TetGenInterface::parseElementsFileHeader(): Could not read number "
481 "of tetrahedra specified in header.");
485 pos_beg = line.find_first_not_of(
" \t", pos_end);
486 pos_end = line.find_first_of(
" \t", pos_beg);
488 line.substr(pos_beg, pos_end - pos_beg));
490 pos_beg = line.find_first_not_of(
" \t", pos_end);
491 pos_end = line.find_first_of(
" \t\n", pos_beg);
492 if (pos_end == std::string::npos)
494 pos_end = line.size();
496 region_attribute = line.substr(pos_beg, pos_end - pos_beg) ==
"1";
502 std::vector<MeshLib::Element*>& elements,
503 std::vector<int>& materials,
504 const std::vector<MeshLib::Node*>& nodes,
506 std::size_t n_nodes_per_tet,
507 bool region_attribute)
const
510 std::vector<std::size_t> ids(n_nodes_per_tet);
512 elements.reserve(n_tets);
513 materials.reserve(n_tets);
516 for (std::size_t k(0); k < n_tets && !ins.fail(); k++)
518 std::getline(ins, line);
521 ERR(
"TetGenInterface::parseElements(): Error reading tetrahedron "
527 std::size_t pos_end = 0;
528 std::size_t pos_beg = line.find_first_not_of(
' ', pos_end);
529 pos_end = line.find_first_of(
" \n", pos_beg);
531 if (line.empty() || pos_beg == pos_end ||
532 line.compare(pos_beg, 1,
"#") == 0)
538 if (pos_beg == std::string::npos || pos_end == std::string::npos)
540 ERR(
"TetGenInterface::parseElements(): Error reading id of "
547 for (std::size_t i(0); i < n_nodes_per_tet; i++)
549 pos_beg = line.find_first_not_of(
' ', pos_end);
550 pos_end = line.find_first_of(
' ', pos_beg);
551 if (pos_end == std::string::npos)
553 pos_end = line.size();
555 if (pos_beg != std::string::npos && pos_end != std::string::npos)
558 line.substr(pos_beg, pos_end - pos_beg)) -
563 ERR(
"TetGenInterface::parseElements(): Error reading node {:d} "
564 "of tetrahedron {:d}.",
573 if (region_attribute)
575 pos_beg = line.find_first_not_of(
' ', pos_end);
576 pos_end = line.find_first_of(
' ', pos_beg);
577 if (pos_end == std::string::npos)
579 pos_end = line.size();
581 if (pos_beg != std::string::npos && pos_end != std::string::npos)
584 line.substr(pos_beg, pos_end - pos_beg));
588 ERR(
"TetGenInterface::parseElements(): Error reading region "
589 "attribute of tetrahedron {:d}.",
596 for (
int n = 0; n < 4; n++)
598 tet_nodes[n] = nodes[ids[n]];
601 materials.push_back(region);
608 const std::string& file_name,
610 const std::string& geo_name,
611 const std::vector<GeoLib::Point>& attribute_points)
613 std::vector<GeoLib::Point*>
const*
const points =
615 std::vector<GeoLib::Surface*>
const*
const surfaces =
618 if (points ==
nullptr)
620 ERR(
"Geometry {:s} not found.", geo_name);
623 if (surfaces ==
nullptr)
625 WARN(
"No surfaces found for geometry {:s}. Writing points only.",
629 std::ofstream out(file_name.c_str(), std::ios::out);
630 out.precision(std::numeric_limits<double>::max_digits10);
632 const std::size_t nPoints(points->size());
633 out << nPoints <<
" 3\n";
635 for (std::size_t i = 0; i < nPoints; ++i)
637 out << i <<
" " << (*(*points)[i])[0] <<
" " << (*(*points)[i])[1]
638 <<
" " << (*(*points)[i])[2] <<
"\n";
641 const std::size_t nSurfaces = (surfaces) ? surfaces->size() : 0;
642 std::size_t nTotalTriangles(0);
643 for (std::size_t i = 0; i < nSurfaces; ++i)
645 nTotalTriangles += (*surfaces)[i]->getNumberOfTriangles();
647 out << nTotalTriangles <<
" 1\n";
649 for (std::size_t i = 0; i < nSurfaces; ++i)
651 const std::size_t nTriangles((*surfaces)[i]->getNumberOfTriangles());
652 const std::size_t marker(i + 1);
654 for (std::size_t j = 0; j < nTriangles; ++j)
657 out <<
"3 " << tri[0] <<
" " << tri[1] <<
" " << tri[2] <<
" "
663 if (attribute_points.empty())
669 const std::size_t nAttributePoints(attribute_points.size());
670 out << nAttributePoints <<
"\n";
671 for (std::size_t i = 0; i < nAttributePoints; ++i)
673 out << i + 1 <<
" " << attribute_points[i][0] <<
" "
674 << attribute_points[i][1] <<
" " << attribute_points[i][2]
675 <<
" " << 10 * attribute_points[i].getID() <<
"\n";
679 "TetGenInterface::writeTetGenSmesh() - {:d} points and {:d} surfaces "
680 "successfully written.",
688 const std::string& file_name,
690 std::vector<MeshLib::Node>& attribute_points)
const
697 const std::vector<MeshLib::Node*>& nodes = mesh.
getNodes();
699 std::ofstream out(file_name.c_str(), std::ios::out);
700 out.precision(std::numeric_limits<double>::max_digits10);
702 const std::size_t nPoints(nodes.size());
703 out << nPoints <<
" 3\n";
705 for (std::size_t i = 0; i < nPoints; ++i)
707 out << i <<
" " << (*nodes[i])[0] <<
" " << (*nodes[i])[1] <<
" "
708 << (*nodes[i])[2] <<
"\n";
723 if (attribute_points.empty())
729 const std::size_t nAttributePoints(attribute_points.size());
730 out << nAttributePoints <<
"\n";
731 for (std::size_t i = 0; i < nAttributePoints; ++i)
733 out << i + 1 <<
" " << attribute_points[i][0] <<
" "
734 << attribute_points[i][1] <<
" " << attribute_points[i][2]
735 <<
" " << 10 * attribute_points[i].getID() <<
"\n";
740 "TetGenInterface::writeTetGenPoly() - {:d} points and {:d} surfaces "
741 "successfully written.",
754 std::size_t
const n_tri =
758 std::size_t
const n_quad =
762 unsigned const nTotalTriangles = n_tri + 2 * n_quad;
763 out << nTotalTriangles <<
" 1\n";
765 const std::vector<MeshLib::Element*>& elements = mesh.
getElements();
770 const std::size_t nElements(elements.size());
771 unsigned element_count(0);
772 for (std::size_t i = 0; i < nElements; ++i)
774 std::string
const matId = mat_ids ? std::to_string((*mat_ids)[i]) :
"";
782 std::vector<MeshLib::Node>& attribute_points)
const
784 const std::vector<MeshLib::Element*>& elements = mesh.
getElements();
785 const std::size_t nElements(elements.size());
786 if (!attribute_points.empty())
788 attribute_points.clear();
793 const std::streamoff before_elems_pos(out.tellp());
794 const unsigned n_spaces(
795 static_cast<unsigned>(std::floor(log(nElements * 8))) + 1);
796 out << std::string(n_spaces,
' ') <<
" 1\n";
797 auto const*
const materialIds =
799 unsigned element_count(0);
800 for (std::size_t i = 0; i < nElements; ++i)
802 if (elements[i]->getDimension() < 3)
807 const unsigned nFaces(elements[i]->getNumberOfNeighbors());
808 std::string
const mat_id_str =
809 materialIds ? std::to_string((*materialIds)[i]) :
"";
810 for (std::size_t j = 0; j < nFaces; ++j)
814 if (neighbor && materialIds &&
815 (*materialIds)[i] <= (*materialIds)[neighbor->
getID()])
820 std::unique_ptr<MeshLib::Element const>
const face(
821 elements[i]->getFace(j));
826 attribute_points.emplace_back(
832 const std::streamoff after_elems_pos(out.tellp());
833 out.seekp(before_elems_pos);
834 out << element_count;
835 out.seekp(after_elems_pos);
840 unsigned& element_count,
841 std::string
const& matId)
846 out <<
"3 " << getNodeIndex(element, 0) <<
" "
847 << getNodeIndex(element, 1) <<
" " << getNodeIndex(element, 2)
848 <<
" " << matId <<
" # " << element_count <<
"\n";
852 out <<
"3 " << getNodeIndex(element, 0) <<
" "
853 << getNodeIndex(element, 1) <<
" " << getNodeIndex(element, 2)
854 <<
" " << matId <<
" # " << element_count <<
"\n";
856 out <<
"3 " << getNodeIndex(element, 0) <<
" "
857 << getNodeIndex(element, 2) <<
" " << getNodeIndex(element, 3)
858 <<
" " << matId <<
" # " << element_count <<
"\n";
Definition of the Element class.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... 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)
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.
const std::vector< Point * > * getPointVec(const std::string &name) const
void addPointVec(std::vector< Point * > &&points, std::string &name, PointVec::NameIdMap &&pnt_id_name_map, double const eps=std::sqrt(std::numeric_limits< double >::epsilon()))
const PointVec * getPointVecObj(const std::string &name) const
const std::vector< Surface * > * getSurfaceVec(const std::string &name) const
Returns the surface vector with the given name as a const.
void addSurfaceVec(std::vector< Surface * > &&sfc, const std::string &name, SurfaceVec::NameIdMap &&sfc_names)
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...
std::map< std::string, std::size_t > NameIdMap
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
virtual MeshElemType getGeomType() const =0
std::size_t getID() const
Returns the ID of the element.
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
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....
bool existsPropertyVector(std::string_view name) const
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
PropertyVector< T > const * getPropertyVector(std::string_view name) const
void simplify(std::string &str)
std::string getFileExtension(const std::string &path)
void cleanupVectorElements(std::vector< T * > &items)
std::string extractBaseNameWithoutExtension(std::string const &pathname)
std::vector< std::string > splitString(std::string const &str)
T str2number(const std::string &str)
auto constructPointsFromNodes(std::vector< MeshLib::Node * > nodes)
MathLib::Point3d getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.