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