31 SmartPointer<Mesh<T>> mesh =
nullptr;
34 std::unordered_map<std::string, std::vector<T>> metaData;
36 unsigned vtk_nodes_for_cell_type[15] = {0, 1, 0, 2, 0, 3, 0, 0,
40 void extractFieldData(vtkDataSet *data) {
41 vtkFieldData *fieldData = data->GetFieldData();
45 for (
int i = 0; i < fieldData->GetNumberOfArrays(); ++i) {
46 vtkDataArray *array = fieldData->GetArray(i);
50 const char *name = array->GetName();
54 int numTuples = array->GetNumberOfTuples();
55 int numComponents = array->GetNumberOfComponents();
57 std::vector<T> values;
58 values.reserve(numTuples * numComponents);
60 for (
int t = 0; t < numTuples; ++t) {
62 array->GetTuple(t, tuple);
63 for (
int c = 0; c < numComponents; ++c)
64 values.push_back(tuple[c]);
67 metaData[name] = std::move(values);
77 : mesh(passedMesh), fileName(std::move(passedFileName)) {}
80 std::string passedFileName)
81 : mesh(passedMesh), fileFormat(passedFormat),
82 fileName(std::move(passedFileName)) {}
92 fileName = std::move(passedFileName);
99 if (mesh ==
nullptr) {
100 Logger::getInstance()
101 .addWarning(
"No mesh was passed to VTKReader. Not reading.")
106 if (fileName.empty()) {
107 Logger::getInstance()
108 .addWarning(
"No file name specified for VTKReader. Not reading.")
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.")
122 auto ending = fileName.substr(dotPos);
123 if (ending ==
".vtk") {
125 }
else if (ending ==
".vtp") {
127 }
else if (ending ==
".vtu") {
130 Logger::getInstance()
131 .addWarning(
"No valid file format found based on the file ending "
132 "passed to VTKReader. Not reading.")
139 switch (fileFormat) {
141 readVTKLegacy(fileName);
143#ifdef VIENNALS_USE_VTK
153 Logger::getInstance()
155 "VTKReader was built without VTK support. Only VTK_LEGACY "
156 "can be used. File not read.")
160 Logger::getInstance()
161 .addWarning(
"No valid file format set for VTKReader. Not reading.")
167#ifdef VIENNALS_USE_VTK
168 void readVTP(
const std::string &filename) {
171 vtkSmartPointer<vtkXMLPolyDataReader> pReader =
172 vtkSmartPointer<vtkXMLPolyDataReader>::New();
173 pReader->SetFileName(filename.c_str());
176 vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
177 polyData = pReader->GetOutput();
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;
186 vtkSmartPointer<vtkCellArray> cellArray =
187 vtkSmartPointer<vtkCellArray>::New();
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);
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);
212 mesh->lines.push_back(cell);
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);
227 mesh->triangles.push_back(cell);
232 vtkSmartPointer<vtkPointData> pointData =
233 vtkSmartPointer<vtkPointData>::New();
234 pointData = polyData->GetPointData();
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);
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]));
262 vtkSmartPointer<vtkCellData> cellData = vtkSmartPointer<vtkCellData>::New();
263 cellData = polyData->GetCellData();
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);
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]));
291 extractFieldData(polyData);
294 void readVTU(
const std::string &filename) {
298 vtkSmartPointer<vtkXMLUnstructuredGridReader> greader =
299 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
300 greader->SetFileName(filename.c_str());
303 vtkSmartPointer<vtkUnstructuredGrid> ugrid =
304 vtkSmartPointer<vtkUnstructuredGrid>::New();
305 ugrid = greader->GetOutput();
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;
316 for (
unsigned i = 0; i < ugrid->GetNumberOfCells(); ++i) {
317 vtkIdList *pointList = vtkIdList::New();
318 ugrid->GetCellPoints(i, pointList);
320 switch (ugrid->GetCellType(i)) {
323 std::array<unsigned, 1> vert{};
324 vert[0] = pointList->GetId(0);
325 mesh->vertices.push_back(vert);
329 std::array<unsigned, 2> elements{};
330 for (
unsigned j = 0; j < 2; ++j) {
331 elements[j] = pointList->GetId(j);
333 mesh->lines.push_back(elements);
337 std::array<unsigned, 3> elements{};
338 for (
unsigned j = 0; j < 3; ++j) {
339 elements[j] = pointList->GetId(j);
341 mesh->triangles.push_back(elements);
345 std::array<unsigned, 4> elements{};
346 for (
unsigned j = 0; j < 4; ++j) {
347 elements[j] = pointList->GetId(j);
349 mesh->tetras.push_back(elements);
353 std::array<unsigned, 8> elements{};
354 for (
unsigned j = 0; j < 8; ++j) {
355 elements[j] = pointList->GetId(j);
357 mesh->hexas.push_back(elements);
363 vtkSmartPointer<vtkPointData> pointData =
364 vtkSmartPointer<vtkPointData>::New();
365 pointData = ugrid->GetPointData();
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);
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]));
393 vtkSmartPointer<vtkCellData> cellData = vtkSmartPointer<vtkCellData>::New();
394 cellData = ugrid->GetCellData();
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);
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]));
422 extractFieldData(ugrid);
427 void readVTKLegacy(
const std::string &filename) {
431 std::ifstream f(filename.c_str());
433 Logger::getInstance().addError(
"Could not open geometry file!");
437 while (std::getline(f, temp)) {
438 if (temp.find(
"DATASET") != std::string::npos)
441 if (temp.find(
"UNSTRUCTURED_GRID") == std::string::npos) {
442 Logger::getInstance().addError(
"DATASET is not an UNSTRUCTURED_GRID!");
446 while (std::getline(f, temp)) {
447 if (temp.find(
"POINTS") != std::string::npos)
450 int num_nodes = atoi(&temp[temp.find(
' ') + 1]);
452 mesh->nodes.resize(num_nodes);
454 for (
int i = 0; i < num_nodes; i++) {
457 for (
double &coord : coords)
459 for (
int j = 0; j < 3; j++) {
460 mesh->nodes[i][j] = coords[j];
472 while (std::getline(f, temp)) {
473 if (temp.find(
"CELLS") == 0)
477 int num_elems = atoi(&temp[temp.find(
' ') + 1]);
479 std::ifstream f_ct(filename.c_str());
484 while (std::getline(f_ct, temp)) {
485 if (temp.find(
"CELL_TYPES") == 0)
488 int num_cell_types = atoi(&temp[temp.find(
' ') + 1]);
490 if (num_elems != num_cell_types) {
491 Logger::getInstance().addError(
492 "Corrupt input geometry! Number of CELLS and CELL_TYPES "
496 bool is_material =
true;
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);
513 std::vector<VectorType<unsigned int, 4>> elements;
514 elements.reserve(num_elems);
516 std::vector<double> materials;
520 unsigned cell_material;
521 for (
int i = 0; i < num_elems; i++) {
525 f_m >> cell_material;
531 unsigned number_nodes = vtk_nodes_for_cell_type[cell_type];
532 if (number_nodes == elems_fake || number_nodes == 0) {
536 std::array<unsigned, 1> elem{};
538 mesh->template getElements<1>().push_back(elem);
539 materials.push_back(cell_material);
543 std::array<unsigned, 2> elem{};
544 for (
unsigned j = 0; j < number_nodes; ++j) {
547 mesh->template getElements<2>().push_back(elem);
548 materials.push_back(cell_material);
553 std::array<unsigned, 3> elem{};
554 for (
unsigned j = 0; j < number_nodes; ++j) {
557 mesh->template getElements<3>().push_back(elem);
558 materials.push_back(cell_material);
564 std::array<unsigned, 4> elem{};
565 for (
unsigned j = 0; j < number_nodes; ++j) {
568 mesh->template getElements<4>().push_back(elem);
569 materials.push_back(cell_material);
575 std::array<unsigned, 3> elem{};
576 for (
unsigned j = 0; j < 3; ++j) {
579 mesh->template getElements<3>().push_back(
581 materials.push_back(cell_material);
584 mesh->template getElements<3>().push_back(elem);
585 materials.push_back(cell_material);
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();
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());
602 f.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
606 mesh->cellData.insertNextScalarData(materials,
"Material");
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]);
616 std::cout <<
"Read cell data: " << num_cell_data << std::endl;
619 std::cout <<
"reading scalar data" << std::endl;
620 while (std::getline(f, temp)) {
621 if (temp.find(
"SCALARS") != std::string::npos)
628 std::string scalarDataName;
630 auto firstS = temp.find(
' ') + 1;
631 auto secondS = temp.find(
' ', firstS + 1);
632 scalarDataName = temp.substr(firstS, secondS - firstS);
634 std::vector<double> scalarData;
637 std::getline(f, temp);
638 if (temp !=
"LOOKUP_TABLE default") {
639 Logger::getInstance()
640 .addWarning(
"Wrong lookup table for VTKLegacy: " + temp)
645 for (
int i = 0; i < num_cell_data; ++i) {
648 scalarData.push_back(data);
651 mesh->cellData.insertNextScalarData(scalarData, scalarDataName);