30 SmartPointer<Mesh<T>> mesh =
nullptr;
34 unsigned vtk_nodes_for_cell_type[15] = {0, 1, 0, 2, 0, 3, 0, 0,
43 : mesh(passedMesh), fileName(std::move(passedFileName)) {}
46 std::string passedFileName)
47 : mesh(passedMesh), fileFormat(passedFormat),
48 fileName(std::move(passedFileName)) {}
58 fileName = std::move(passedFileName);
63 if (mesh ==
nullptr) {
65 .addWarning(
"No mesh was passed to VTKReader. Not reading.")
70 if (fileName.empty()) {
72 .addWarning(
"No file name specified for VTKReader. Not reading.")
78 auto dotPos = fileName.rfind(
'.');
79 if (dotPos == std::string::npos) {
81 .addWarning(
"No valid file format found based on the file ending "
82 "passed to VTKReader. Not reading.")
86 auto ending = fileName.substr(dotPos);
87 if (ending ==
".vtk") {
89 }
else if (ending ==
".vtp") {
91 }
else if (ending ==
".vtu") {
95 .addWarning(
"No valid file format found based on the file ending "
96 "passed to VTKReader. Not reading.")
103 switch (fileFormat) {
105 readVTKLegacy(fileName);
107#ifdef VIENNALS_USE_VTK
117 Logger::getInstance()
119 "VTKReader was built without VTK support. Only VTK_LEGACY "
120 "can be used. File not read.")
124 Logger::getInstance()
125 .addWarning(
"No valid file format set for VTKReader. Not reading.")
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.")
141 vtkSmartPointer<vtkXMLPolyDataReader> pReader =
142 vtkSmartPointer<vtkXMLPolyDataReader>::New();
143 pReader->SetFileName(filename.c_str());
146 vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
147 polyData = pReader->GetOutput();
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;
156 vtkSmartPointer<vtkCellArray> cellArray =
157 vtkSmartPointer<vtkCellArray>::New();
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);
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);
182 mesh->lines.push_back(cell);
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);
197 mesh->triangles.push_back(cell);
202 vtkSmartPointer<vtkPointData> pointData =
203 vtkSmartPointer<vtkPointData>::New();
204 pointData = polyData->GetPointData();
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);
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]));
232 vtkSmartPointer<vtkCellData> cellData = vtkSmartPointer<vtkCellData>::New();
233 cellData = polyData->GetCellData();
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);
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]));
261 void readVTU(
const std::string &filename) {
262 if (mesh ==
nullptr) {
263 Logger::getInstance()
264 .addWarning(
"No mesh was passed to VTKReader.")
271 vtkSmartPointer<vtkXMLUnstructuredGridReader> greader =
272 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
273 greader->SetFileName(filename.c_str());
276 vtkSmartPointer<vtkUnstructuredGrid> ugrid =
277 vtkSmartPointer<vtkUnstructuredGrid>::New();
278 ugrid = greader->GetOutput();
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;
289 for (
unsigned i = 0; i < ugrid->GetNumberOfCells(); ++i) {
290 vtkIdList *pointList = vtkIdList::New();
291 ugrid->GetCellPoints(i, pointList);
293 switch (ugrid->GetCellType(i)) {
296 std::array<unsigned, 1> vert{};
297 vert[0] = pointList->GetId(0);
298 mesh->vertices.push_back(vert);
302 std::array<unsigned, 2> elements{};
303 for (
unsigned j = 0; j < 2; ++j) {
304 elements[j] = pointList->GetId(j);
306 mesh->lines.push_back(elements);
310 std::array<unsigned, 3> elements{};
311 for (
unsigned j = 0; j < 3; ++j) {
312 elements[j] = pointList->GetId(j);
314 mesh->triangles.push_back(elements);
318 std::array<unsigned, 4> elements{};
319 for (
unsigned j = 0; j < 4; ++j) {
320 elements[j] = pointList->GetId(j);
322 mesh->tetras.push_back(elements);
326 std::array<unsigned, 8> elements{};
327 for (
unsigned j = 0; j < 8; ++j) {
328 elements[j] = pointList->GetId(j);
330 mesh->hexas.push_back(elements);
336 vtkSmartPointer<vtkPointData> pointData =
337 vtkSmartPointer<vtkPointData>::New();
338 pointData = ugrid->GetPointData();
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);
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]));
366 vtkSmartPointer<vtkCellData> cellData = vtkSmartPointer<vtkCellData>::New();
367 cellData = ugrid->GetCellData();
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);
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]));
397 void readVTKLegacy(
const std::string &filename) {
398 if (mesh ==
nullptr) {
399 Logger::getInstance()
400 .addWarning(
"No mesh was passed to VTKReader.")
407 std::ifstream f(filename.c_str());
409 Logger::getInstance().addError(
"Could not open geometry file!");
413 while (std::getline(f, temp)) {
414 if (temp.find(
"DATASET") != std::string::npos)
417 if (temp.find(
"UNSTRUCTURED_GRID") == std::string::npos) {
418 Logger::getInstance().addError(
"DATASET is not an UNSTRUCTURED_GRID!");
422 while (std::getline(f, temp)) {
423 if (temp.find(
"POINTS") != std::string::npos)
426 int num_nodes = atoi(&temp[temp.find(
' ') + 1]);
428 mesh->nodes.resize(num_nodes);
430 for (
int i = 0; i < num_nodes; i++) {
433 for (
double &coord : coords)
435 for (
int j = 0; j < 3; j++) {
436 mesh->nodes[i][j] = coords[j];
448 while (std::getline(f, temp)) {
449 if (temp.find(
"CELLS") == 0)
453 int num_elems = atoi(&temp[temp.find(
' ') + 1]);
455 std::ifstream f_ct(filename.c_str());
460 while (std::getline(f_ct, temp)) {
461 if (temp.find(
"CELL_TYPES") == 0)
464 int num_cell_types = atoi(&temp[temp.find(
' ') + 1]);
466 if (num_elems != num_cell_types) {
467 Logger::getInstance().addError(
468 "Corrupt input geometry! Number of CELLS and CELL_TYPES "
472 bool is_material =
true;
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);
489 std::vector<VectorType<unsigned int, 4>> elements;
490 elements.reserve(num_elems);
492 std::vector<double> materials;
496 unsigned cell_material;
497 for (
int i = 0; i < num_elems; i++) {
501 f_m >> cell_material;
507 unsigned number_nodes = vtk_nodes_for_cell_type[cell_type];
508 if (number_nodes == elems_fake || number_nodes == 0) {
512 std::array<unsigned, 1> elem{};
514 mesh->template getElements<1>().push_back(elem);
515 materials.push_back(cell_material);
519 std::array<unsigned, 2> elem{};
520 for (
unsigned j = 0; j < number_nodes; ++j) {
523 mesh->template getElements<2>().push_back(elem);
524 materials.push_back(cell_material);
529 std::array<unsigned, 3> elem{};
530 for (
unsigned j = 0; j < number_nodes; ++j) {
533 mesh->template getElements<3>().push_back(elem);
534 materials.push_back(cell_material);
540 std::array<unsigned, 4> elem{};
541 for (
unsigned j = 0; j < number_nodes; ++j) {
544 mesh->template getElements<4>().push_back(elem);
545 materials.push_back(cell_material);
551 std::array<unsigned, 3> elem{};
552 for (
unsigned j = 0; j < 3; ++j) {
555 mesh->template getElements<3>().push_back(
557 materials.push_back(cell_material);
560 mesh->template getElements<3>().push_back(elem);
561 materials.push_back(cell_material);
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();
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());
578 f.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
582 mesh->cellData.insertNextScalarData(materials,
"Material");
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]);
592 std::cout <<
"Read cell data: " << num_cell_data << std::endl;
595 std::cout <<
"reading scalar data" << std::endl;
596 while (std::getline(f, temp)) {
597 if (temp.find(
"SCALARS") != std::string::npos)
604 std::string scalarDataName;
606 auto firstS = temp.find(
' ') + 1;
607 auto secondS = temp.find(
' ', firstS + 1);
608 scalarDataName = temp.substr(firstS, secondS - firstS);
610 std::vector<double> scalarData;
613 std::getline(f, temp);
614 if (temp !=
"LOOKUP_TABLE default") {
615 Logger::getInstance()
616 .addWarning(
"Wrong lookup table for VTKLegacy: " + temp)
621 for (
int i = 0; i < num_cell_data; ++i) {
624 scalarData.push_back(data);
627 mesh->cellData.insertNextScalarData(scalarData, scalarDataName);