ViennaLS
Loading...
Searching...
No Matches
lsWriteVisualizationMesh.hpp
Go to the documentation of this file.
1#pragma once
2
3#ifdef VIENNALS_USE_VTK // this class needs vtk support
4
5#include <vtkAppendFilter.h>
6#include <vtkAppendPolyData.h>
7#include <vtkCellData.h>
8#include <vtkDataSetTriangleFilter.h>
9#include <vtkFloatArray.h>
10#include <vtkGeometryFilter.h>
11#include <vtkIncrementalOctreePointLocator.h>
12#include <vtkIntArray.h>
13#include <vtkPointData.h>
14#include <vtkProbeFilter.h>
15#include <vtkRectilinearGrid.h>
16#include <vtkSmartPointer.h>
17#include <vtkTableBasedClipDataSet.h>
18#include <vtkTriangleFilter.h>
19#include <vtkUnstructuredGrid.h>
20#include <vtkXMLPolyDataWriter.h>
21#include <vtkXMLUnstructuredGridWriter.h>
22
23#include <hrleDenseIterator.hpp>
24
25#include <lsDomain.hpp>
26#include <lsMaterialMap.hpp>
28#include <unordered_map>
29#include <utility>
30
31// #define LS_TO_VISUALIZATION_DEBUG
32#ifdef LS_TO_VISUALIZATION_DEBUG
33#include <vtkXMLRectilinearGridWriter.h>
34#endif
35
36namespace viennals {
37
38using namespace viennacore;
39
46template <class T, int D> class WriteVisualizationMesh {
47 typedef typename Domain<T, D>::DomainType hrleDomainType;
48 using LevelSetsType = std::vector<SmartPointer<Domain<T, D>>>;
49 LevelSetsType levelSets;
50 SmartPointer<MaterialMap> materialMap = nullptr;
51 std::string fileName;
52 bool extractVolumeMesh = true;
53 bool extractHullMesh = false;
54 bool bottomRemoved = false;
55 double LSEpsilon = 1e-2;
56 std::unordered_map<std::string, std::vector<T>> metaData;
57
61 static void removeDuplicatePoints(vtkSmartPointer<vtkPolyData> &polyData,
62 const double tolerance) {
63
64 vtkSmartPointer<vtkPolyData> newPolyData =
65 vtkSmartPointer<vtkPolyData>::New();
66 vtkSmartPointer<vtkIncrementalOctreePointLocator> ptInserter =
67 vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
68 ptInserter->SetTolerance(tolerance);
69
70 vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
71
72 // get bounds
73 double gridBounds[6];
74 polyData->GetBounds(gridBounds);
75
76 // start point insertion from original points
77 std::vector<vtkIdType> newPointIds;
78 newPointIds.reserve(polyData->GetNumberOfPoints());
79 ptInserter->InitPointInsertion(newPoints, gridBounds);
80
81 // make new point list
82 for (vtkIdType pointId = 0; pointId < polyData->GetNumberOfPoints();
83 ++pointId) {
84 vtkIdType globalPtId = 0;
85 ptInserter->InsertUniquePoint(polyData->GetPoint(pointId), globalPtId);
86 newPointIds.push_back(globalPtId);
87 }
88
89 // now set the new points to the unstructured grid
90 newPolyData->SetPoints(newPoints);
91
92 // go through all cells and change point ids to match the new ids
93 vtkSmartPointer<vtkCellArray> oldCells = polyData->GetPolys();
94 vtkSmartPointer<vtkCellArray> newCells =
95 vtkSmartPointer<vtkCellArray>::New();
96
97 vtkSmartPointer<vtkIdList> cellPoints = vtkIdList::New();
98 oldCells->InitTraversal();
99 while (oldCells->GetNextCell(cellPoints)) {
100 for (vtkIdType pointId = 0; pointId < cellPoints->GetNumberOfIds();
101 ++pointId) {
102 cellPoints->SetId(pointId, newPointIds[cellPoints->GetId(pointId)]);
103 }
104 // insert same cell with new points
105 newCells->InsertNextCell(cellPoints);
106 }
107
108 newPolyData->SetPolys(newCells);
109
110 // conserve all cell data
111 // TODO transfer point data as well (do with "InsertTuples (vtkIdList
112 // *dstIds, vtkIdList *srcIds, vtkAbstractArray *source) override)" of
113 // vtkDataArray class)
114 newPolyData->GetCellData()->ShallowCopy(polyData->GetCellData());
115
116 // set ugrid to the newly created grid
117 polyData = newPolyData;
118 }
119
122 static void removeDuplicatePoints(vtkSmartPointer<vtkUnstructuredGrid> &ugrid,
123 const double tolerance) {
124
125 vtkSmartPointer<vtkUnstructuredGrid> newGrid =
126 vtkSmartPointer<vtkUnstructuredGrid>::New();
127 vtkSmartPointer<vtkIncrementalOctreePointLocator> ptInserter =
128 vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
129 ptInserter->SetTolerance(tolerance);
130
131 vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
132
133 // get bounds
134 double gridBounds[6];
135 ugrid->GetBounds(gridBounds);
136
137 // start point insertion from original points
138 std::vector<vtkIdType> newPointIds;
139 newPointIds.reserve(ugrid->GetNumberOfPoints());
140 ptInserter->InitPointInsertion(newPoints, gridBounds);
141
142 // make new point list
143 for (vtkIdType pointId = 0; pointId < ugrid->GetNumberOfPoints();
144 ++pointId) {
145 vtkIdType globalPtId = 0;
146 ptInserter->InsertUniquePoint(ugrid->GetPoint(pointId), globalPtId);
147 newPointIds.push_back(globalPtId);
148 }
149
150 // now set the new points to the unstructured grid
151 newGrid->SetPoints(newPoints);
152
153 // go through all cells and change point ids to match the new ids
154 for (vtkIdType cellId = 0; cellId < ugrid->GetNumberOfCells(); ++cellId) {
155 vtkSmartPointer<vtkIdList> cellPoints = vtkSmartPointer<vtkIdList>::New();
156 ugrid->GetCellPoints(cellId, cellPoints);
157 for (vtkIdType pointId = 0; pointId < cellPoints->GetNumberOfIds();
158 ++pointId) {
159 cellPoints->SetId(pointId, newPointIds[cellPoints->GetId(pointId)]);
160 }
161 // insert same cell with new points
162 newGrid->InsertNextCell(ugrid->GetCell(cellId)->GetCellType(),
163 cellPoints);
164 }
165
166 // conserve all cell data
167 // TODO transfer point data as well (do with "InsertTuples (vtkIdList
168 // *dstIds, vtkIdList *srcIds, vtkAbstractArray *source) override)" of
169 // vtkDataArray class)
170 newGrid->GetCellData()->ShallowCopy(ugrid->GetCellData());
171
172 // set ugrid to the newly created grid
173 ugrid = newGrid;
174 }
175
177 static void
178 removeDegenerateTetras(vtkSmartPointer<vtkUnstructuredGrid> &ugrid) {
179 vtkSmartPointer<vtkUnstructuredGrid> newGrid =
180 vtkSmartPointer<vtkUnstructuredGrid>::New();
181
182 // need to copy material numbers
183 vtkSmartPointer<vtkIntArray> materialNumberArray =
184 vtkSmartPointer<vtkIntArray>::New();
185 materialNumberArray->SetNumberOfComponents(1);
186 materialNumberArray->SetName("Material");
187
188 // see if material is defined
189 int arrayIndex;
190 vtkDataArray *matArray =
191 ugrid->GetCellData()->GetArray("Material", arrayIndex);
192 const int &materialArrayIndex = arrayIndex;
193
194 // go through all cells and delete those with duplicate entries
195 for (vtkIdType cellId = 0; cellId < ugrid->GetNumberOfCells(); ++cellId) {
196 vtkSmartPointer<vtkIdList> cellPoints = vtkSmartPointer<vtkIdList>::New();
197 ugrid->GetCellPoints(cellId, cellPoints);
198 bool isDuplicate = false;
199 for (vtkIdType pointId = 0; pointId < cellPoints->GetNumberOfIds();
200 ++pointId) {
201 for (vtkIdType nextId = pointId + 1;
202 nextId < cellPoints->GetNumberOfIds(); ++nextId) {
203 // if they are the same, remove the cell
204 if (cellPoints->GetId(pointId) == cellPoints->GetId(nextId))
205 isDuplicate = true;
206 }
207 }
208 if (!isDuplicate) {
209 // insert same cell if no degenerate points
210 newGrid->InsertNextCell(ugrid->GetCell(cellId)->GetCellType(),
211 cellPoints);
212 // if material was defined before, use it now
213 if (materialArrayIndex >= 0)
214 materialNumberArray->InsertNextValue(matArray->GetTuple1(cellId));
215 }
216 }
217
218 // just take the old points and point data
219 newGrid->SetPoints(ugrid->GetPoints());
220 newGrid->GetPointData()->ShallowCopy(ugrid->GetPointData());
221 // set material cell data
222 newGrid->GetCellData()->SetScalars(materialNumberArray);
223
224 ugrid = newGrid;
225 }
226
227 // This function takes a levelset and converts it to a vtkRectilinearGrid
228 // The full domain contains values, which are capped at numLayers * gridDelta
229 // gridExtraPoints layers of grid points are added to the domain according to
230 // boundary conditions
231 template <int gridExtraPoints = 0>
232 vtkSmartPointer<vtkRectilinearGrid>
233 LS2RectiLinearGrid(SmartPointer<Domain<T, D>> levelSet, const double LSOffset,
234 int infiniteMinimum = std::numeric_limits<int>::max(),
235 int infiniteMaximum = -std::numeric_limits<int>::max()) {
236
237 auto &grid = levelSet->getGrid();
238 auto &domain = levelSet->getDomain();
239 double gridDelta = grid.getGridDelta();
240 int numLayers = levelSet->getLevelSetWidth();
241
242 vtkSmartPointer<vtkFloatArray>
243 coords[3]; // needs to be 3 because vtk only knows 3D
244 int gridMin = 0, gridMax = 0;
245 // int openJumpDirection = -1; // direction above the open boundary
246 // direction
247
248 // fill grid with offset depending on orientation
249 for (unsigned i = 0; i < D; ++i) {
250 coords[i] = vtkSmartPointer<vtkFloatArray>::New();
251
252 if (grid.getBoundaryConditions(i) ==
253 Domain<T, D>::BoundaryType::INFINITE_BOUNDARY) {
254 // add one to gridMin and gridMax for numerical stability
255 gridMin = std::min(domain.getMinRunBreak(i), infiniteMinimum) -
256 1; // choose the smaller number so that for first levelset the
257 // overall minimum can be chosen
258 gridMax = std::max(domain.getMaxRunBreak(i), infiniteMaximum) + 1;
259
260 // openJumpDirection = i + 1;
261 } else {
262 gridMin = grid.getMinGridPoint(i) - gridExtraPoints;
263 gridMax = grid.getMaxGridPoint(i) + gridExtraPoints;
264 }
265
266 for (int x = gridMin; x <= gridMax; ++x) {
267 coords[i]->InsertNextValue(x * gridDelta);
268 }
269 }
270
271 // if we work in 2D, just add 1 grid point at origin
272 if (D == 2) {
273 coords[2] = vtkSmartPointer<vtkFloatArray>::New();
274 coords[2]->InsertNextValue(0);
275 }
276
277 vtkSmartPointer<vtkRectilinearGrid> rgrid =
278 vtkSmartPointer<vtkRectilinearGrid>::New();
279
280 rgrid->SetDimensions(coords[0]->GetNumberOfTuples(),
281 coords[1]->GetNumberOfTuples(),
282 coords[2]->GetNumberOfTuples());
283 rgrid->SetXCoordinates(coords[0]);
284 rgrid->SetYCoordinates(coords[1]);
285 rgrid->SetZCoordinates(coords[2]);
286
287 // bool fixBorderPoints = (gridExtraPoints != 0);
288
289 // need to save the current position one dimension above open boundary
290 // direction, so we can register a jump in the open boundary direction when
291 // it occurs, so we can fix the LS value as follows: if remove_bottom==true,
292 // we need to flip the sign otherwise the sign stays the same
293 // int currentOpenIndex =
294 // it.getIndices()[(openJumpDirection < D) ? openJumpDirection : 0];
295
296 // Make array to store signed distance function
297 auto const numGridPoints = rgrid->GetNumberOfPoints();
298 vtkSmartPointer<vtkFloatArray> signedDistances =
299 vtkSmartPointer<vtkFloatArray>::New();
300 signedDistances->SetNumberOfComponents(1);
301 signedDistances->SetNumberOfTuples(numGridPoints);
302 signedDistances->SetName("SignedDistances");
303
304#pragma omp parallel
305 {
306 // use dense iterator to got to every index location
307 viennahrle::ConstDenseIterator<typename Domain<T, D>::DomainType> it(
308 levelSet->getDomain());
309
310#pragma omp for
311 for (vtkIdType pointId = 0; pointId < numGridPoints; ++pointId) {
312 // iterate until all grid points have a signed distance value
313
314 double p[3];
315 rgrid->GetPoint(pointId, p);
316 // create index vector
317 viennahrle::Index<D> indices(grid.globalCoordinates2GlobalIndices(p));
318
319 // write the corresponding LSValue
320 T value;
321
322 // if indices are outside of domain map to local indices
323 if (grid.isOutsideOfDomain(indices)) {
324 indices = grid.globalIndices2LocalIndices(indices);
325 // fixBorderPoints = true;
326 // signedDistances->InsertNextValue(
327 // signedDistances->GetDataTypeValueMax());
328 }
329
330 it.goToIndices(indices);
331 if (it.getValue() == Domain<T, D>::POS_VALUE) {
332 value = numLayers;
333 } else if (it.getValue() == Domain<T, D>::NEG_VALUE) {
334 value = -numLayers;
335 } else {
336 value = it.getValue() + LSOffset;
337 }
338
339 // if (removeBottom) {
340 // // if we jump from one end of the domain to the other and are not
341 // // already in the new run, we need to fix the sign of the run
342 // if (currentOpenIndex != indices[openJumpDirection]) {
343 // currentOpenIndex = indices[openJumpDirection];
344 // if (indices >= it.getIndices()) {
345 // value = -value;
346 // }
347 // }
348 // }
349
350 signedDistances->SetValue(pointId, value * gridDelta);
351 // signedDistances->InsertNextValue(value * gridDelta);
352 }
353
354 // // advance iterator until it is at correct point
355 // while (compare(it_l.end_indices(), indices) < 0) {
356 // it_l.next();
357 // if (it_l.is_finished())
358 // break;
359 // }
360 // // now move iterator with pointId to get to next point
361 // switch (compare(it_l.end_indices(), indices)) {
362 // case 0:
363 // it_l.next();
364 // default:
365 // ++pointId;
366 // }
367 // }
368 }
369
370 // now need to go through again to fix border points, this is done by
371 // mapping existing points onto the points outside the domain according
372 // to the correct boundary conditions
373 // if (fixBorderPoints) {
374 // vtkIdType pointId = 0;
375 // while ((pointId < rgrid->GetNumberOfPoints())) {
376 // if (signedDistances->GetValue(pointId) ==
377 // signedDistances->GetDataTypeValueMax()) {
378 // double p[3];
379 // rgrid->GetPoint(pointId, p);
380
381 // // create index vector
382 // viennahrle::Index<D> indices(
383 // grid.globalCoordinates2GlobalIndices(p));
384
385 // // vector for mapped point inside domain
386 // viennahrle::Index<D> localIndices =
387 // grid.globalIndices2LocalIndices(indices);
388
389 // // now find ID of point we need to take value from
390 // int originalPointId = 0;
391 // for (int i = D - 1; i >= 0; --i) {
392 // originalPointId *=
393 // coords[i]->GetNumberOfTuples(); // extent in direction
394 // originalPointId += localIndices[i] - indices[i];
395 // }
396 // originalPointId += pointId;
397
398 // // now put value of mapped point in global point
399 // signedDistances->SetValue(pointId,
400 // signedDistances->GetValue(originalPointId));
401 // }
402 // ++pointId;
403 // }
404 // }
405
406 // Add the SignedDistances to the grid
407 rgrid->GetPointData()->SetScalars(signedDistances);
408
409 return rgrid;
410 }
411
412 void addMetaDataToVTK(vtkDataSet *data) const {
413 if (metaData.empty()) {
414 return;
415 }
416
417 // add metadata to field data
418 vtkSmartPointer<vtkFieldData> fieldData = data->GetFieldData();
419 for (const auto &meta : metaData) {
420 if (meta.second.empty())
421 continue; // skip empty metadata
422
423 vtkSmartPointer<vtkFloatArray> metaDataArray =
424 vtkSmartPointer<vtkFloatArray>::New();
425 metaDataArray->SetName(meta.first.c_str());
426 metaDataArray->SetNumberOfValues(meta.second.size());
427 for (size_t i = 0; i < meta.second.size(); ++i) {
428 metaDataArray->SetValue(i, meta.second[i]);
429 }
430 fieldData->AddArray(metaDataArray);
431 }
432 }
433
434public:
435 WriteVisualizationMesh() = default;
436
437 WriteVisualizationMesh(SmartPointer<Domain<T, D>> levelSet) {
438 levelSets.push_back(levelSet);
439 }
440
442 void insertNextLevelSet(SmartPointer<Domain<T, D>> levelSet) {
443 levelSets.push_back(levelSet);
444 }
445
448 void setFileName(std::string passedFileName) {
449 fileName = std::move(passedFileName);
450 }
451
453 void setExtractHullMesh(bool passedExtractHullMesh) {
454 extractHullMesh = passedExtractHullMesh;
455 }
456
458 void setExtractVolumeMesh(bool passedExtractVolumeMesh) {
459 extractVolumeMesh = passedExtractVolumeMesh;
460 }
461
462 void setMaterialMap(SmartPointer<MaterialMap> passedMaterialMap) {
463 materialMap = passedMaterialMap;
464 }
465
466 void setWrappingLayerEpsilon(double epsilon) { LSEpsilon = epsilon; }
467
468 void setMetaData(
469 const std::unordered_map<std::string, std::vector<T>> &passedMetaData) {
470 metaData = passedMetaData;
471 }
472
473 void addMetaData(const std::string &key, T value) {
474 metaData[key] = std::vector<T>{value};
475 }
476
477 void addMetaData(const std::string &key, const std::vector<T> &values) {
478 metaData[key] = values;
479 }
480
481 void addMetaData(
482 const std::unordered_map<std::string, std::vector<T>> &newMetaData) {
483 for (const auto &pair : newMetaData) {
484 metaData[pair.first] = pair.second;
485 }
486 }
487
488 void apply() {
489 // check if level sets have enough layers
490 for (unsigned i = 0; i < levelSets.size(); ++i) {
491 if (levelSets[i]->getLevelSetWidth() < 2) {
492 Logger::getInstance()
493 .addWarning(
494 "WriteVisualizationMesh: Level Set " + std::to_string(i) +
495 " should have a width greater than 1! Conversion might fail!")
496 .print();
497 }
498 }
499
500 const double gridDelta = levelSets[0]->getGrid().getGridDelta();
501
502 // store volume for each material
503 std::vector<vtkSmartPointer<vtkUnstructuredGrid>> materialMeshes;
504 std::vector<unsigned> materialIds;
505
506 int totalMinimum = std::numeric_limits<int>::max();
507 int totalMaximum = -std::numeric_limits<int>::max();
508 for (auto &it : levelSets) {
509 if (it->getNumberOfPoints() == 0) {
510 continue;
511 }
512 auto &grid = it->getGrid();
513 auto &domain = it->getDomain();
514 for (unsigned i = 0; i < D; ++i) {
515 if (grid.getBoundaryConditions(i) ==
516 Domain<T, D>::BoundaryType::INFINITE_BOUNDARY) {
517 totalMinimum = std::min(totalMinimum, domain.getMinRunBreak(i));
518 totalMaximum = std::max(totalMaximum, domain.getMaxRunBreak(i));
519 }
520 }
521 }
522
523 // create volume mesh for largest LS
524 // Use vtkClipDataSet to slice the grid
525 vtkSmartPointer<vtkTableBasedClipDataSet> clipper =
526 vtkSmartPointer<vtkTableBasedClipDataSet>::New();
527 auto topGrid = vtkSmartPointer<vtkRectilinearGrid>::New();
528 // if (bottomRemoved) {
529 topGrid =
530 LS2RectiLinearGrid(levelSets.back(), 0, totalMinimum, totalMaximum);
531 // } else {
532 // topGrid = LS2RectiLinearGrid<false>(levelSets.back(), 0,
533 // totalMinimum,
534 // totalMaximum);
535 // }
536#ifdef LS_TO_VISUALIZATION_DEBUG
537 {
538 auto gwriter = vtkSmartPointer<vtkXMLRectilinearGridWriter>::New();
539 gwriter->SetFileName("./grid_0.vtr");
540 gwriter->SetInputData(topGrid);
541 gwriter->Write();
542 std::cout << "Wrote grid 0" << std::endl;
543 }
544#endif
545 clipper->SetInputData(topGrid);
546 clipper->InsideOutOn();
547 clipper->SetValue(0.0);
548 clipper->GenerateClippedOutputOn(); // TODO remove
549 clipper->Update();
550
551#ifdef LS_TO_VISUALIZATION_DEBUG
552 {
553 auto gwriter = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
554 gwriter->SetFileName("./clipped.vtu");
555 gwriter->SetInputData(clipper->GetClippedOutput());
556 gwriter->Write();
557 std::cout << "Wrote clipped" << std::endl;
558 }
559#endif
560
561 const bool useMaterialMap = materialMap != nullptr;
562 materialMeshes.emplace_back(clipper->GetOutput());
563 materialIds.push_back(useMaterialMap ? materialMap->getMaterialId(0) : 0);
564
565#ifdef LS_TO_VISUALIZATION_DEBUG
566 {
567 auto gwriter = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
568 gwriter->SetFileName("./probed_0.vtu");
569 gwriter->SetInputData(materialMeshes.front());
570 gwriter->Write();
571 }
572#endif
573
574 unsigned counter = 1;
575
576 // now cut large volume mesh with all the smaller ones
577 for (typename LevelSetsType::const_reverse_iterator it =
578 ++levelSets.rbegin();
579 it != levelSets.rend(); ++it) {
580 if (it->get()->getNumberOfPoints() == 0)
581 continue; // ignore empty levelSets
582
583 // create grid of next LS with slight offset and project into current
584 // mesh
585 vtkSmartPointer<vtkRectilinearGrid> rgrid =
586 vtkSmartPointer<vtkRectilinearGrid>::New();
587 // if (bottomRemoved) {
588 rgrid = LS2RectiLinearGrid<1>(*it, -LSEpsilon * counter, totalMinimum,
589 totalMaximum);
590 // } else {
591 // rgrid = LS2RectiLinearGrid<false, 1>(*it, -LSEpsilon * counter,
592 // totalMinimum, totalMaximum);
593 // }
594
595#ifdef LS_TO_VISUALIZATION_DEBUG
596 {
597 vtkSmartPointer<vtkXMLRectilinearGridWriter> gwriter =
598 vtkSmartPointer<vtkXMLRectilinearGridWriter>::New();
599 gwriter->SetFileName(
600 ("./grid_" + std::to_string(counter) + ".vtr").c_str());
601 gwriter->SetInputData(rgrid);
602 gwriter->Write();
603 std::cout << "Wrote grid " << counter << std::endl;
604 }
605#endif
606
607 // now transfer implicit values to mesh points
608 vtkSmartPointer<vtkProbeFilter> probeFilter =
609 vtkSmartPointer<vtkProbeFilter>::New();
610 probeFilter->SetInputData(materialMeshes.back()); // last element
611 probeFilter->SetSourceData(rgrid);
612 probeFilter->Update();
613
614#ifdef LS_TO_VISUALIZATION_DEBUG
615 {
616 vtkSmartPointer<vtkXMLUnstructuredGridWriter> gwriter =
617 vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
618 gwriter->SetFileName(
619 ("./probed_" + std::to_string(counter) + ".vtu").c_str());
620 gwriter->SetInputData(probeFilter->GetOutput());
621 gwriter->Write();
622 std::cout << "Wrote unstructured grid " << counter << std::endl;
623 }
624#endif
625
626 // now clip the mesh and save the clipped as the 1st layer and use the
627 // inverse for the next layer clipping Use vtkTabelBasedClipDataSet to
628 // slice the grid
629 vtkSmartPointer<vtkTableBasedClipDataSet> insideClipper =
630 vtkSmartPointer<vtkTableBasedClipDataSet>::New();
631 insideClipper->SetInputConnection(probeFilter->GetOutputPort());
632 insideClipper->GenerateClippedOutputOn();
633 insideClipper->Update();
634
635 materialMeshes.back() = insideClipper->GetOutput();
636 materialMeshes.emplace_back(insideClipper->GetClippedOutput());
637 auto material = counter;
638 if (useMaterialMap)
639 material = materialMap->getMaterialId(counter);
640 materialIds.push_back(material);
641
642 ++counter;
643 }
644
645 vtkSmartPointer<vtkAppendFilter> appendFilter =
646 vtkSmartPointer<vtkAppendFilter>::New();
647
648 vtkSmartPointer<vtkAppendPolyData> hullAppendFilter =
649 vtkSmartPointer<vtkAppendPolyData>::New();
650
651 for (unsigned i = 0; i < materialMeshes.size(); ++i) {
652
653 // write material number in mesh
654 vtkSmartPointer<vtkIntArray> materialNumberArray =
655 vtkSmartPointer<vtkIntArray>::New();
656 materialNumberArray->SetNumberOfComponents(1);
657 materialNumberArray->SetName("Material");
658 for (unsigned j = 0;
659 j <
660 materialMeshes[materialMeshes.size() - 1 - i]->GetNumberOfCells();
661 ++j) {
662 materialNumberArray->InsertNextValue(materialIds[i]);
663 }
664 materialMeshes[materialMeshes.size() - 1 - i]->GetCellData()->SetScalars(
665 materialNumberArray);
666
667 // delete all point data, so it is not in ouput
668 // TODO this includes signed distance information which could be
669 // conserved for debugging also includes wheter a cell was vaild for
670 // cutting by the grid
671 vtkSmartPointer<vtkPointData> pointData =
672 materialMeshes[materialMeshes.size() - 1 - i]->GetPointData();
673 const int numberOfArrays = pointData->GetNumberOfArrays();
674 for (int j = 0; j < numberOfArrays; ++j) {
675 pointData->RemoveArray(0); // remove first array until none are left
676 }
677
678 // if hull mesh should be exported, create hull for each layer and put
679 // them together
680 if (extractHullMesh) {
681 vtkSmartPointer<vtkGeometryFilter> geoFilter =
682 vtkSmartPointer<vtkGeometryFilter>::New();
683 geoFilter->SetInputData(materialMeshes[materialMeshes.size() - 1 - i]);
684 geoFilter->Update();
685 hullAppendFilter->AddInputData(geoFilter->GetOutput());
686 }
687
688 appendFilter->AddInputData(materialMeshes[materialMeshes.size() - 1 - i]);
689 }
690
691 // do not need tetrahedral volume mesh if we do not print volume
692 auto volumeVTK = vtkSmartPointer<vtkUnstructuredGrid>::New();
693 auto hullVTK = vtkSmartPointer<vtkPolyData>::New();
694 if (extractVolumeMesh) {
695 appendFilter->Update();
696
697 // remove degenerate points and remove cells which collapse to zero
698 // volume then
699 volumeVTK = appendFilter->GetOutput();
700#ifdef LS_TO_VISUALIZATION_DEBUG
701 {
702 std::cout << "Before duplicate removal: " << std::endl;
703 std::cout << "Points: " << volumeVTK->GetNumberOfPoints() << std::endl;
704 std::cout << "Cells: " << volumeVTK->GetNumberOfCells() << std::endl;
705 vtkSmartPointer<vtkXMLUnstructuredGridWriter> gwriter =
706 vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
707 gwriter->SetFileName("before_removal.vtu");
708 gwriter->SetInputData(appendFilter->GetOutput());
709 gwriter->Update();
710 }
711#endif
712
713 // use 1/1000th of grid spacing for contraction of two similar points,
714 // so that tetrahedralisation works correctly
715 removeDuplicatePoints(volumeVTK, 1e-3 * gridDelta);
716
717#ifdef LS_TO_VISUALIZATION_DEBUG
718 {
719 std::cout << "After duplicate removal: " << std::endl;
720 std::cout << "Points: " << volumeVTK->GetNumberOfPoints() << std::endl;
721 std::cout << "Cells: " << volumeVTK->GetNumberOfCells() << std::endl;
722 vtkSmartPointer<vtkXMLUnstructuredGridWriter> gwriter =
723 vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
724 gwriter->SetFileName("after_removal.vtu");
725 gwriter->SetInputData(volumeVTK);
726 gwriter->Update();
727 }
728#endif
729
730 // change all 3D cells into tetras and all 2D cells to triangles
731 vtkSmartPointer<vtkDataSetTriangleFilter> triangleFilter =
732 vtkSmartPointer<vtkDataSetTriangleFilter>::New();
733 triangleFilter->SetInputData(volumeVTK);
734 triangleFilter->Update();
735 volumeVTK = triangleFilter->GetOutput();
736
737 // now that only tetras are left, remove the ones with degenerate points
738 removeDegenerateTetras(volumeVTK);
739
740 // add meta data to volume mesh
741 addMetaDataToVTK(volumeVTK);
742
743 auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
744 writer->SetFileName((fileName + "_volume.vtu").c_str());
745 writer->SetInputData(volumeVTK);
746 writer->Write();
747 }
748
749 // Now make hull mesh if necessary
750 if (extractHullMesh) {
751 hullAppendFilter->Update();
752 hullVTK = hullAppendFilter->GetOutput();
753 // use 1/1000th of grid spacing for contraction of two similar points
754 removeDuplicatePoints(hullVTK, 1e-3 * gridDelta);
755
756 vtkSmartPointer<vtkTriangleFilter> hullTriangleFilter =
757 vtkSmartPointer<vtkTriangleFilter>::New();
758 hullTriangleFilter->SetInputData(hullVTK);
759 hullTriangleFilter->Update();
760
761 hullVTK = hullTriangleFilter->GetOutput();
762
763 // add meta data to hull mesh
764 addMetaDataToVTK(hullVTK);
765
766 auto writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
767 writer->SetFileName((fileName + "_hull.vtp").c_str());
768 writer->SetInputData(hullVTK);
769 writer->Write();
770 }
771 }
772}; // namespace viennals
773
774// add all template specialisations for this class
775PRECOMPILE_PRECISION_DIMENSION(WriteVisualizationMesh)
776
777} // namespace viennals
778
779#endif // VIENNALS_USE_VTK
constexpr int D
Definition Epitaxy.cpp:9
double T
Definition Epitaxy.cpp:10
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
float gridDelta
Definition AirGapDeposition.py:21
writer
Definition AirGapDeposition.py:89
int counter
Definition Deposition.py:71
Definition lsAdvect.hpp:36