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