ViennaLS
Loading...
Searching...
No Matches
lsVTKWriter.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 <vtkCellArray.h>
15#include <vtkCellData.h>
16#include <vtkFloatArray.h>
17#include <vtkPointData.h>
18#include <vtkPoints.h>
19#include <vtkPolyData.h>
20#include <vtkSmartPointer.h>
21#include <vtkXMLPolyDataWriter.h>
22
23#include <vtkUnstructuredGrid.h>
24#include <vtkXMLUnstructuredGridWriter.h>
25#endif // VIENNALS_USE_VTK
26
27namespace viennals {
28
29using namespace viennacore;
30
32template <class T> class VTKWriter {
33 SmartPointer<Mesh<T>> mesh = nullptr;
35 std::string fileName;
36
37#ifdef VIENNALS_USE_VTK
38 template <class In, class Out>
39 void addDataFromMesh(const In &inData, Out outData) const {
40 // now add pointData
41 for (unsigned i = 0; i < inData.getScalarDataSize(); ++i) {
42 vtkSmartPointer<vtkFloatArray> pointData =
43 vtkSmartPointer<vtkFloatArray>::New();
44 pointData->SetNumberOfComponents(1);
45 pointData->SetName(inData.getScalarDataLabel(i).c_str());
46 auto scalars = *(inData.getScalarData(i));
47 for (unsigned j = 0; j < inData.getScalarData(i)->size(); ++j) {
48 pointData->InsertNextValue(scalars[j]);
49 }
50 outData->AddArray(pointData);
51 }
52
53 // now add vector data
54 for (unsigned i = 0; i < inData.getVectorDataSize(); ++i) {
55 vtkSmartPointer<vtkFloatArray> vectorData =
56 vtkSmartPointer<vtkFloatArray>::New();
57 vectorData->SetNumberOfComponents(3);
58 vectorData->SetName(inData.getVectorDataLabel(i).c_str());
59 auto vectors = *(inData.getVectorData(i));
60 for (unsigned j = 0; j < inData.getVectorData(i)->size(); ++j) {
61 vectorData->InsertNextTuple3(vectors[j][0], vectors[j][1],
62 vectors[j][2]);
63 }
64 outData->AddArray(vectorData);
65 }
66 }
67#endif // VIENNALS_USE_VTK
68
69public:
70 VTKWriter() = default;
71
72 VTKWriter(SmartPointer<Mesh<T>> passedMesh) : mesh(passedMesh) {}
73
74 VTKWriter(SmartPointer<Mesh<T>> passedMesh, std::string passedFileName)
75 : mesh(passedMesh), fileName(std::move(passedFileName)) {}
76
77 VTKWriter(SmartPointer<Mesh<T>> passedMesh, FileFormatEnum passedFormat,
78 std::string passedFileName)
79 : mesh(passedMesh), fileFormat(passedFormat),
80 fileName(std::move(passedFileName)) {}
81
82 void setMesh(SmartPointer<Mesh<T>> passedMesh) { mesh = passedMesh; }
83
85 void setFileFormat(FileFormatEnum passedFormat) { fileFormat = passedFormat; }
86
88 void setFileName(std::string passedFileName) {
89 fileName = std::move(passedFileName);
90 }
91
92 void apply() {
93 // check mesh
94 if (mesh == nullptr) {
95 Logger::getInstance()
96 .addWarning("No mesh was passed to VTKWriter. Not writing.")
97 .print();
98 return;
99 }
100 // check filename
101 if (fileName.empty()) {
102 Logger::getInstance()
103 .addWarning("No file name specified for VTKWriter. Not writing.")
104 .print();
105 return;
106 }
107
108 if (fileFormat == FileFormatEnum::VTK_AUTO) {
109 auto dotPos = fileName.rfind('.');
110 if (dotPos == std::string::npos) {
111 fileFormat = FileFormatEnum::VTP;
112 } else {
113 auto ending = fileName.substr(dotPos);
114 if (ending == ".vtk") {
115 fileFormat = FileFormatEnum::VTK_LEGACY;
116 } else if (ending == ".vtp") {
117 fileFormat = FileFormatEnum::VTP;
118 } else if (ending == ".vtu") {
119 fileFormat = FileFormatEnum::VTU;
120 } else {
121 Logger::getInstance()
122 .addWarning("No valid file format found based on the file ending "
123 "passed to VTKWriter. Not writing.")
124 .print();
125 return;
126 }
127 }
128 }
129
130 // check file format
131 switch (fileFormat) {
133 writeVTKLegacy(fileName);
134 break;
135#ifdef VIENNALS_USE_VTK
137 writeVTP(fileName);
138 break;
140 writeVTU(fileName);
141 break;
142#else
145 Logger::getInstance()
146 .addWarning("VTKWriter was built without VTK support. Falling back "
147 "to VTK_LEGACY.")
148 .print();
149 writeVTKLegacy(fileName);
150 break;
151#endif
152 default:
153 Logger::getInstance()
154 .addWarning("No valid file format set for VTKWriter. Not writing.")
155 .print();
156 }
157 }
158
159private:
160#ifdef VIENNALS_USE_VTK
161 void writeVTP(std::string filename) const {
162 if (mesh == nullptr) {
163 Logger::getInstance()
164 .addWarning("No mesh was passed to VTKWriter.")
165 .print();
166 return;
167 }
168
169 if (filename.find(".vtp") != filename.size() - 4)
170 filename += ".vtp";
171 vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
172
173 // Points
174 vtkSmartPointer<vtkPoints> polyPoints = vtkSmartPointer<vtkPoints>::New();
175 for (auto it = mesh->getNodes().begin(); it != mesh->getNodes().end();
176 ++it) {
177 polyPoints->InsertNextPoint((*it)[0], (*it)[1], (*it)[2]);
178 }
179 polyData->SetPoints(polyPoints);
180
181 // Vertices
182 if (mesh->vertices.size() > 0) {
183 vtkSmartPointer<vtkCellArray> polyCells =
184 vtkSmartPointer<vtkCellArray>::New();
185 for (auto it = mesh->vertices.begin(); it != mesh->vertices.end(); ++it) {
186 polyCells->InsertNextCell(1);
187 polyCells->InsertCellPoint((*it)[0]);
188 }
189 polyData->SetVerts(polyCells);
190 }
191
192 // Lines
193 if (mesh->lines.size() > 0) {
194 vtkSmartPointer<vtkCellArray> polyCells =
195 vtkSmartPointer<vtkCellArray>::New();
196 for (auto it = mesh->lines.begin(); it != mesh->lines.end(); ++it) {
197 polyCells->InsertNextCell(2);
198 for (unsigned i = 0; i < 2; ++i) {
199 polyCells->InsertCellPoint((*it)[i]);
200 }
201 }
202 polyData->SetLines(polyCells);
203 }
204
205 // Triangles
206 if (mesh->triangles.size() > 0) {
207 vtkSmartPointer<vtkCellArray> polyCells =
208 vtkSmartPointer<vtkCellArray>::New();
209 for (auto it = mesh->triangles.begin(); it != mesh->triangles.end();
210 ++it) {
211 polyCells->InsertNextCell(3);
212 for (unsigned i = 0; i < 3; ++i) {
213 polyCells->InsertCellPoint((*it)[i]);
214 }
215 }
216 polyData->SetPolys(polyCells);
217 }
218
219 addDataFromMesh(mesh->pointData, polyData->GetPointData());
220 addDataFromMesh(mesh->cellData, polyData->GetCellData());
221
222 vtkSmartPointer<vtkXMLPolyDataWriter> pwriter =
223 vtkSmartPointer<vtkXMLPolyDataWriter>::New();
224 pwriter->SetFileName(filename.c_str());
225 pwriter->SetInputData(polyData);
226 pwriter->Write();
227 }
228
229 void writeVTU(std::string filename) const {
230 if (mesh == nullptr) {
231 Logger::getInstance()
232 .addWarning("No mesh was passed to VTKWriter.")
233 .print();
234 return;
235 }
236
237 if (filename.find(".vtu") != filename.size() - 4)
238 filename += ".vtu";
239
240 vtkSmartPointer<vtkUnstructuredGrid> uGrid =
241 vtkSmartPointer<vtkUnstructuredGrid>::New();
242
243 // Points
244 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
245 for (auto it = mesh->getNodes().begin(); it != mesh->getNodes().end();
246 ++it) {
247 points->InsertNextPoint((*it)[0], (*it)[1], (*it)[2]);
248 }
249 uGrid->SetPoints(points);
250
251 // Now set all cells
252 vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
253 std::vector<int> cellTypes;
254 cellTypes.reserve(mesh->vertices.size() + mesh->lines.size() +
255 mesh->triangles.size() + mesh->tetras.size() +
256 mesh->hexas.size());
257
258 // Vertices
259 if (mesh->vertices.size() > 0) {
260 for (auto it = mesh->vertices.begin(); it != mesh->vertices.end(); ++it) {
261 cells->InsertNextCell(1);
262 cells->InsertCellPoint((*it)[0]);
263 cellTypes.push_back(1); // vtk Vertex
264 }
265 }
266
267 // Lines
268 if (mesh->lines.size() > 0) {
269 for (auto it = mesh->lines.begin(); it != mesh->lines.end(); ++it) {
270 cells->InsertNextCell(2);
271 for (unsigned i = 0; i < 2; ++i) {
272 cells->InsertCellPoint((*it)[i]);
273 }
274 cellTypes.push_back(3); // vtk Line
275 }
276 }
277
278 // Triangles
279 if (mesh->triangles.size() > 0) {
280 for (auto it = mesh->triangles.begin(); it != mesh->triangles.end();
281 ++it) {
282 cells->InsertNextCell(3);
283 for (unsigned i = 0; i < 3; ++i) {
284 cells->InsertCellPoint((*it)[i]);
285 }
286 cellTypes.push_back(5); // vtk Triangle
287 }
288 }
289
290 // Tetras
291 if (mesh->tetras.size() > 0) {
292 for (auto it = mesh->tetras.begin(); it != mesh->tetras.end(); ++it) {
293 cells->InsertNextCell(4);
294 for (unsigned i = 0; i < 4; ++i) {
295 cells->InsertCellPoint((*it)[i]);
296 }
297 cellTypes.push_back(10); // vtk Tetra
298 }
299 }
300
301 // Hexas
302 if (mesh->hexas.size() > 0) {
303 for (auto it = mesh->hexas.begin(); it != mesh->hexas.end(); ++it) {
304 cells->InsertNextCell(8);
305 for (unsigned i = 0; i < 8; ++i) {
306 cells->InsertCellPoint((*it)[i]);
307 }
308 cellTypes.push_back(12); // vtk Hexahedron
309 }
310 }
311
312 // set cells
313 uGrid->SetCells(&(cellTypes[0]), cells);
314
315 addDataFromMesh(mesh->pointData, uGrid->GetPointData());
316 addDataFromMesh(mesh->cellData, uGrid->GetCellData());
317
318 // // now add pointData
319 // for (unsigned i = 0; i < mesh->cellData.getScalarDataSize(); ++i) {
320 // vtkSmartPointer<vtkFloatArray> pointData =
321 // vtkSmartPointer<vtkFloatArray>::New();
322 // pointData->SetNumberOfComponents(1);
323 // pointData->SetName(mesh->cellData.getScalarDataLabel(i).c_str());
324 // auto scalars = *(mesh->cellData.getScalarData(i));
325 // for (unsigned j = 0; j < scalars.size(); ++j) {
326 // pointData->InsertNextValue(scalars[j]);
327 // }
328 // uGrid->GetCellData()->AddArray(pointData);
329 // }
330
331 // // now add vector data
332 // for (unsigned i = 0; i < mesh->cellData.getVectorDataSize(); ++i) {
333 // vtkSmartPointer<vtkFloatArray> vectorData =
334 // vtkSmartPointer<vtkFloatArray>::New();
335 // vectorData->SetNumberOfComponents(3);
336 // vectorData->SetName(mesh->cellData.getVectorDataLabel(i).c_str());
337 // auto vectors = *(mesh->cellData.getVectorData(i));
338 // for (unsigned j = 0; j < vectors.size(); ++j) {
339 // vectorData->InsertNextTuple3(vectors[j][0], vectors[j][1],
340 // vectors[j][2]);
341 // }
342 // uGrid->GetCellData()->AddArray(vectorData);
343 // }
344
345 vtkSmartPointer<vtkXMLUnstructuredGridWriter> owriter =
346 vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
347 owriter->SetFileName(filename.c_str());
348 owriter->SetInputData(uGrid);
349 owriter->Write();
350 }
351
352#endif // VIENNALS_USE_VTK
353
354 void writeVTKLegacy(const std::string &filename) {
355 if (mesh == nullptr) {
356 Logger::getInstance()
357 .addWarning("No mesh was passed to VTKWriter.")
358 .print();
359 return;
360 }
361
362 std::ofstream f(filename.c_str());
363
364 f << "# vtk DataFile Version 2.0" << std::endl;
365 f << ((!mesh->lines.empty()) ? 2 : 3) << "D Surface" << std::endl;
366 f << "ASCII" << std::endl;
367 f << "DATASET UNSTRUCTURED_GRID" << std::endl;
368 f << "POINTS " << mesh->nodes.size() << " float" << std::endl;
369
370 // print nodes
371 for (unsigned int i = 0; i < mesh->nodes.size(); i++) {
372 for (int j = 0; j < 3; j++)
373 f << static_cast<float>(mesh->nodes[i][j]) << " ";
374 f << std::endl;
375 }
376
377 const unsigned numberOfCells = mesh->vertices.size() + mesh->lines.size() +
378 mesh->triangles.size() +
379 mesh->tetras.size() + mesh->hexas.size();
380 const unsigned cellDataSize =
381 2 * mesh->vertices.size() + 3 * mesh->lines.size() +
382 4 * mesh->triangles.size() + 5 * mesh->tetras.size() +
383 9 * mesh->hexas.size();
384
385 f << "CELLS " << numberOfCells << " " << cellDataSize << std::endl;
386
387 // print elements
388 for (unsigned int i = 0; i < mesh->vertices.size(); i++) {
389 f << 1 << " " << mesh->vertices[i][0] << std::endl;
390 }
391 for (unsigned int i = 0; i < mesh->lines.size(); i++) {
392 f << 2 << " ";
393 for (int j = 0; j < 2; j++)
394 f << mesh->lines[i][j] << " ";
395 f << std::endl;
396 }
397
398 for (unsigned int i = 0; i < mesh->triangles.size(); i++) {
399 f << 3 << " ";
400 for (int j = 0; j < 3; j++)
401 f << mesh->triangles[i][j] << " ";
402 f << std::endl;
403 }
404
405 for (unsigned int i = 0; i < mesh->tetras.size(); i++) {
406 f << 4 << " ";
407 for (int j = 0; j < 4; j++)
408 f << mesh->tetras[i][j] << " ";
409 f << std::endl;
410 }
411
412 for (unsigned int i = 0; i < mesh->hexas.size(); i++) {
413 f << 8 << " ";
414 for (int j = 0; j < 8; j++)
415 f << mesh->hexas[i][j] << " ";
416 f << std::endl;
417 }
418
419 f << "CELL_TYPES " << numberOfCells << std::endl;
420 for (unsigned i = 0; i < mesh->vertices.size(); ++i)
421 f << 1 << std::endl;
422
423 for (unsigned i = 0; i < mesh->lines.size(); ++i)
424 f << 3 << std::endl;
425
426 for (unsigned i = 0; i < mesh->triangles.size(); ++i)
427 f << 5 << std::endl;
428
429 for (unsigned i = 0; i < mesh->tetras.size(); ++i)
430 f << 10 << std::endl;
431
432 for (unsigned i = 0; i < mesh->hexas.size(); ++i)
433 f << 12 << std::endl;
434
435 // WRITE POINT DATA
436 if (mesh->pointData.getScalarDataSize() ||
437 mesh->pointData.getVectorDataSize()) {
438 Logger::getInstance()
439 .addWarning("Point data output not supported for legacy VTK output. "
440 "Point data is ignored.")
441 .print();
442 }
443
444 // WRITE SCALAR DATA
445 if (mesh->cellData.getScalarDataSize()) {
446 f << "CELL_DATA " << mesh->cellData.getScalarData(0)->size() << std::endl;
447 for (unsigned i = 0; i < mesh->cellData.getScalarDataSize(); ++i) {
448 auto scalars = *(mesh->cellData.getScalarData(i));
449 f << "SCALARS " << mesh->cellData.getScalarDataLabel(i) << " float"
450 << std::endl;
451 f << "LOOKUP_TABLE default" << std::endl;
452 for (unsigned j = 0; j < scalars.size(); ++j) {
453 f << ((std::abs(scalars[j]) < 1e-6) ? 0.0 : scalars[j]) << std::endl;
454 }
455 }
456 }
457
458 // WRITE VECTOR DATA
459 if (mesh->cellData.getVectorDataSize()) {
460 if (!mesh->cellData.getScalarDataSize())
461 f << "CELL_DATA " << mesh->cellData.getVectorData(0)->size()
462 << std::endl;
463 for (unsigned i = 0; i < mesh->cellData.getVectorDataSize(); ++i) {
464 auto vectors = *(mesh->cellData.getVectorData(i));
465 f << "VECTORS " << mesh->cellData.getVectorDataLabel(i) << " float"
466 << std::endl;
467 for (unsigned j = 0; j < vectors.size(); ++j) {
468 for (unsigned k = 0; k < 3; ++k) {
469 f << vectors[j][k] << " ";
470 }
471 f << std::endl;
472 }
473 }
474 }
475
476 f.close();
477 }
478};
479
480} // 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
void setFileName(std::string passedFileName)
set file name for file to write
Definition lsVTKWriter.hpp:88
VTKWriter(SmartPointer< Mesh< T > > passedMesh, FileFormatEnum passedFormat, std::string passedFileName)
Definition lsVTKWriter.hpp:77
VTKWriter(SmartPointer< Mesh< T > > passedMesh, std::string passedFileName)
Definition lsVTKWriter.hpp:74
void apply()
Definition lsVTKWriter.hpp:92
void setMesh(SmartPointer< Mesh< T > > passedMesh)
Definition lsVTKWriter.hpp:82
void setFileFormat(FileFormatEnum passedFormat)
set file format for file to write. Defaults to VTK_LEGACY.
Definition lsVTKWriter.hpp:85
VTKWriter(SmartPointer< Mesh< T > > passedMesh)
Definition lsVTKWriter.hpp:72
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