OGS
TetGenInterface.cpp
Go to the documentation of this file.
1
15#include "TetGenInterface.h"
16
17#include <cstddef>
18#include <fstream>
19#include <range/v3/numeric/accumulate.hpp>
20#include <string>
21
22#include "BaseLib/FileTools.h"
23#include "BaseLib/Logging.h"
24#include "BaseLib/StringTools.h"
25#include "GeoLib/Triangle.h"
28#include "MeshLib/Mesh.h"
29#include "MeshLib/Node.h"
31
32namespace FileIO
33{
35
36auto constructPointsFromNodes(std::vector<MeshLib::Node*> nodes)
37{
38 std::vector<GeoLib::Point*> points;
39 points.reserve(nodes.size());
40 std::transform(nodes.begin(), nodes.end(), std::back_inserter(points),
41 [](auto const* node)
42 { return new GeoLib::Point(*node, node->getID()); });
43 return points;
44}
45
46bool TetGenInterface::readTetGenGeometry(std::string const& geo_fname,
47 GeoLib::GEOObjects& geo_objects)
48{
49 std::ifstream poly_stream(geo_fname.c_str());
50
51 if (!poly_stream)
52 {
53 ERR("TetGenInterface::readTetGenGeometry() failed to open {:s}",
54 geo_fname);
55 return false;
56 }
57 std::string ext(BaseLib::getFileExtension(geo_fname));
58 if (ext != ".smesh")
59 {
60 ERR("TetGenInterface::readTetGenGeometry() - unknown file type (only "
61 "*.smesh is supported).");
62 return false;
63 }
64
65 std::vector<MeshLib::Node*> nodes;
66 if (!readNodesFromStream(poly_stream, nodes))
67 {
68 // remove nodes read until now
70 return false;
71 }
72 auto points = constructPointsFromNodes(nodes);
74
75 std::string geo_name(BaseLib::extractBaseNameWithoutExtension(geo_fname));
76 geo_objects.addPointVec(std::move(points), geo_name);
77 const std::vector<std::size_t>& id_map(
78 geo_objects.getPointVecObj(geo_name)->getIDMap());
79
80 std::vector<GeoLib::Surface*> surfaces;
81 if (!parseSmeshFacets(poly_stream, surfaces,
82 *geo_objects.getPointVec(geo_name), id_map))
83 {
84 // remove surfaces read until now but keep the points
86 }
87 geo_objects.addSurfaceVec(std::move(surfaces), geo_name,
89
90 return true;
91}
92
93std::size_t TetGenInterface::getNFacets(std::ifstream& input)
94{
95 std::string line;
96 while (!input.fail())
97 {
98 std::getline(input, line);
99 if (input.fail())
100 {
101 ERR("TetGenInterface::getNFacets(): Error reading number of "
102 "facets.");
103 return 0;
104 }
105
106 BaseLib::simplify(line);
107 if (line.empty() || line.compare(0, 1, "#") == 0)
108 {
109 continue;
110 }
111
112 const std::list<std::string> fields = BaseLib::splitString(line, ' ');
113 auto it = fields.begin();
114 const auto nFacets(BaseLib::str2number<std::size_t>(*it));
115 if (fields.size() > 1)
116 {
117 _boundary_markers = BaseLib::str2number<std::size_t>(*(++it)) != 0;
118 }
119 return nFacets;
120 }
121 return 0;
122}
123
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)
129{
130 const std::size_t nFacets(this->getNFacets(input));
131 std::string line;
132 surfaces.reserve(nFacets);
133 std::list<std::string>::const_iterator it;
134
135 const unsigned offset = (_zero_based_idx) ? 0 : 1;
136 std::vector<std::size_t> idx_map;
137
138 std::size_t k(0);
139 while (k < nFacets && !input.fail())
140 {
141 std::getline(input, line);
142 if (input.fail())
143 {
144 ERR("TetGenInterface::parseFacets(): Error reading facet {:d}.", k);
145 return false;
146 }
147
148 BaseLib::simplify(line);
149 if (line.empty() || line.compare(0, 1, "#") == 0)
150 {
151 continue;
152 }
153
154 // read facets
155 const std::list<std::string> point_fields =
156 BaseLib::splitString(line, ' ');
157 it = point_fields.begin();
158 const auto nPoints = BaseLib::str2number<std::size_t>(*it);
159 if (nPoints != 3)
160 {
161 ERR("Smesh-files are currently only supported for triangle "
162 "meshes.");
163 return false;
164 }
165 const std::size_t point_field_size =
166 (_boundary_markers) ? nPoints + 1 : nPoints;
167 if (point_fields.size() > point_field_size)
168 {
169 std::vector<std::size_t> point_ids;
170 for (std::size_t j(0); j < nPoints; ++j)
171 {
172 point_ids.push_back(
173 pnt_id_map[BaseLib::str2number<std::size_t>(*(++it)) -
174 offset]);
175 }
176
177 const std::size_t sfc_marker =
178 (_boundary_markers) ? BaseLib::str2number<std::size_t>(*(++it))
179 : 0;
180 const std::size_t idx =
181 std::find(idx_map.begin(), idx_map.end(), sfc_marker) -
182 idx_map.begin();
183 if (idx >= surfaces.size())
184 {
185 idx_map.push_back(sfc_marker);
186 surfaces.push_back(new GeoLib::Surface(points));
187 }
188 surfaces[idx]->addTriangle(point_ids[0], point_ids[1],
189 point_ids[2]);
190 }
191 else
192 {
193 ERR("TetGenInterface::parseFacets(): Error reading points for "
194 "facet {:d}.",
195 k);
196 return false;
197 }
198 ++k;
199 }
200 // here the poly-file potentially defines a number of points to mark holes
201 // within the volumes defined by the facets, these are ignored for now here
202 // the poly-file potentially defines a number of region attributes, these
203 // are ignored for now
204
205 auto const nTotalTriangles = ranges::accumulate(
206 surfaces, std::size_t{0}, {},
207 [](auto const* surface) { return surface->getNumberOfTriangles(); });
208 if (nTotalTriangles == nFacets)
209 {
210 return true;
211 }
212
213 ERR("TetGenInterface::parseFacets(): Number of expected total triangles "
214 "({:d}) does not match number of found triangles ({:d}).",
215 surfaces.size(),
216 nTotalTriangles);
217 return false;
218}
219
220MeshLib::Mesh* TetGenInterface::readTetGenMesh(std::string const& nodes_fname,
221 std::string const& ele_fname)
222{
223 std::ifstream ins_nodes(nodes_fname.c_str());
224 std::ifstream ins_ele(ele_fname.c_str());
225
226 if (!ins_nodes || !ins_ele)
227 {
228 if (!ins_nodes)
229 {
230 ERR("TetGenInterface::readTetGenMesh failed to open {:s}",
231 nodes_fname);
232 }
233 if (!ins_ele)
234 {
235 ERR("TetGenInterface::readTetGenMesh failed to open {:s}",
236 ele_fname);
237 }
238 return nullptr;
239 }
240
241 std::vector<MeshLib::Node*> nodes;
242 if (!readNodesFromStream(ins_nodes, nodes))
243 {
244 // remove nodes read until now
246 return nullptr;
247 }
248
249 std::vector<MeshLib::Element*> elements;
250 std::vector<int> materials;
251 if (!readElementsFromStream(ins_ele, elements, materials, nodes))
252 {
253 BaseLib::cleanupVectorElements(nodes, elements);
254 return nullptr;
255 }
256
257 MeshLib::Properties properties;
258 // Transmit material values if there is any material value != 0
259 if (std::any_of(materials.cbegin(), materials.cend(),
260 [](int m) { return m != 0; }))
261 {
262 auto* const mat_props = properties.createNewPropertyVector<int>(
263 "MaterialIDs", MeshLib::MeshItemType::Cell);
264 mat_props->reserve(elements.size());
265 std::copy(materials.cbegin(),
266 materials.cend(),
267 std::back_inserter(*mat_props));
268 }
269
270 const std::string mesh_name(
272 return new MeshLib::Mesh(mesh_name, nodes, elements, properties);
273}
274
276 std::vector<MeshLib::Node*>& nodes)
277{
278 std::string line;
279 std::getline(ins, line);
280 std::size_t n_nodes;
281 std::size_t dim;
282 std::size_t n_attributes;
283 bool boundary_markers;
284
285 while (!ins.fail())
286 {
287 BaseLib::simplify(line);
288 if (line.empty() || line.compare(0, 1, "#") == 0)
289 {
290 // this line is a comment - skip
291 std::getline(ins, line);
292 continue;
293 }
294 // read header line
295 bool header_okay = parseNodesFileHeader(line, n_nodes, dim,
296 n_attributes, boundary_markers);
297 if (!header_okay)
298 {
299 return false;
300 }
301 if (!parseNodes(ins, nodes, n_nodes, dim))
302 {
303 return false;
304 }
305 return true;
306 }
307 return false;
308}
309
310bool TetGenInterface::parseNodesFileHeader(std::string const& line,
311 std::size_t& n_nodes,
312 std::size_t& dim,
313 std::size_t& n_attributes,
314 bool& boundary_markers)
315{
316 std::list<std::string> pnt_header = BaseLib::splitString(line, ' ');
317 if (pnt_header.empty())
318 {
319 ERR("TetGenInterface::parseNodesFileHeader(): could not read number of "
320 "nodes specified in header.");
321 return false;
322 }
323 auto it = pnt_header.begin();
324 n_nodes = BaseLib::str2number<std::size_t>(*it);
325 dim = (pnt_header.size() == 1) ? 3
326 : BaseLib::str2number<std::size_t>(*(++it));
327
328 if (pnt_header.size() < 4)
329 {
330 n_attributes = 0;
331 boundary_markers = false;
332 return true;
333 }
334
335 n_attributes = BaseLib::str2number<std::size_t>(*(++it));
336 boundary_markers = *(++it) == "1";
337
338 return true;
339}
340
341bool TetGenInterface::parseNodes(std::ifstream& ins,
342 std::vector<MeshLib::Node*>& nodes,
343 std::size_t n_nodes,
344 std::size_t dim)
345{
346 std::string line;
347 nodes.reserve(n_nodes);
348
349 std::size_t k(0);
350 while (k < n_nodes && !ins.fail())
351 {
352 std::vector<double> coordinates(dim);
353 std::getline(ins, line);
354 if (ins.fail())
355 {
356 ERR("TetGenInterface::parseNodes(): Error reading node {:d}.", k);
357 return false;
358 }
359
360 std::size_t id;
361 std::size_t pos_end = 0;
362 std::size_t pos_beg = line.find_first_not_of(' ', pos_end);
363 pos_end = line.find_first_of(" \n", pos_beg);
364
365 if (line.empty() || pos_beg == pos_end ||
366 line.compare(pos_beg, 1, "#") == 0)
367 {
368 continue;
369 }
370
371 if (pos_beg != std::string::npos && pos_end != std::string::npos)
372 {
373 id = BaseLib::str2number<std::size_t>(
374 line.substr(pos_beg, pos_end - pos_beg));
375 if (k == 0 && id == 0)
376 {
377 _zero_based_idx = true;
378 }
379 }
380 else
381 {
382 ERR("TetGenInterface::parseNodes(): Error reading ID of node {:d}.",
383 k);
384 return false;
385 }
386 // read coordinates
387 const unsigned offset = (_zero_based_idx) ? 0 : 1;
388 for (std::size_t i(0); i < dim; i++)
389 {
390 pos_beg = line.find_first_not_of(' ', pos_end);
391 pos_end = line.find_first_of(" \n", pos_beg);
392 if (pos_end == std::string::npos)
393 {
394 pos_end = line.size();
395 }
396 if (pos_beg != std::string::npos)
397 {
398 coordinates[i] = BaseLib::str2number<double>(
399 line.substr(pos_beg, pos_end - pos_beg));
400 }
401 else
402 {
403 ERR("TetGenInterface::parseNodes(): error reading coordinate "
404 "{:d} of node {:d}.",
405 i,
406 k);
407 return false;
408 }
409 }
410
411 nodes.push_back(new MeshLib::Node(coordinates.data(), id - offset));
412 // read attributes and boundary markers ... - at the moment we do not
413 // use this information
414 ++k;
415 }
416
417 return true;
418}
419
421 std::ifstream& ins,
422 std::vector<MeshLib::Element*>& elements,
423 std::vector<int>& materials,
424 const std::vector<MeshLib::Node*>& nodes) const
425{
426 std::string line;
427 std::getline(ins, line);
428 std::size_t n_tets;
429 std::size_t n_nodes_per_tet;
430 bool region_attributes;
431
432 while (!ins.fail())
433 {
434 BaseLib::simplify(line);
435 if (line.empty() || line.compare(0, 1, "#") == 0)
436 {
437 // this line is a comment - skip
438 std::getline(ins, line);
439 continue;
440 }
441
442 // read header line
443 bool header_okay = parseElementsFileHeader(
444 line, n_tets, n_nodes_per_tet, region_attributes);
445 if (!header_okay)
446 {
447 return false;
448 }
449 if (!parseElements(ins,
450 elements,
451 materials,
452 nodes,
453 n_tets,
454 n_nodes_per_tet,
455 region_attributes))
456 {
457 return false;
458 }
459 return true;
460 }
461 return false;
462}
463
465 std::size_t& n_tets,
466 std::size_t& n_nodes_per_tet,
467 bool& region_attribute)
468{
469 std::size_t pos_beg;
470 std::size_t pos_end;
471
472 // number of tetrahedras
473 pos_beg = line.find_first_not_of(' ');
474 pos_end = line.find_first_of(' ', pos_beg);
475 if (pos_beg != std::string::npos && pos_end != std::string::npos)
476 {
477 n_tets = BaseLib::str2number<std::size_t>(
478 line.substr(pos_beg, pos_end - pos_beg));
479 }
480 else
481 {
482 ERR("TetGenInterface::parseElementsFileHeader(): Could not read number "
483 "of tetrahedra specified in header.");
484 return false;
485 }
486 // nodes per tet - either 4 or 10
487 pos_beg = line.find_first_not_of(" \t", pos_end);
488 pos_end = line.find_first_of(" \t", pos_beg);
489 n_nodes_per_tet = BaseLib::str2number<std::size_t>(
490 line.substr(pos_beg, pos_end - pos_beg));
491 // region attribute at tetrahedra?
492 pos_beg = line.find_first_not_of(" \t", pos_end);
493 pos_end = line.find_first_of(" \t\n", pos_beg);
494 if (pos_end == std::string::npos)
495 {
496 pos_end = line.size();
497 }
498 region_attribute = line.substr(pos_beg, pos_end - pos_beg) == "1";
499
500 return true;
501}
502
503bool TetGenInterface::parseElements(std::ifstream& ins,
504 std::vector<MeshLib::Element*>& elements,
505 std::vector<int>& materials,
506 const std::vector<MeshLib::Node*>& nodes,
507 std::size_t n_tets,
508 std::size_t n_nodes_per_tet,
509 bool region_attribute) const
510{
511 std::string line;
512 std::vector<std::size_t> ids(n_nodes_per_tet);
513
514 elements.reserve(n_tets);
515 materials.reserve(n_tets);
516
517 const unsigned offset = (_zero_based_idx) ? 0 : 1;
518 for (std::size_t k(0); k < n_tets && !ins.fail(); k++)
519 {
520 std::getline(ins, line);
521 if (ins.fail())
522 {
523 ERR("TetGenInterface::parseElements(): Error reading tetrahedron "
524 "{:d}.",
525 k);
526 return false;
527 }
528
529 std::size_t pos_end = 0;
530 std::size_t pos_beg = line.find_first_not_of(' ', pos_end);
531 pos_end = line.find_first_of(" \n", pos_beg);
532
533 if (line.empty() || pos_beg == pos_end ||
534 line.compare(pos_beg, 1, "#") == 0)
535 {
536 k--;
537 continue;
538 }
539
540 if (pos_beg == std::string::npos || pos_end == std::string::npos)
541 {
542 ERR("TetGenInterface::parseElements(): Error reading id of "
543 "tetrahedron {:d}.",
544 k);
545 return false;
546 }
547
548 // read node ids
549 for (std::size_t i(0); i < n_nodes_per_tet; i++)
550 {
551 pos_beg = line.find_first_not_of(' ', pos_end);
552 pos_end = line.find_first_of(' ', pos_beg);
553 if (pos_end == std::string::npos)
554 {
555 pos_end = line.size();
556 }
557 if (pos_beg != std::string::npos && pos_end != std::string::npos)
558 {
559 ids[i] = BaseLib::str2number<std::size_t>(
560 line.substr(pos_beg, pos_end - pos_beg)) -
561 offset;
562 }
563 else
564 {
565 ERR("TetGenInterface::parseElements(): Error reading node {:d} "
566 "of tetrahedron {:d}.",
567 i,
568 k);
569 return false;
570 }
571 }
572
573 // read region attribute - this is something like material group
574 int region(0);
575 if (region_attribute)
576 {
577 pos_beg = line.find_first_not_of(' ', pos_end);
578 pos_end = line.find_first_of(' ', pos_beg);
579 if (pos_end == std::string::npos)
580 {
581 pos_end = line.size();
582 }
583 if (pos_beg != std::string::npos && pos_end != std::string::npos)
584 {
585 region = BaseLib::str2number<int>(
586 line.substr(pos_beg, pos_end - pos_beg));
587 }
588 else
589 {
590 ERR("TetGenInterface::parseElements(): Error reading region "
591 "attribute of tetrahedron {:d}.",
592 k);
593 return false;
594 }
595 }
596 // insert new element into vector
597 auto** tet_nodes = new MeshLib::Node*[4];
598 for (int n = 0; n < 4; n++)
599 {
600 tet_nodes[n] = nodes[ids[n]];
601 }
602 elements.push_back(new MeshLib::Tet(tet_nodes));
603 materials.push_back(region);
604 }
605
606 return true;
607}
608
610 const std::string& file_name,
611 const GeoLib::GEOObjects& geo_objects,
612 const std::string& geo_name,
613 const std::vector<GeoLib::Point>& attribute_points)
614{
615 std::vector<GeoLib::Point*> const* const points =
616 geo_objects.getPointVec(geo_name);
617 std::vector<GeoLib::Surface*> const* const surfaces =
618 geo_objects.getSurfaceVec(geo_name);
619
620 if (points == nullptr)
621 {
622 ERR("Geometry {:s} not found.", geo_name);
623 return false;
624 }
625 if (surfaces == nullptr)
626 {
627 WARN("No surfaces found for geometry {:s}. Writing points only.",
628 geo_name);
629 }
630
631 std::ofstream out(file_name.c_str(), std::ios::out);
632 out.precision(std::numeric_limits<double>::digits10);
633 // the points header
634 const std::size_t nPoints(points->size());
635 out << nPoints << " 3\n";
636 // the point list
637 for (std::size_t i = 0; i < nPoints; ++i)
638 {
639 out << i << " " << (*(*points)[i])[0] << " " << (*(*points)[i])[1]
640 << " " << (*(*points)[i])[2] << "\n";
641 }
642 // the surfaces header
643 const std::size_t nSurfaces = (surfaces) ? surfaces->size() : 0;
644 std::size_t nTotalTriangles(0);
645 for (std::size_t i = 0; i < nSurfaces; ++i)
646 {
647 nTotalTriangles += (*surfaces)[i]->getNumberOfTriangles();
648 }
649 out << nTotalTriangles << " 1\n";
650
651 for (std::size_t i = 0; i < nSurfaces; ++i)
652 {
653 const std::size_t nTriangles((*surfaces)[i]->getNumberOfTriangles());
654 const std::size_t marker(i + 1); // must NOT be 0!
655 // the poly list
656 for (std::size_t j = 0; j < nTriangles; ++j)
657 {
658 const GeoLib::Triangle& tri = *(*(*surfaces)[i])[j];
659 out << "3 " << tri[0] << " " << tri[1] << " " << tri[2] << " "
660 << marker << "\n";
661 }
662 }
663 out << "0\n"; // the polygon holes list
664 // the region attributes list
665 if (attribute_points.empty())
666 {
667 out << "0\n";
668 }
669 else
670 {
671 const std::size_t nAttributePoints(attribute_points.size());
672 out << nAttributePoints << "\n";
673 for (std::size_t i = 0; i < nAttributePoints; ++i)
674 {
675 out << i + 1 << " " << attribute_points[i][0] << " "
676 << attribute_points[i][1] << " " << attribute_points[i][2]
677 << " " << 10 * attribute_points[i].getID() << "\n";
678 }
679 }
680 INFO(
681 "TetGenInterface::writeTetGenSmesh() - {:d} points and {:d} surfaces "
682 "successfully written.",
683 nPoints,
684 nSurfaces);
685 out.close();
686 return true;
687}
688
690 const std::string& file_name,
691 const MeshLib::Mesh& mesh,
692 std::vector<MeshLib::Node>& attribute_points) const
693{
694 if (mesh.getDimension() == 1)
695 {
696 return false;
697 }
698
699 const std::vector<MeshLib::Node*>& nodes = mesh.getNodes();
700
701 std::ofstream out(file_name.c_str(), std::ios::out);
702 out.precision(std::numeric_limits<double>::digits10);
703 // the points header
704 const std::size_t nPoints(nodes.size());
705 out << nPoints << " 3\n";
706 // the point list
707 for (std::size_t i = 0; i < nPoints; ++i)
708 {
709 out << i << " " << (*nodes[i])[0] << " " << (*nodes[i])[1] << " "
710 << (*nodes[i])[2] << "\n";
711 }
712
713 if (mesh.getDimension() == 2)
714 {
715 write2dElements(out, mesh);
716 }
717 else
718 {
719 write3dElements(out, mesh, attribute_points);
720 }
721
722 out << "0\n"; // the polygon holes list
723
724 // the region attributes list
725 if (attribute_points.empty())
726 {
727 out << "0\n";
728 }
729 else
730 {
731 const std::size_t nAttributePoints(attribute_points.size());
732 out << nAttributePoints << "\n";
733 for (std::size_t i = 0; i < nAttributePoints; ++i)
734 {
735 out << i + 1 << " " << attribute_points[i][0] << " "
736 << attribute_points[i][1] << " " << attribute_points[i][2]
737 << " " << 10 * attribute_points[i].getID() << "\n";
738 }
739 }
740
741 INFO(
742 "TetGenInterface::writeTetGenPoly() - {:d} points and {:d} surfaces "
743 "successfully written.",
744 nPoints,
745 mesh.getNumberOfElements());
746 out.close();
747 return true;
748}
749
750void TetGenInterface::write2dElements(std::ofstream& out,
751 const MeshLib::Mesh& mesh) const
752{
753 // the surfaces header
754 auto const& types =
756 std::size_t const n_tri =
757 (types.find(MeshLib::MeshElemType::TRIANGLE) != types.end())
759 : 0;
760 std::size_t const n_quad =
761 (types.find(MeshLib::MeshElemType::QUAD) != types.end())
763 : 0;
764 unsigned const nTotalTriangles = n_tri + 2 * n_quad;
765 out << nTotalTriangles << " 1\n";
766
767 const std::vector<MeshLib::Element*>& elements = mesh.getElements();
768 MeshLib::PropertyVector<int> const* const mat_ids =
769 mesh.getProperties().existsPropertyVector<int>("MaterialIDs")
770 ? mesh.getProperties().getPropertyVector<int>("MaterialIDs")
771 : nullptr;
772 const std::size_t nElements(elements.size());
773 unsigned element_count(0);
774 for (std::size_t i = 0; i < nElements; ++i)
775 {
776 std::string const matId = mat_ids ? std::to_string((*mat_ids)[i]) : "";
777 this->writeElementToFacets(out, *elements[i], element_count, matId);
778 }
779}
780
782 std::ofstream& out,
783 const MeshLib::Mesh& mesh,
784 std::vector<MeshLib::Node>& attribute_points) const
785{
786 const std::vector<MeshLib::Element*>& elements = mesh.getElements();
787 const std::size_t nElements(elements.size());
788 if (!attribute_points.empty())
789 {
790 attribute_points.clear();
791 }
792
793 // get position where number of facets need to be written and figure out
794 // worst case of chars that are needed
795 const std::streamoff before_elems_pos(out.tellp());
796 const unsigned n_spaces(
797 static_cast<unsigned>(std::floor(log(nElements * 8))) + 1);
798 out << std::string(n_spaces, ' ') << " 1\n";
799 auto const* const materialIds =
800 mesh.getProperties().getPropertyVector<int>("MaterialIDs");
801 unsigned element_count(0);
802 for (std::size_t i = 0; i < nElements; ++i)
803 {
804 if (elements[i]->getDimension() < 3)
805 {
806 continue;
807 }
808
809 const unsigned nFaces(elements[i]->getNumberOfNeighbors());
810 std::string const mat_id_str =
811 materialIds ? std::to_string((*materialIds)[i]) : "";
812 for (std::size_t j = 0; j < nFaces; ++j)
813 {
814 MeshLib::Element const* const neighbor(elements[i]->getNeighbor(j));
815
816 if (neighbor && materialIds &&
817 (*materialIds)[i] <= (*materialIds)[neighbor->getID()])
818 {
819 continue;
820 }
821
822 std::unique_ptr<MeshLib::Element const> const face(
823 elements[i]->getFace(j));
824 this->writeElementToFacets(out, *face, element_count, mat_id_str);
825 }
826 if (materialIds)
827 {
828 attribute_points.emplace_back(
829 MeshLib::getCenterOfGravity(*elements[i]).data(),
830 (*materialIds)[i]);
831 }
832 }
833 // add number of facets at correct position and jump back
834 const std::streamoff after_elems_pos(out.tellp());
835 out.seekp(before_elems_pos);
836 out << element_count;
837 out.seekp(after_elems_pos);
838}
839
841 const MeshLib::Element& element,
842 unsigned& element_count,
843 std::string const& matId)
844{
845 element_count++;
847 {
848 out << "3 " << getNodeIndex(element, 0) << " "
849 << getNodeIndex(element, 1) << " " << getNodeIndex(element, 2)
850 << " " << matId << " # " << element_count << "\n";
851 }
852 else if (element.getGeomType() == MeshLib::MeshElemType::QUAD)
853 {
854 out << "3 " << getNodeIndex(element, 0) << " "
855 << getNodeIndex(element, 1) << " " << getNodeIndex(element, 2)
856 << " " << matId << " # " << element_count << "\n";
857 element_count++;
858 out << "3 " << getNodeIndex(element, 0) << " "
859 << getNodeIndex(element, 2) << " " << getNodeIndex(element, 3)
860 << " " << matId << " # " << element_count << "\n";
861 }
862}
863
864} // end namespace FileIO
Definition of the Element class.
Filename manipulation routines.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:40
Definition of the MeshInformation class.
Definition of the Mesh class.
Definition of the Node class.
Definition of string helper functions.
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 &region_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.
Definition: GEOObjects.h:61
const std::vector< Point * > * getPointVec(const std::string &name) const
Definition: GEOObjects.cpp:72
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()))
Definition: GEOObjects.cpp:47
const PointVec * getPointVecObj(const std::string &name) const
Definition: GEOObjects.cpp:85
const std::vector< Surface * > * getSurfaceVec(const std::string &name) const
Returns the surface vector with the given name as a const.
Definition: GEOObjects.cpp:270
void addSurfaceVec(std::vector< Surface * > &&sfc, const std::string &name, SurfaceVec::NameIdMap &&sfc_names)
Definition: GEOObjects.cpp:240
const std::vector< std::size_t > & getIDMap() const
Definition: PointVec.h:96
A Surface is represented by Triangles. It consists of a reference to a vector of (pointers to) points...
Definition: Surface.h:35
std::map< std::string, std::size_t > NameIdMap
Definition: TemplateVec.h:43
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
Definition: Triangle.h:27
virtual MeshElemType getGeomType() const =0
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:102
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:105
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:84
Properties & getProperties()
Definition: Mesh.h:130
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:93
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition: Properties.h:36
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
static std::map< MeshLib::MeshElemType, unsigned > getNumberOfElementTypes(const MeshLib::Mesh &mesh)
void simplify(std::string &str)
Definition: StringTools.cpp:77
std::string getFileExtension(const std::string &path)
Definition: FileTools.cpp:193
void cleanupVectorElements(std::vector< T * > &items)
Definition: Algorithm.h:241
std::string extractBaseNameWithoutExtension(std::string const &pathname)
Definition: FileTools.cpp:187
std::vector< std::string > splitString(std::string const &str)
Definition: StringTools.cpp:29
auto constructPointsFromNodes(std::vector< MeshLib::Node * > nodes)
MathLib::Point3d getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.
Definition: Element.cpp:124