ViennaLS
Loading...
Searching...
No Matches
lsVTKReader.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <fstream>
4#include <string>
5#include <unordered_map>
6
7#include <lsFileFormats.hpp>
8#include <lsMesh.hpp>
9
10#include <utility>
11#include <vcLogger.hpp>
12#include <vcSmartPointer.hpp>
13
14#ifdef VIENNALS_USE_VTK
15#include <vtkCellData.h>
16#include <vtkIdList.h>
17#include <vtkPointData.h>
18#include <vtkPolyData.h>
19#include <vtkSmartPointer.h>
20#include <vtkUnstructuredGrid.h>
21#include <vtkXMLPolyDataReader.h>
22#include <vtkXMLUnstructuredGridReader.h>
23#endif // VIENNALS_USE_VTK
24
25namespace viennals {
26
27using namespace viennacore;
28
30template <class T = double> class VTKReader {
31 SmartPointer<Mesh<T>> mesh = nullptr;
33 std::string fileName;
34 std::unordered_map<std::string, std::vector<T>> metaData;
35
36 unsigned vtk_nodes_for_cell_type[15] = {0, 1, 0, 2, 0, 3, 0, 0,
37 4, 4, 4, 8, 8, 6, 5};
38
39 // accepts either vtkPolyData or vtkUnstructuredGrid
40 void extractFieldData(vtkDataSet *data) {
41 vtkFieldData *fieldData = data->GetFieldData();
42 if (!fieldData)
43 return;
44
45 for (int i = 0; i < fieldData->GetNumberOfArrays(); ++i) {
46 vtkDataArray *array = fieldData->GetArray(i);
47 if (!array)
48 continue;
49
50 const char *name = array->GetName();
51 if (!name)
52 continue;
53
54 int numTuples = array->GetNumberOfTuples();
55 int numComponents = array->GetNumberOfComponents();
56
57 std::vector<T> values;
58 values.reserve(numTuples * numComponents);
59
60 for (int t = 0; t < numTuples; ++t) {
61 T tuple[9]; // safe default (up to 9 components)
62 array->GetTuple(t, tuple);
63 for (int c = 0; c < numComponents; ++c)
64 values.push_back(tuple[c]);
65 }
66
67 metaData[name] = std::move(values);
68 }
69 }
70
71public:
72 VTKReader() = default;
73
74 VTKReader(SmartPointer<Mesh<T>> passedMesh) : mesh(passedMesh) {}
75
76 VTKReader(SmartPointer<Mesh<T>> passedMesh, std::string passedFileName)
77 : mesh(passedMesh), fileName(std::move(passedFileName)) {}
78
79 VTKReader(SmartPointer<Mesh<>> passedMesh, FileFormatEnum passedFormat,
80 std::string passedFileName)
81 : mesh(passedMesh), fileFormat(passedFormat),
82 fileName(std::move(passedFileName)) {}
83
85 void setMesh(SmartPointer<Mesh<>> passedMesh) { mesh = passedMesh; }
86
88 void setFileFormat(FileFormatEnum passedFormat) { fileFormat = passedFormat; }
89
91 void setFileName(std::string passedFileName) {
92 fileName = std::move(passedFileName);
93 }
94
95 auto &getMetaData() { return metaData; }
96
97 void apply() {
98 // check mesh
99 if (mesh == nullptr) {
100 Logger::getInstance()
101 .addWarning("No mesh was passed to VTKReader. Not reading.")
102 .print();
103 return;
104 }
105 // check filename
106 if (fileName.empty()) {
107 Logger::getInstance()
108 .addWarning("No file name specified for VTKReader. Not reading.")
109 .print();
110 return;
111 }
112
113 if (fileFormat == FileFormatEnum::VTK_AUTO) {
114 auto dotPos = fileName.rfind('.');
115 if (dotPos == std::string::npos) {
116 Logger::getInstance()
117 .addWarning("No valid file format found based on the file ending "
118 "passed to VTKReader. Not reading.")
119 .print();
120 return;
121 }
122 auto ending = fileName.substr(dotPos);
123 if (ending == ".vtk") {
124 fileFormat = FileFormatEnum::VTK_LEGACY;
125 } else if (ending == ".vtp") {
126 fileFormat = FileFormatEnum::VTP;
127 } else if (ending == ".vtu") {
128 fileFormat = FileFormatEnum::VTU;
129 } else {
130 Logger::getInstance()
131 .addWarning("No valid file format found based on the file ending "
132 "passed to VTKReader. Not reading.")
133 .print();
134 return;
135 }
136 }
137
138 // check file format
139 switch (fileFormat) {
141 readVTKLegacy(fileName);
142 break;
143#ifdef VIENNALS_USE_VTK
145 readVTP(fileName);
146 break;
148 readVTU(fileName);
149 break;
150#else
153 Logger::getInstance()
154 .addWarning(
155 "VTKReader was built without VTK support. Only VTK_LEGACY "
156 "can be used. File not read.")
157 .print();
158#endif
159 default:
160 Logger::getInstance()
161 .addWarning("No valid file format set for VTKReader. Not reading.")
162 .print();
163 }
164 }
165
166private:
167#ifdef VIENNALS_USE_VTK
168 void readVTP(const std::string &filename) {
169
170 mesh->clear();
171 vtkSmartPointer<vtkXMLPolyDataReader> pReader =
172 vtkSmartPointer<vtkXMLPolyDataReader>::New();
173 pReader->SetFileName(filename.c_str());
174 pReader->Update();
175
176 vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
177 polyData = pReader->GetOutput();
178
179 mesh->nodes.resize(polyData->GetNumberOfPoints());
180 for (unsigned i = 0; i < mesh->nodes.size(); ++i) {
181 std::array<double, 3> coords{};
182 polyData->GetPoint(i, coords.data());
183 mesh->nodes[i] = coords;
184 }
185
186 vtkSmartPointer<vtkCellArray> cellArray =
187 vtkSmartPointer<vtkCellArray>::New();
188 // get vertices
189 {
190 mesh->vertices.reserve(polyData->GetNumberOfVerts());
191 cellArray = polyData->GetVerts();
192 cellArray->InitTraversal();
193 vtkIdList *pointList = vtkIdList::New();
194 while (cellArray->GetNextCell(pointList)) {
195 std::array<unsigned, 1> cell{};
196 cell[0] = pointList->GetId(0);
197 mesh->vertices.push_back(cell);
198 }
199 }
200
201 // get lines
202 {
203 mesh->lines.reserve(polyData->GetNumberOfLines());
204 cellArray = polyData->GetLines();
205 cellArray->InitTraversal();
206 vtkIdList *pointList = vtkIdList::New();
207 while (cellArray->GetNextCell(pointList)) {
208 std::array<unsigned, 2> cell{};
209 for (unsigned i = 0; i < 2; ++i) {
210 cell[i] = pointList->GetId(i);
211 }
212 mesh->lines.push_back(cell);
213 }
214 }
215
216 // get triangles
217 {
218 mesh->triangles.reserve(polyData->GetNumberOfPolys());
219 cellArray = polyData->GetPolys();
220 cellArray->InitTraversal();
221 vtkIdList *pointList = vtkIdList::New();
222 while (cellArray->GetNextCell(pointList)) {
223 std::array<unsigned, 3> cell{};
224 for (unsigned i = 0; i < 3; ++i) {
225 cell[i] = pointList->GetId(i);
226 }
227 mesh->triangles.push_back(cell);
228 }
229 }
230
231 // get point data
232 vtkSmartPointer<vtkPointData> pointData =
233 vtkSmartPointer<vtkPointData>::New();
234 pointData = polyData->GetPointData();
235
236 for (int i = 0; i < static_cast<int>(pointData->GetNumberOfArrays()); ++i) {
237 if (vtkDataArray *dataArray = pointData->GetArray(i);
238 dataArray->GetNumberOfComponents() == 1) {
239 mesh->pointData.insertNextScalarData(
241 std::string(pointData->GetArrayName(i)));
242 auto &scalars = *(mesh->pointData.getScalarData(i));
243 scalars.resize(pointData->GetNumberOfTuples());
244 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
245 scalars[j] = dataArray->GetTuple1(j);
246 }
247 } else if (dataArray->GetNumberOfComponents() == 3) {
248 mesh->pointData.insertNextVectorData(
250 std::string(pointData->GetArrayName(i)));
251 auto &vectors = *(mesh->pointData.getVectorData(i));
252 vectors.resize(pointData->GetNumberOfTuples());
253 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
254 std::array<double, 3> vector{};
255 dataArray->GetTuple(j, &(vector[0]));
256 vectors[j] = vector;
257 }
258 }
259 }
260
261 // get cell data
262 vtkSmartPointer<vtkCellData> cellData = vtkSmartPointer<vtkCellData>::New();
263 cellData = polyData->GetCellData();
264
265 for (int i = 0; i < cellData->GetNumberOfArrays(); ++i) {
266 vtkDataArray *dataArray = cellData->GetArray(i);
267 if (cellData->GetNumberOfComponents() == 1) {
268 mesh->cellData.insertNextScalarData(
270 std::string(cellData->GetArrayName(i)));
271 auto &scalars = *(mesh->cellData.getScalarData(i));
272 scalars.resize(cellData->GetNumberOfTuples());
273 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
274 scalars[j] = dataArray->GetTuple1(j);
275 }
276 } else if (cellData->GetNumberOfComponents() == 3) {
277 mesh->cellData.insertNextVectorData(
279 std::string(cellData->GetArrayName(i)));
280 auto &vectors = *(mesh->cellData.getVectorData(i));
281 vectors.resize(cellData->GetNumberOfTuples());
282 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
283 std::array<double, 3> vector{};
284 dataArray->GetTuple(j, &(vector[0]));
285 vectors[j] = vector;
286 }
287 }
288 }
289
290 // read meta data
291 extractFieldData(polyData);
292 }
293
294 void readVTU(const std::string &filename) {
295
296 mesh->clear();
297
298 vtkSmartPointer<vtkXMLUnstructuredGridReader> greader =
299 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
300 greader->SetFileName(filename.c_str());
301 greader->Update();
302
303 vtkSmartPointer<vtkUnstructuredGrid> ugrid =
304 vtkSmartPointer<vtkUnstructuredGrid>::New();
305 ugrid = greader->GetOutput();
306
307 // get all points
308 mesh->nodes.resize(ugrid->GetNumberOfPoints());
309 for (unsigned i = 0; i < mesh->nodes.size(); ++i) {
310 std::array<double, 3> coords{};
311 ugrid->GetPoint(i, &(coords[0]));
312 mesh->nodes[i] = coords;
313 }
314
315 // get cells
316 for (unsigned i = 0; i < ugrid->GetNumberOfCells(); ++i) {
317 vtkIdList *pointList = vtkIdList::New();
318 ugrid->GetCellPoints(i, pointList);
319
320 switch (ugrid->GetCellType(i)) {
321 case 1: // vert
322 {
323 std::array<unsigned, 1> vert{};
324 vert[0] = pointList->GetId(0);
325 mesh->vertices.push_back(vert);
326 } break;
327 case 3: // line
328 {
329 std::array<unsigned, 2> elements{};
330 for (unsigned j = 0; j < 2; ++j) {
331 elements[j] = pointList->GetId(j);
332 }
333 mesh->lines.push_back(elements);
334 } break;
335 case 5: // triangle
336 {
337 std::array<unsigned, 3> elements{};
338 for (unsigned j = 0; j < 3; ++j) {
339 elements[j] = pointList->GetId(j);
340 }
341 mesh->triangles.push_back(elements);
342 } break;
343 case 10: // tetra
344 {
345 std::array<unsigned, 4> elements{};
346 for (unsigned j = 0; j < 4; ++j) {
347 elements[j] = pointList->GetId(j);
348 }
349 mesh->tetras.push_back(elements);
350 } break;
351 case 12: // hexa
352 {
353 std::array<unsigned, 8> elements{};
354 for (unsigned j = 0; j < 8; ++j) {
355 elements[j] = pointList->GetId(j);
356 }
357 mesh->hexas.push_back(elements);
358 } break;
359 }
360 }
361
362 // get point data
363 vtkSmartPointer<vtkPointData> pointData =
364 vtkSmartPointer<vtkPointData>::New();
365 pointData = ugrid->GetPointData();
366
367 for (int i = 0; i < pointData->GetNumberOfArrays(); ++i) {
368 if (vtkDataArray *dataArray = pointData->GetArray(i);
369 dataArray->GetNumberOfComponents() == 1) {
370 mesh->pointData.insertNextScalarData(
372 std::string(pointData->GetArrayName(i)));
373 auto &scalars = *(mesh->pointData.getScalarData(i));
374 scalars.resize(pointData->GetNumberOfTuples());
375 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
376 scalars[j] = dataArray->GetTuple1(j);
377 }
378 } else if (dataArray->GetNumberOfComponents() == 3) {
379 mesh->pointData.insertNextVectorData(
381 std::string(pointData->GetArrayName(i)));
382 auto &vectors = *(mesh->pointData.getVectorData(i));
383 vectors.resize(pointData->GetNumberOfTuples());
384 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
385 std::array<double, 3> vector{};
386 dataArray->GetTuple(j, &(vector[0]));
387 vectors[j] = vector;
388 }
389 }
390 }
391
392 // get cell data
393 vtkSmartPointer<vtkCellData> cellData = vtkSmartPointer<vtkCellData>::New();
394 cellData = ugrid->GetCellData();
395
396 for (int i = 0; i < static_cast<int>(cellData->GetNumberOfArrays()); ++i) {
397 vtkDataArray *dataArray = cellData->GetArray(i);
398 if (cellData->GetNumberOfComponents() == 1) {
399 mesh->cellData.insertNextScalarData(
401 std::string(cellData->GetArrayName(i)));
402 auto &scalars = *(mesh->cellData.getScalarData(i));
403 scalars.resize(cellData->GetNumberOfTuples());
404 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
405 scalars[j] = dataArray->GetTuple1(j);
406 }
407 } else if (cellData->GetNumberOfComponents() == 3) {
408 mesh->cellData.insertNextVectorData(
410 std::string(cellData->GetArrayName(i)));
411 auto &vectors = *(mesh->cellData.getVectorData(i));
412 vectors.resize(cellData->GetNumberOfTuples());
413 for (unsigned j = 0; j < dataArray->GetNumberOfTuples(); ++j) {
414 std::array<double, 3> vector{};
415 dataArray->GetTuple(j, &(vector[0]));
416 vectors[j] = vector;
417 }
418 }
419 }
420
421 // read meta data
422 extractFieldData(ugrid);
423 }
424
425#endif // VIENNALS_USE_VTK
426
427 void readVTKLegacy(const std::string &filename) {
428
429 mesh->clear();
430 // open geometry file
431 std::ifstream f(filename.c_str());
432 if (!f)
433 Logger::getInstance().addError("Could not open geometry file!");
434 std::string temp;
435
436 // Check if geometry is an unstructured grid as required
437 while (std::getline(f, temp)) {
438 if (temp.find("DATASET") != std::string::npos)
439 break;
440 }
441 if (temp.find("UNSTRUCTURED_GRID") == std::string::npos) {
442 Logger::getInstance().addError("DATASET is not an UNSTRUCTURED_GRID!");
443 }
444
445 // Find POINTS in file to know number of nodes to read in
446 while (std::getline(f, temp)) {
447 if (temp.find("POINTS") != std::string::npos)
448 break;
449 }
450 int num_nodes = atoi(&temp[temp.find(' ') + 1]);
451
452 mesh->nodes.resize(num_nodes);
453
454 for (int i = 0; i < num_nodes; i++) {
455 double coords[3];
456
457 for (double &coord : coords)
458 f >> coord;
459 for (int j = 0; j < 3; j++) {
460 mesh->nodes[i][j] = coords[j];
461 // int shift_size = shift.size();
462 // if (shift_size > j)
463 // mesh->nodes[i][j] += shift[j]; // Assign desired shift
464 // if (InputTransformationSigns[j])
465 // mesh->nodes[i][j] = -mesh->nodes[i][j]; // Assign sign
466 // transformation, if needed
467 // mesh->nodes[i][j] *= scale; // Scale the geometry according to
468 // parameters file
469 }
470 }
471
472 while (std::getline(f, temp)) {
473 if (temp.find("CELLS") == 0)
474 break;
475 }
476
477 int num_elems = atoi(&temp[temp.find(' ') + 1]);
478
479 std::ifstream f_ct(filename.c_str()); // stream to read cell CELL_TYPES
480 std::ifstream f_m(
481 filename.c_str()); // stream for material numbers if they exist
482
483 // advance to cell types and check if there are the right number
484 while (std::getline(f_ct, temp)) {
485 if (temp.find("CELL_TYPES") == 0)
486 break;
487 }
488 int num_cell_types = atoi(&temp[temp.find(' ') + 1]);
489 // need a cell_type for each cell
490 if (num_elems != num_cell_types) {
491 Logger::getInstance().addError(
492 "Corrupt input geometry! Number of CELLS and CELL_TYPES "
493 "is different!");
494 }
495
496 bool is_material = true;
497 // advance to material if it is specified
498 while (std::getline(f_m, temp)) {
499 if (temp.find("CELL_DATA") != std::string::npos) {
500 std::getline(f_m, temp);
501 if ((temp.find("SCALARS material") != std::string::npos) ||
502 (temp.find("SCALARS Material") != std::string::npos)) {
503 std::getline(f_m, temp);
504 break;
505 }
506 }
507 }
508 if (f_m.eof()) {
509 is_material = false;
510 }
511
512 // maximum cell to read in is a tetra with 4 points
513 std::vector<VectorType<unsigned int, 4>> elements;
514 elements.reserve(num_elems);
515
516 std::vector<double> materials;
517
518 unsigned elems_fake;
519 unsigned cell_type;
520 unsigned cell_material;
521 for (int i = 0; i < num_elems; i++) {
522 f >> elems_fake;
523 f_ct >> cell_type;
524 if (is_material)
525 f_m >> cell_material;
526 else
527 cell_material =
528 1; // if there are no materials specified make all the same
529
530 // check if the correct number of nodes for cell_type is given
531 unsigned number_nodes = vtk_nodes_for_cell_type[cell_type];
532 if (number_nodes == elems_fake || number_nodes == 0) {
533 // check for different types to subdivide them into supported types
534 switch (cell_type) {
535 case 1: {
536 std::array<unsigned, 1> elem{};
537 f >> elem[0];
538 mesh->template getElements<1>().push_back(elem);
539 materials.push_back(cell_material);
540 break;
541 }
542 case 3: {
543 std::array<unsigned, 2> elem{};
544 for (unsigned j = 0; j < number_nodes; ++j) {
545 f >> elem[j];
546 }
547 mesh->template getElements<2>().push_back(elem);
548 materials.push_back(cell_material);
549 break;
550 }
551 case 5: // triangle for 2D
552 {
553 std::array<unsigned, 3> elem{};
554 for (unsigned j = 0; j < number_nodes; ++j) {
555 f >> elem[j];
556 }
557 mesh->template getElements<3>().push_back(elem);
558 materials.push_back(cell_material);
559 break;
560 }
561
562 case 10: // tetra for 3D
563 {
564 std::array<unsigned, 4> elem{};
565 for (unsigned j = 0; j < number_nodes; ++j) {
566 f >> elem[j];
567 }
568 mesh->template getElements<4>().push_back(elem);
569 materials.push_back(cell_material);
570 break;
571 }
572
573 case 9: // this is a quad, so just plit it into two triangles
574 {
575 std::array<unsigned, 3> elem{};
576 for (unsigned j = 0; j < 3; ++j) {
577 f >> elem[j];
578 }
579 mesh->template getElements<3>().push_back(
580 elem); // push the first three nodes as a triangle
581 materials.push_back(cell_material);
582
583 f >> elem[1]; // replace middle element to create other triangle
584 mesh->template getElements<3>().push_back(elem);
585 materials.push_back(cell_material);
586 break;
587 }
588
589 default:
590 std::ostringstream oss;
591 oss << "VTK Cell type " << cell_type
592 << " is not supported. Cell ignored..." << std::endl;
593 Logger::getInstance().addWarning(oss.str()).print();
594 }
595 } else {
596 std::ostringstream oss;
597 oss << "INVALID CELL TYPE! Expected number of nodes: " << number_nodes
598 << ", Found number of nodes: " << elems_fake
599 << "; Ignoring element...";
600 Logger::getInstance().addError(oss.str());
601 // ignore rest of lines
602 f.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
603 }
604 }
605
606 mesh->cellData.insertNextScalarData(materials, "Material");
607
608 // Now read Cell Data
609 int num_cell_data = 0;
610 while (std::getline(f, temp)) {
611 if (temp.find("CELL_DATA") != std::string::npos) {
612 num_cell_data = atoi(&temp[temp.find(' ') + 1]);
613 break;
614 }
615 }
616 std::cout << "Read cell data: " << num_cell_data << std::endl;
617
618 while (!f.eof()) {
619 std::cout << "reading scalar data" << std::endl;
620 while (std::getline(f, temp)) {
621 if (temp.find("SCALARS") != std::string::npos)
622 break;
623 }
624 if (f.eof()) {
625 break;
626 }
627
628 std::string scalarDataName;
629 {
630 auto firstS = temp.find(' ') + 1;
631 auto secondS = temp.find(' ', firstS + 1);
632 scalarDataName = temp.substr(firstS, secondS - firstS);
633 }
634 std::vector<double> scalarData;
635
636 // consume one line, which defines the lookup table
637 std::getline(f, temp);
638 if (temp != "LOOKUP_TABLE default") {
639 Logger::getInstance()
640 .addWarning("Wrong lookup table for VTKLegacy: " + temp)
641 .print();
642 }
643
644 // now read scalar values
645 for (int i = 0; i < num_cell_data; ++i) {
646 double data;
647 f >> data;
648 scalarData.push_back(data);
649 }
650
651 mesh->cellData.insertNextScalarData(scalarData, scalarDataName);
652 }
653
654 f_ct.close();
655 f_m.close();
656 f.close();
657 }
658};
659
660} // namespace viennals
double T
Definition Epitaxy.cpp:10
This class holds an explicit mesh, which is always given in 3 dimensions. If it describes a 2D mesh,...
Definition lsMesh.hpp:22
std::vector< Vec3D< T > > VectorDataType
Definition lsPointData.hpp:24
std::vector< T > ScalarDataType
Definition lsPointData.hpp:23
void setFileName(std::string passedFileName)
set file name for file to read
Definition lsVTKReader.hpp:91
void setMesh(SmartPointer< Mesh<> > passedMesh)
set the mesh the file should be read into
Definition lsVTKReader.hpp:85
VTKReader(SmartPointer< Mesh< T > > passedMesh, std::string passedFileName)
Definition lsVTKReader.hpp:76
VTKReader(SmartPointer< Mesh< T > > passedMesh)
Definition lsVTKReader.hpp:74
VTKReader(SmartPointer< Mesh<> > passedMesh, FileFormatEnum passedFormat, std::string passedFileName)
Definition lsVTKReader.hpp:79
void setFileFormat(FileFormatEnum passedFormat)
set file format for file to read. Defaults to VTK_LEGACY.
Definition lsVTKReader.hpp:88
void apply()
Definition lsVTKReader.hpp:97
auto & getMetaData()
Definition lsVTKReader.hpp:95
Definition lsAdvect.hpp:36
FileFormatEnum
Definition lsFileFormats.hpp:4
@ VTK_AUTO
Definition lsFileFormats.hpp:8
@ VTK_LEGACY
Definition lsFileFormats.hpp:5
@ VTP
Definition lsFileFormats.hpp:6
@ VTU
Definition lsFileFormats.hpp:7